diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index de11f63..04e18c8 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -9,3 +9,4 @@ | JoeCartonKelly-MO | Joseph Carton-Kelly | Met Office | 2026-04-15 | | yg460-cam | Yao Ge | University of Cambridge | 2026-04-17 | | theabro | Nathan Luke Abraham | NCAS & University of Cambridge | 2026-03-19 | +| mn416 | Matthew Naylor | University of Cambridge, 2026-07-06 | diff --git a/CONTRIBUTORS_LIST.txt b/CONTRIBUTORS_LIST.txt index 220f842..5378b7c 100644 --- a/CONTRIBUTORS_LIST.txt +++ b/CONTRIBUTORS_LIST.txt @@ -26,6 +26,7 @@ Florent Malavelle Met Office Andy Malcolm Met Office Cameron McErlich University of Canterbury Jane Mulcahy Met Office +Matthew Naylor University of Cambridge Charlotte Norris Met Office Amy Peace Met Office Leighton Regayre Met Office and University of Leeds diff --git a/src/control/core/interface/ukca_config_specification_mod.F90 b/src/control/core/interface/ukca_config_specification_mod.F90 index a32b743..dc79d5a 100644 --- a/src/control/core/interface/ukca_config_specification_mod.F90 +++ b/src/control/core/interface/ukca_config_specification_mod.F90 @@ -185,6 +185,8 @@ MODULE ukca_config_specification_mod ! (2-50, 3 recommended) INTEGER :: ukca_chem_seg_size ! Grid points per segment for the ! column-based chemical solver + INTEGER :: ukca_chem_full_chunk_size(3) ! Grid points per 3D chunk in the + ! full-domain-based chemical solver REAL :: max_z_for_offline_chem ! Maximum height at which to integrate ! chemistry with the explicit B-E ! Offline Oxidants scheme @@ -874,6 +876,7 @@ SUBROUTINE init_ukca_configuration() ukca_config%i_ukca_quasinewton_start = imdi ukca_config%i_ukca_quasinewton_end = imdi ukca_config%ukca_chem_seg_size = imdi +ukca_config%ukca_chem_full_chunk_size(:) = [-1, -1, -1] ukca_config%max_z_for_offline_chem = rmdi ukca_config%nlev_above_trop_o3_env = imdi ukca_config%nlev_ch4_stratloss = imdi @@ -1090,6 +1093,7 @@ SUBROUTINE ukca_get_config( & dts0, nit, & i_ukca_quasinewton_start, i_ukca_quasinewton_end, & ukca_chem_seg_size, & + ukca_chem_full_chunk_size, & nlev_above_trop_o3_env, & nlev_ch4_stratloss, & i_ukca_topboundary, & @@ -1284,6 +1288,7 @@ SUBROUTINE ukca_get_config( & INTEGER, OPTIONAL, INTENT(OUT) :: i_ukca_quasinewton_start INTEGER, OPTIONAL, INTENT(OUT) :: i_ukca_quasinewton_end INTEGER, OPTIONAL, INTENT(OUT) :: ukca_chem_seg_size +INTEGER, OPTIONAL, INTENT(OUT) :: ukca_chem_full_chunk_size(3) INTEGER, OPTIONAL, INTENT(OUT) :: nlev_above_trop_o3_env INTEGER, OPTIONAL, INTENT(OUT) :: nlev_ch4_stratloss INTEGER, OPTIONAL, INTENT(OUT) :: i_ukca_topboundary @@ -1572,6 +1577,8 @@ SUBROUTINE ukca_get_config( & i_ukca_quasinewton_end = ukca_config%i_ukca_quasinewton_end IF (PRESENT(ukca_chem_seg_size)) & ukca_chem_seg_size = ukca_config%ukca_chem_seg_size +IF (PRESENT(ukca_chem_full_chunk_size)) & + ukca_chem_full_chunk_size = ukca_config%ukca_chem_full_chunk_size IF (PRESENT(max_z_for_offline_chem)) & max_z_for_offline_chem = ukca_config%max_z_for_offline_chem IF (PRESENT(nlev_above_trop_o3_env)) & diff --git a/src/control/core/top_level/ukca_main1-ukca_main1.F90 b/src/control/core/top_level/ukca_main1-ukca_main1.F90 index 044776f..d5ce0df 100644 --- a/src/control/core/top_level/ukca_main1-ukca_main1.F90 +++ b/src/control/core/top_level/ukca_main1-ukca_main1.F90 @@ -581,6 +581,9 @@ SUBROUTINE ukca_main1(error_code_ptr, timestep_number, current_time, & TYPE(autotune_type), ALLOCATABLE, SAVE :: autotune_state #endif +! Full-domain chunk size in each dimension +INTEGER :: full_chunk_x, full_chunk_y, full_chunk_z + CHARACTER(LEN=*), PARAMETER :: RoutineName='UKCA_MAIN1' !- End of header @@ -601,10 +604,19 @@ SUBROUTINE ukca_main1(error_code_ptr, timestep_number, current_time, & rows = ukca_config%rows model_levels = ukca_config%model_levels +! Determine chunk size for full-domain 3D chunking +full_chunk_x = ukca_config%ukca_chem_full_chunk_size(1) +full_chunk_y = ukca_config%ukca_chem_full_chunk_size(2) +full_chunk_z = ukca_config%ukca_chem_full_chunk_size(3) +IF (full_chunk_x <= 0 .OR. full_chunk_x > row_length) full_chunk_x = row_length +IF (full_chunk_y <= 0 .OR. full_chunk_y > rows) full_chunk_y = rows +IF (full_chunk_z <= 0 .OR. full_chunk_z > model_levels) & + full_chunk_z = model_levels + theta_field_size = row_length * rows tot_n_pnts = theta_field_size * model_levels IF (ukca_config%l_ukca_asad_full) THEN - n_pnts = tot_n_pnts + n_pnts = full_chunk_x * full_chunk_y * full_chunk_z ELSE IF (ukca_config%l_ukca_asad_columns) THEN n_pnts = ukca_config%ukca_chem_seg_size ELSE @@ -2340,6 +2352,7 @@ SUBROUTINE ukca_main1(error_code_ptr, timestep_number, current_time, & CALL ukca_chemistry_ctl_full( & row_length, rows, model_levels, & theta_field_size, tot_n_pnts, & + full_chunk_x, full_chunk_y, full_chunk_z, & n_chem_tracers+n_aero_tracers, & istore_h2so4, & p_theta_levels, & diff --git a/src/control/core/top_level/ukca_setup_mod.F90 b/src/control/core/top_level/ukca_setup_mod.F90 index 9a347c1..d678b4a 100644 --- a/src/control/core/top_level/ukca_setup_mod.F90 +++ b/src/control/core/top_level/ukca_setup_mod.F90 @@ -81,6 +81,7 @@ SUBROUTINE ukca_setup(error_code, & i_ukca_quasinewton_start, & i_ukca_quasinewton_end, & ukca_chem_seg_size, & + ukca_chem_full_chunk_size, & nlev_above_trop_o3_env, & nlev_ch4_stratloss, & i_ukca_topboundary, & @@ -414,6 +415,7 @@ SUBROUTINE ukca_setup(error_code, & INTEGER, OPTIONAL, INTENT(IN) :: i_ukca_quasinewton_start INTEGER, OPTIONAL, INTENT(IN) :: i_ukca_quasinewton_end INTEGER, OPTIONAL, INTENT(IN) :: ukca_chem_seg_size +INTEGER, OPTIONAL, INTENT(IN) :: ukca_chem_full_chunk_size(3) INTEGER, OPTIONAL, INTENT(IN) :: nlev_above_trop_o3_env INTEGER, OPTIONAL, INTENT(IN) :: nlev_ch4_stratloss INTEGER, OPTIONAL, INTENT(IN) :: i_ukca_topboundary @@ -838,6 +840,13 @@ SUBROUTINE ukca_setup(error_code, & END IF + ! Full-domain-based run configuration + IF (ukca_config%l_ukca_asad_full) THEN + ukca_config%ukca_chem_full_chunk_size(:) = [-1, -1, -1] + IF (PRESENT(ukca_chem_full_chunk_size)) & + ukca_config%ukca_chem_full_chunk_size(:) = ukca_chem_full_chunk_size(:) + END IF + ! Configuration specific to explicit B-E Offline Oxidants scheme IF (ukca_config%i_ukca_chem == i_ukca_chem_offline_be) THEN IF (PRESENT(max_z_for_offline_chem)) & diff --git a/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 b/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 index bc84636..a1ee2a6 100644 --- a/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 +++ b/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 @@ -37,6 +37,7 @@ MODULE ukca_chemistry_ctl_full_mod SUBROUTINE ukca_chemistry_ctl_full( & row_length, rows, model_levels, theta_field_size, tot_n_pnts, & + chunk_n_x, chunk_n_y, chunk_n_z, & ntracers, & istore_h2so4, & pres, temp, q, & @@ -106,60 +107,64 @@ SUBROUTINE ukca_chemistry_ctl_full( & INTEGER, INTENT(IN) :: model_levels ! size of UKCA z dimension INTEGER, INTENT(IN) :: theta_field_size ! no. of points in horizontal INTEGER, INTENT(IN) :: tot_n_pnts ! no. of points in full domain +INTEGER, INTENT(IN) :: chunk_n_x ! chunk size in x dimension +INTEGER, INTENT(IN) :: chunk_n_y ! chunk size in y dimension +INTEGER, INTENT(IN) :: chunk_n_z ! chunk size in z dimension INTEGER, INTENT(IN) :: ntracers ! no. of tracers INTEGER, INTENT(IN) :: uph2so4inaer ! flag for H2SO4 updating INTEGER, INTENT(IN) :: istore_h2so4 ! location of H2SO4 in f array -INTEGER, INTENT(IN) :: nlev_with_ddep(theta_field_size) ! No levs in bl - -REAL, INTENT(IN) :: pres(tot_n_pnts) ! pressure -REAL, INTENT(IN) :: temp(tot_n_pnts) ! actual temperature -REAL, INTENT(IN) :: volume(tot_n_pnts) ! cell volume -REAL, INTENT(IN) :: H_plus(tot_n_pnts) ! pH array -REAL, INTENT(IN) :: qcf(tot_n_pnts) -REAL, INTENT(IN) :: qcl(tot_n_pnts) -REAL, INTENT(IN) :: cloud_frac(tot_n_pnts) -REAL, INTENT(IN) :: so4_sa(tot_n_pnts) ! aerosol surface area -REAL, INTENT(IN) :: zdryrt(theta_field_size,jpdd) ! dry dep rate -REAL, INTENT(IN) :: zwetrt(tot_n_pnts,jpdw) ! wet dep rate -REAL, INTENT(IN) :: photol_rates(tot_n_pnts,jppj) -REAL, INTENT(IN) :: co2_interactive(tot_n_pnts) -REAL, INTENT(OUT) :: shno3(tot_n_pnts) -REAL, INTENT(IN OUT) :: q(tot_n_pnts) ! water vapour -REAL, INTENT(IN OUT) :: tracer(tot_n_pnts,ntracers) ! tracer MMR +INTEGER, INTENT(IN) :: nlev_with_ddep(row_length,rows) ! No levs in bl + +REAL, INTENT(IN) :: pres(row_length,rows,model_levels) ! pressure +REAL, INTENT(IN) :: temp(row_length,rows,model_levels) ! actual temperature +REAL, INTENT(IN) :: volume(row_length,rows,model_levels) ! cell volume +REAL, INTENT(IN) :: H_plus(row_length,rows,model_levels) ! pH array +REAL, INTENT(IN) :: qcf(row_length,rows,model_levels) +REAL, INTENT(IN) :: qcl(row_length,rows,model_levels) +REAL, INTENT(IN) :: cloud_frac(row_length,rows,model_levels) +REAL, INTENT(IN) :: so4_sa(row_length,rows,model_levels) ! aerosol surface area +REAL, INTENT(IN) :: zdryrt(row_length,rows,jpdd) ! dry dep rate +REAL, INTENT(IN) :: zwetrt(row_length,rows,model_levels,jpdw) ! wet dep rate +REAL, INTENT(IN) :: photol_rates(row_length,rows,model_levels,jppj) +REAL, INTENT(IN) :: co2_interactive(row_length,rows,model_levels) +REAL, INTENT(OUT) :: shno3(row_length,rows,model_levels) +REAL, INTENT(IN OUT) :: q(row_length,rows,model_levels) ! water vapour +REAL, INTENT(IN OUT) :: tracer(row_length,rows,model_levels, & + ntracers) ! tracer MMR ! SO2 increments -REAL, INTENT(IN OUT) :: delSO2_wet_H2O2(tot_n_pnts) -REAL, INTENT(IN OUT) :: delSO2_wet_O3(tot_n_pnts) -REAL, INTENT(IN OUT) :: delh2so4_chem(tot_n_pnts) +REAL, INTENT(IN OUT) :: delSO2_wet_H2O2(row_length,rows,model_levels) +REAL, INTENT(IN OUT) :: delSO2_wet_O3(row_length,rows,model_levels) +REAL, INTENT(IN OUT) :: delh2so4_chem(row_length,rows,model_levels) ! Atmospheric Burden of CH4 -REAL, INTENT(IN OUT) :: atm_ch4_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_ch4_mol(row_length,rows,model_levels) ! Atmospheric Burden of CO -REAL, INTENT(IN OUT) :: atm_co_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_co_mol(row_length,rows,model_levels) ! Atmospheric Burden of Nitrous Oxide (N2O) -REAL, INTENT(IN OUT) :: atm_n2o_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_n2o_mol(row_length,rows,model_levels) ! Atmospheric Burden of CFC-12 -REAL, INTENT(IN OUT) :: atm_cf2cl2_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_cf2cl2_mol(row_length,rows,model_levels) ! Atmospheric Burden of CFC-11 -REAL, INTENT(IN OUT) :: atm_cfcl3_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_cfcl3_mol(row_length,rows,model_levels) ! Atmospheric Burden of CH3Br -REAL, INTENT(IN OUT) :: atm_mebr_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_mebr_mol(row_length,rows,model_levels) ! Atmospheric Burden of H2 -REAL, INTENT(IN OUT) :: atm_h2_mol(tot_n_pnts) +REAL, INTENT(IN OUT) :: atm_h2_mol(row_length,rows,model_levels) ! Non transported prognostics TYPE(ntp_type), INTENT(IN OUT) :: all_ntp(dim_ntp) ! Mask to limit formation of Nat below specified height -LOGICAL, INTENT(IN) :: have_nat(tot_n_pnts) +LOGICAL, INTENT(IN) :: have_nat(row_length,rows,model_levels) ! Stratospheric mask -LOGICAL, INTENT(IN) :: L_stratosphere(tot_n_pnts) +LOGICAL, INTENT(IN) :: L_stratosphere(row_length,rows,model_levels) ! Flag for determining if this is the first chemistry call LOGICAL, INTENT(IN) :: firstcall @@ -184,27 +189,61 @@ SUBROUTINE ukca_chemistry_ctl_full( & CHARACTER(LEN=errormessagelength) :: cmessage ! Error message CHARACTER(LEN=10) :: prods(2) ! Products CHARACTER(LEN=10) :: prods3(3) ! Products -LOGICAL :: ddmask(theta_field_size) ! mask +LOGICAL :: ddmask(row_length,rows) ! mask -REAL, ALLOCATABLE :: ystore(:) ! array for H2SO4 when updated in MODE -REAL :: zftr(tot_n_pnts,jpcspf) ! 1-D array of chemically active species +REAL, ALLOCATABLE :: ystore(:,:,:) ! array for H2SO4 when updated in MODE +REAL :: zftr(row_length,rows,model_levels,jpcspf) + ! array of chemically active species ! including RO2 species, in VMR -REAL :: zq(tot_n_pnts) ! 1-D water vapour vmr -REAL :: co2_1d(tot_n_pnts) ! 1-D CO2 -REAL :: zprt1d(tot_n_pnts,jppj) ! 1-D photolysis rates for ASAD -REAL :: zdryrt2(tot_n_pnts,jpdd) ! dry dep rate -REAL :: rc_het(tot_n_pnts,2) ! heterog rates for trop chem +REAL :: zq(row_length,rows,model_levels) ! water vapour vmr +REAL :: co2_tmp(row_length,rows,model_levels) ! CO2 +REAL :: zprt(row_length,rows,model_levels,jppj) ! photolysis rates for ASAD +REAL :: zdryrt2(row_length,rows,model_levels,jpdd) ! dry dep rate +REAL :: rc_het(row_length,rows,model_levels,2) ! heterog rates for trop chem ! 1-D mask for stratosphere -LOGICAL :: stratflag(tot_n_pnts) +LOGICAL :: stratflag(row_length,rows,model_levels) ! Full ntp array -REAL :: ntp_data(tot_n_pnts,dim_ntp) +REAL :: ntp_data(row_length,rows,model_levels,dim_ntp) INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle +! Full-domain versions of ASAD modules variables +REAL :: full_dpd(row_length,rows,model_levels,jpspec) +REAL :: full_dpw(row_length,rows,model_levels,jpspec) +REAL :: full_fpsc1(row_length,rows,model_levels) +REAL :: full_fpsc2(row_length,rows,model_levels) +REAL :: full_prk(row_length,rows,model_levels,jpnr) +REAL :: full_y(row_length,rows,model_levels,jpspec) +REAL :: full_sph2o(row_length,rows,model_levels) +REAL :: full_za(row_length,rows,model_levels) +REAL :: full_rk(row_length,rows,model_levels,jpnr) +REAL :: full_tnd(row_length,rows,model_levels) +REAL :: full_sphno3(row_length,rows,model_levels) +REAL :: full_zftr(row_length,rows,model_levels,jpcspf) + +! Chunked versions of ASAD modules variables +REAL :: chunk_dpd(chunk_n_x,chunk_n_y,chunk_n_z,jpspec) +REAL :: chunk_dpw(chunk_n_x,chunk_n_y,chunk_n_z,jpspec) +REAL :: chunk_fpsc1(chunk_n_x,chunk_n_y,chunk_n_z) +REAL :: chunk_fpsc2(chunk_n_x,chunk_n_y,chunk_n_z) +REAL :: chunk_prk(chunk_n_x,chunk_n_y,chunk_n_z,jpnr) +REAL :: chunk_y(chunk_n_x,chunk_n_y,chunk_n_z,jpspec) +REAL :: chunk_rk(chunk_n_x,chunk_n_y,chunk_n_z,jpnr) +REAL :: chunk_tnd(chunk_n_x,chunk_n_y,chunk_n_z) +REAL :: chunk_sphno3(chunk_n_x,chunk_n_y,chunk_n_z) +REAL :: chunk_zftr(chunk_n_x,chunk_n_y,chunk_n_z,jpcspf) + +! Chunk start, end, iterator, and offset in each dimension +INTEGER :: xs, xe, xi, xo +INTEGER :: ys, ye, yi, yo +INTEGER :: zs, ze, zi, zo +INTEGER :: chunk_n_pnts +LOGICAL :: chunking_enabled + CHARACTER(LEN=*), PARAMETER :: RoutineName='UKCA_CHEMISTRY_CTL_FULL' @@ -222,35 +261,35 @@ SUBROUTINE ukca_chemistry_ctl_full( & END IF ! Fill stratospheric flag indicator and ntp_data array -stratflag(:) = (L_stratosphere(:)) +stratflag(:,:,:) = L_stratosphere(:,:,:) !$OMP PARALLEL DEFAULT(NONE) PRIVATE(l) & !$OMP SHARED(all_ntp, ntp_data, tot_n_pnts) !$OMP DO SCHEDULE(DYNAMIC) DO l = 1, dim_ntp IF (all_ntp(l)%l_required) THEN - ntp_data(:,l) = RESHAPE(all_ntp(l)%data_3d(:,:,:),[tot_n_pnts]) + ntp_data(:,:,:,l) = all_ntp(l)%data_3d(:,:,:) END IF END DO !$OMP END DO !$OMP END PARALLEL -! Put tracer mmr into 1-D array for use in ASAD chemical solver -zq(:) = q(:)/c_h2o +! Store tracer mmr into for use in ASAD chemical solver +zq(:,:,:) = q(:,:,:)/c_h2o -! Map photolysis rates onto 1-D array. +! Store photolysis rates IF (ukca_config%l_ukca_offline) THEN ! Offline chemistry has no photolysis - zprt1d(:,:) = 0.0 + zprt(:,:,:,:) = 0.0 ELSE - zprt1d(:,:) = photol_rates(:,:) + zprt(:,:,:,:) = photol_rates(:,:,:,:) END IF IF (ANY(speci(:) == 'CO2 ')) THEN ! Copy the CO2 concentration into the asad module as VMR IF (ukca_config%l_chem_environ_co2_fld) THEN - co2_1d(:) = co2_interactive(:)/c_co2 + co2_tmp(:,:,:) = co2_interactive(:,:,:)/c_co2 ELSE - co2_1d(:) = rmdi + co2_tmp(:,:,:) = rmdi END IF END IF ! CO2 as species @@ -258,29 +297,27 @@ SUBROUTINE ukca_chemistry_ctl_full( & IF (ukca_config%l_ukca_trophet) THEN ! N2O5 l = name2ntpindex('het_n2o5 ') - rc_het(:,1) = ntp_data(:,l) + rc_het(:,:,:,1) = ntp_data(:,:,:,l) ! HO2+HO2 l = name2ntpindex('het_ho2 ') - rc_het(:,2) = ntp_data(:,l) + rc_het(:,:,:,2) = ntp_data(:,:,:,l) ELSE - rc_het(:,:) = 0.0 + rc_het(:,:,:,:) = 0.0 END IF -zdryrt2(:,:) = 0.0e0 +zdryrt2(:,:,:,:) = 0.0 IF (ukca_config%l_ukca_intdd) THEN ! Interactive scheme extracts from levels in boundary layer DO k = 1,model_levels - kcs = (k - 1) * theta_field_size + 1 - kce = k * theta_field_size - ddmask(:) = (k <= nlev_with_ddep(:)) DO l=1,jpdd - WHERE (ddmask(:)) - zdryrt2(kcs:kce,l) = zdryrt(:,l) + ddmask(:,:) = k <= nlev_with_ddep(:,:) + WHERE (ddmask(:,:)) + zdryrt2(:,:,k,l) = zdryrt(:,:,l) END WHERE END DO END DO ELSE ! non-interactive - zdryrt2(1:theta_field_size,:) = zdryrt(:,:) + zdryrt2(:,:,1,:) = zdryrt(:,:,:) END IF ! Convert mmr into vmr for tracers. @@ -296,7 +333,7 @@ SUBROUTINE ukca_chemistry_ctl_full( & IF (advt(jtr) == speci(js)) THEN jspf = jspf+1 ! Map data from tracer 3D array into zftr - zftr(:,jspf) = tracer(:,jtr)/c_species(jtr) + full_zftr(:,:,:,jspf) = tracer(:,:,:,jtr)/c_species(jtr) END IF END DO ! Close advected tracer do loop @@ -314,7 +351,7 @@ SUBROUTINE ukca_chemistry_ctl_full( & ! Map data from all_ntp 3D array into zftr IF (nadvt(jna) == spro2(jro2)) THEN jspf = jspf+1 - zftr(:,jspf) = ntp_data(:,l) / c_na_species(jna) + full_zftr(:,:,:,jspf) = ntp_data(:,:,:,l) / c_na_species(jna) ELSE errcode = jro2 WRITE(umMessage,'(A)') '** ERROR in ukca_chemistry_ctl_full' @@ -348,69 +385,167 @@ SUBROUTINE ukca_chemistry_ctl_full( & END IF IF (.NOT. ALLOCATED(ystore) .AND. uph2so4inaer == 1) & - ALLOCATE(ystore(tot_n_pnts)) + ALLOCATE(ystore(row_length,rows,model_levels)) ! Copy water vapour and ice field into 1-D arrays +full_sph2o(:,:,:) = 0.0 IF (ukca_config%l_ukca_het_psc) THEN - sph2o(:) = qcf(:)/c_h2o + full_sph2o(:,:,:) = qcf(:,:,:) / c_h2o END IF ! Fill SO4 surface area -za(:) = so4_sa(:) +full_za(:,:,:) = so4_sa(:,:,:) + +! Number of grid points in every chunk +chunk_n_pnts = chunk_n_x * chunk_n_y * chunk_n_z +chunking_enabled = chunk_n_pnts /= tot_n_pnts IF (uph2so4inaer == 1) THEN ! H2SO4 will be updated in MODE, so store old value here IF (ukca_config%l_fix_ukca_h2so4_ystore) THEN ! primary array passed is zftr, so save this, NOT y - ystore(:) = zftr(:,istore_h2so4) - ELSE - ystore(:) = y(:,nn_h2so4) + ystore(:,:,:) = full_zftr(:,:,:,istore_h2so4) + ELSE IF (.NOT. chunking_enabled) THEN + ! Preserve non-fixed behaviour when chunking is disabled + ystore(:,:,:) = RESHAPE(y(:,nn_h2so4), [row_length,rows,model_levels]) END IF END IF -! Call ASAD routines to do chemistry integration. -! -! Note: unlike in the other chemistry_ctl variants, here we call asad_cdrive -! once, passing in all data points from the current MPI process. - -CALL asad_cdrive(zftr, pres, temp, zq, co2_1d, cloud_frac, qcl, ix, jy, k, & - zdryrt2, zwetrt, rc_het, zprt1d, tot_n_pnts, have_nat, & - stratflag, H_plus) +! This next section of code calls ASAD routines to do chemistry integration +! with 3D chunking. Unlike column-mode chunking, the chunk size is +! parameterised in all three dimensions (for maximum flexibility) and +! we do not resize ASAD arrays when the chunk size does not divide evenly +! into the full-domain size. Instead, we retract the starting point of the +! final chunk in any dimension to obtain a full-chunk worth of data. This +! leads to some redundant computation (some grid points form part of two +! different chunks) but, for small chunk sizes, this has a neglible impact +! on performance. As different chunks can be processed in parallel, care is +! taken to avoid data races where two different threads are computing the +! same grid point at the same time. + +! 3D chunking loop +!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(DYNAMIC) COLLAPSE(3) & +!$OMP PRIVATE(zi, ze, zs, zo, yi, ye, ys, yo, xi, xe, xs, xo, & +!$OMP chunk_zftr, chunk_dpd, chunk_dpw, chunk_fpsc1, chunk_fpsc2, & +!$OMP chunk_prk, chunk_y, chunk_sphno3, chunk_tnd, chunk_rk) & +!$OMP SHARED(chunk_n_pnts, chunk_n_x, chunk_n_y, chunk_n_z, model_levels, & +!$OMP rows, row_length, jpspec, full_sph2o, full_za, full_y, full_dpd, & +!$OMP full_dpw, full_fpsc1, full_fpsc2, full_prk, full_sphno3, & +!$OMP full_tnd, full_rk, pres, temp, zq, co2_tmp, cloud_frac, qcl, & +!$OMP zdryrt2, zwetrt, zprt, rc_het, have_nat, stratflag, H_plus, & +!$OMP zftr, full_zftr, ix, jy, k) +DO zi = 1, model_levels, chunk_n_z + DO yi = 1, rows, chunk_n_y + DO xi = 1, row_length, chunk_n_x + ! Chunk start, end, and offset in Z dimension + ze = MIN(model_levels, zi + (chunk_n_z - 1)) + zs = ze - chunk_n_z + 1 + zo = zi - zs + 1 + + ! Chunk start, end, and offset in Y dimension + ye = MIN(rows, yi + (chunk_n_y - 1)) + ys = ye - chunk_n_y + 1 + yo = yi - ys + 1 + + ! Chunk start, end, and offset in X dimension + xe = MIN(row_length, xi + (chunk_n_x - 1)) + xs = xe - chunk_n_x + 1 + xo = xi - xs + 1 + + ! Copy from full-domain arrays to ASAD module variables + sph2o(:) = RESHAPE(full_sph2o(xs:xe,ys:ye,zs:ze), SHAPE(sph2o)) + za(:) = RESHAPE(full_za(xs:xe,ys:ye,zs:ze), SHAPE(za)) + chunk_zftr(:,:,:,:) = full_zftr(xs:xe,ys:ye,zs:ze,:) + + CALL asad_cdrive( & + chunk_zftr, & + pres(xs:xe,ys:ye,zs:ze), & + temp(xs:xe,ys:ye,zs:ze), & + zq(xs:xe,ys:ye,zs:ze), & + co2_tmp(xs:xe,ys:ye,zs:ze), & + cloud_frac(xs:xe,ys:ye,zs:ze), & + qcl(xs:xe,ys:ye,zs:ze), & + ix, & + jy, & + k, & + zdryrt2(xs:xe,ys:ye,zs:ze,:), & + zwetrt(xs:xe,ys:ye,zs:ze,:), & + rc_het(xs:xe,ys:ye,zs:ze,:), & + zprt(xs:xe,ys:ye,zs:ze,:), & + chunk_n_pnts, & + have_nat(xs:xe,ys:ye,zs:ze), & + stratflag(xs:xe,ys:ye,zs:ze), & + H_plus(xs:xe,ys:ye,zs:ze)) + + ! Copy from ASAD module variables to full-domain arrays + zftr(xi:xe,yi:ye,zi:ze,:) = chunk_zftr(xo:,yo:,zo:,:) + + chunk_dpd(:,:,:,:) = RESHAPE(dpd, SHAPE(chunk_dpd)) + full_dpd(xi:xe,yi:ye,zi:ze,:) = chunk_dpd(xo:,yo:,zo:,:) + + chunk_dpw(:,:,:,:) = RESHAPE(dpw, SHAPE(chunk_dpw)) + full_dpw(xi:xe,yi:ye,zi:ze,:) = chunk_dpw(xo:,yo:,zo:,:) + + chunk_fpsc1(:,:,:) = RESHAPE(fpsc1, SHAPE(chunk_fpsc1)) + full_fpsc1(xi:xe,yi:ye,zi:ze) = chunk_fpsc1(xo:,yo:,zo:) + + chunk_fpsc2(:,:,:) = RESHAPE(fpsc2, SHAPE(chunk_fpsc2)) + full_fpsc2(xi:xe,yi:ye,zi:ze) = chunk_fpsc2(xo:,yo:,zo:) + + chunk_prk(:,:,:,:) = RESHAPE(prk, SHAPE(chunk_prk)) + full_prk(xi:xe,yi:ye,zi:ze,:) = chunk_prk(xo:,yo:,zo:,:) + + chunk_y(:,:,:,:) = RESHAPE(y, SHAPE(chunk_y)) + full_y(xi:xe,yi:ye,zi:ze,:) = chunk_y(xo:,yo:,zo:,:) + + chunk_sphno3(:,:,:) = RESHAPE(sphno3, SHAPE(chunk_sphno3)) + full_sphno3(xi:xe,yi:ye,zi:ze) = chunk_sphno3(xo:,yo:,zo:) + + chunk_tnd(:,:,:) = RESHAPE(tnd, SHAPE(chunk_tnd)) + full_tnd(xi:xe,yi:ye,zi:ze) = chunk_tnd(xo:,yo:,zo:) + + chunk_rk(:,:,:,:) = RESHAPE(rk, SHAPE(chunk_rk)) + full_rk(xi:xe,yi:ye,zi:ze,:) = chunk_rk(xo:,yo:,zo:,:) + END DO + END DO +END DO +!$OMP END PARALLEL DO IF (ukca_config%l_ukca_het_psc) THEN ! Save MMR of NAT PSC particles into 3-D array for PSC sedimentation. ! Note that sphno3 is NAT in number density of HNO3. - IF (ANY(sphno3(:) > 0.0)) THEN - shno3(:) = sphno3(:)/tnd(:)*c_hono2 + IF (ANY(full_sphno3(:,:,:) > 0.0)) THEN + shno3(:,:,:) = full_sphno3(:,:,:)/full_tnd(:,:,:)*c_hono2 ELSE - shno3(:) = 0.0 + shno3(:,:,:) = 0.0 END IF END IF IF (ukca_config%l_ukca_chem .AND. ukca_config%l_ukca_nr_aqchem) THEN ! Calculate chemical fluxes for MODE IF (ihso3_h2o2 > 0) THEN - delSO2_wet_H2O2(:) = delSO2_wet_H2O2(:) + & - rk(:,ihso3_h2o2)*y(:,nn_so2)*y(:,nn_h2o2)*cdt_diag + delSO2_wet_H2O2(:,:,:) = delSO2_wet_H2O2(:,:,:) + & + full_rk(:,:,:,ihso3_h2o2) * full_y(:,:,:,nn_so2) * & + full_y(:,:,:,nn_h2o2) * cdt_diag END IF IF (ihso3_o3 > 0) THEN - delSO2_wet_O3(:) = delSO2_wet_O3(:) + & - rk(:,ihso3_o3)*y(:,nn_so2)*y(:,nn_o3)*cdt_diag + delSO2_wet_O3(:,:,:) = delSO2_wet_O3(:,:,:) + & + full_rk(:,:,:,ihso3_o3)*full_y(:,:,:,nn_so2)*full_y(:,:,:,nn_o3)*cdt_diag END IF IF (iso3_o3 > 0) THEN - delSO2_wet_O3(:) = delSO2_wet_O3(:) + & - rk(:,iso3_o3)*y(:,nn_so2)*y(:,nn_o3)*cdt_diag + delSO2_wet_O3(:,:,:) = delSO2_wet_O3(:,:,:) + & + full_rk(:,:,:,iso3_o3)*full_y(:,:,:,nn_so2)*full_y(:,:,:,nn_o3)*cdt_diag END IF ! net H2SO4 production - note that this is affected by ! l_fix_ukca_h2so4_ystore above. Y value is concentration ! from chemistry prior to zftr being over-written below IF (iso2_oh > 0 .AND. ih2so4_hv > 0) THEN - delh2so4_chem(:) = delh2so4_chem(:) + & - (rk(:,iso2_oh)*y(:,nn_so2)*y(:,nn_oh) - rk(:,ih2so4_hv)*y(:,nn_h2so4))* & - cdt_diag + delh2so4_chem(:,:,:) = delh2so4_chem(:,:,:) + & + (full_rk(:,:,:,iso2_oh) * full_y(:,:,:,nn_so2) * full_y(:,:,:,nn_oh) - & + full_rk(:,:,:,ih2so4_hv) * full_y(:,:,:,nn_h2so4)) * cdt_diag ELSE IF (iso2_oh > 0) THEN - delh2so4_chem(:) = delh2so4_chem(:) + & - rk(:,iso2_oh)*y(:,nn_so2)*y(:,nn_oh)*cdt_diag + delh2so4_chem(:,:,:) = delh2so4_chem(:,:,:) + full_rk(:,:,:,iso2_oh) * & + full_y(:,:,:,nn_so2) * full_y(:,:,:,nn_oh) * cdt_diag END IF IF (uph2so4inaer == 1) THEN @@ -420,11 +555,12 @@ SUBROUTINE ukca_chemistry_ctl_full( & ! calculate delh2so4_chem as the difference in H2SO4 over chemistry ! zftr is already in VMR, so divide by diagnostic chemistry timestep to ! give as vmr/s - delh2so4_chem(:) = (zftr(:,istore_h2so4) - ystore(:)) / cdt_diag + delh2so4_chem(:,:,:) = (zftr(:,:,:,istore_h2so4) - ystore(:,:,:)) / & + cdt_diag ! primary array passed is zftr, so copy back to this, NOT y - zftr(:,istore_h2so4) = ystore(:) + zftr(:,:,:,istore_h2so4) = ystore(:,:,:) ELSE - y(:,nn_h2so4) = ystore(:) + full_y(:,:,:,nn_h2so4) = ystore(:,:,:) END IF END IF END IF @@ -434,13 +570,13 @@ SUBROUTINE ukca_chemistry_ctl_full( & ((L_asad_use_flux_rxns .OR. L_asad_use_rxn_rates) .OR. & (L_asad_use_wetdep .OR. L_asad_use_drydep))) THEN CALL asad_chemical_diagnostics(row_length,rows,model_levels,tot_n_pnts, & - dpd,dpw,prk,y,jy,ix,klevel,volume,ierr) + full_dpd,full_dpw,full_prk,full_y,jy,ix,klevel,volume,ierr) END IF ! PSC diagnostics IF (L_asad_use_chem_diags .AND. L_asad_use_psc_diagnostic) THEN CALL asad_psc_diagnostic(row_length,rows,model_levels,tot_n_pnts, & - fpsc1,fpsc2,jy,ix,klevel,ierr) + full_fpsc1,full_fpsc2,jy,ix,klevel,ierr) END IF ! Set SS species concentrations for output (stratospheric configurations) @@ -448,13 +584,13 @@ SUBROUTINE ukca_chemistry_ctl_full( & ! O1D mmr IF (O1D_in_ss) THEN l = name2ntpindex('O(1D) ') - ntp_data(:,l) = y(:,nn_o1d)/tnd(:)*c_o1d + ntp_data(:,:,:,l) = full_y(:,:,:,nn_o1d)/full_tnd(:,:,:)*c_o1d END IF ! O3P mmr IF (O3P_in_ss) THEN l = name2ntpindex('O(3P) ') - ntp_data(:,l) = y(:,nn_o3p)/tnd(:)*c_o3p + ntp_data(:,:,:,l) = full_y(:,:,:,nn_o3p)/full_tnd(:,:,:)*c_o3p END IF ! First copy the concentrations from the zftr array to the @@ -466,55 +602,62 @@ SUBROUTINE ukca_chemistry_ctl_full( & !$OMP SHARED(advt, atm_cf2cl2_mol, atm_cfcl3_mol, atm_ch4_mol, atm_co_mol, & !$OMP atm_h2_mol, atm_mebr_mol, atm_n2o_mol, avogadro, jpcspf, & !$OMP n_cf2cl2, n_cfcl3, n_ch4, n_co, n_h2, n_mebr, n_n2o, specf, & -!$OMP volume, zftr) +!$OMP volume, zftr, full_tnd) !$OMP DO SCHEDULE(DYNAMIC) DO jspf = 1, jpcspf IF (n_ch4 > 0) THEN IF (specf(jspf) == advt(n_ch4)) THEN - atm_ch4_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_ch4_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF ! CO IF (n_co > 0) THEN IF (specf(jspf) == advt(n_co)) THEN - atm_co_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_co_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF ! N2O IF (n_n2o > 0) THEN IF (specf(jspf) == advt(n_n2o)) THEN - atm_n2o_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_n2o_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF ! CFC-12 IF (n_cf2cl2 > 0) THEN IF (specf(jspf) == advt(n_cf2cl2)) THEN - atm_cf2cl2_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_cf2cl2_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF ! CFC-11 IF (n_cfcl3 > 0) THEN IF (specf(jspf) == advt(n_cfcl3)) THEN - atm_cfcl3_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_cfcl3_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF ! CH3Br IF (n_mebr > 0) THEN IF (specf(jspf) == advt(n_mebr)) THEN - atm_mebr_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_mebr_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF ! H2 IF (n_h2 > 0) THEN IF (specf(jspf) == advt(n_h2)) THEN - atm_h2_mol(:) = zftr(:,jspf)*tnd(:)*volume(:)*1.0e6/avogadro + atm_h2_mol(:,:,:) = zftr(:,:,:,jspf) * full_tnd(:,:,:) * & + volume(:,:,:) * 1.0e6/avogadro END IF END IF @@ -531,7 +674,7 @@ SUBROUTINE ukca_chemistry_ctl_full( & DO jtr=1,jpctr IF (advt(jtr) == speci(js)) THEN jspf = jspf+1 - tracer(:,jtr) = zftr(:,jspf) * c_species(jtr) + tracer(:,:,:,jtr) = zftr(:,:,:,jspf) * c_species(jtr) END IF END DO @@ -547,7 +690,7 @@ SUBROUTINE ukca_chemistry_ctl_full( & ! Find location of species in nadvt array jna = nlnaro2(jro2) - ntp_data(:,l) = zftr(:,jspf) * c_na_species(jna) + ntp_data(:,:,:,l) = zftr(:,:,:,jspf) * c_na_species(jna) END IF ! Close IF RO2 species END DO ! Close loop through RO2 species END IF ! Close IF RO2_NTP @@ -559,13 +702,18 @@ SUBROUTINE ukca_chemistry_ctl_full( & !$OMP DO SCHEDULE(DYNAMIC) DO l = 1, dim_ntp IF (all_ntp(l)%l_required) THEN - all_ntp(l)%data_3d(:,:,:) = RESHAPE(ntp_data(:,l), & - [row_length,rows,model_levels]) + all_ntp(l)%data_3d(:,:,:) = ntp_data(:,:,:,l) END IF END DO !$OMP END DO !$OMP END PARALLEL +! Preserve non-fixed behaviour when chunking is disabled +IF (.NOT. ukca_config%l_fix_ukca_h2so4_ystore .AND. & + .NOT. chunking_enabled) THEN + y(:,:) = RESHAPE(full_y, [tot_n_pnts,jpspec]) +END IF + IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE ukca_chemistry_ctl_full