From 170b4beedd5b2e421bba8889e88eecfc9cd007f0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 23 Dec 2020 07:42:20 -0700 Subject: [PATCH 01/22] Fixes inconsistency in the calculation of tendency * Also reducing line lenghts > 120 --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 570b4b9ad8..ccb42ddecf 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -222,8 +222,9 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (G%mask2dT(i,j)>0.) then tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) + if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then - tendency(i,j,k) = (tracer%t(i,j,k)-tracer_old(i,j,k)) * Idt + tendency(i,j,k) = (tracer%t(i,j,k)-tracer_old(i,j,k)) * (h(i,j,k) + GV%H_subroundoff) * Idt endif endif enddo ; enddo ; enddo @@ -579,13 +580,14 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ type(lbd_CS), pointer :: CS !< Lateral diffusion control structure !! the boundary layer ! Local variables - real, dimension(:), allocatable :: dz_top !< The LBD z grid to be created [L ~ m] - real, dimension(:), allocatable :: phi_L_z !< Tracer values in the ztop grid (left) [conc] - real, dimension(:), allocatable :: phi_R_z !< Tracer values in the ztop grid (right) [conc] - real, dimension(:), allocatable :: F_layer_z !< Diffusive flux at U- or V-point in the ztop grid [H L2 conc ~> m3 conc] + real, dimension(:), allocatable :: dz_top !< The LBD z grid to be created [L ~ m] + real, dimension(:), allocatable :: phi_L_z !< Tracer values in the ztop grid (left) [conc] + real, dimension(:), allocatable :: phi_R_z !< Tracer values in the ztop grid (right) [conc] + real, dimension(:), allocatable :: F_layer_z !< Diffusive flux at U- or V-point in the ztop grid + !! [H L2 conc ~> m3 conc] real, dimension(ke) :: h_vel !< Thicknesses at u- and v-points in the native grid - !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2] - real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. real :: htot !< Total column thickness [H ~> m or kg m-2] From 99dcfbb1d691c2f644b127c49252bc7f8307e14b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 28 Dec 2020 12:19:40 -0700 Subject: [PATCH 02/22] Calculate tendency using fluxes --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index ccb42ddecf..e7ece53c0a 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -224,7 +224,8 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then - tendency(i,j,k) = (tracer%t(i,j,k)-tracer_old(i,j,k)) * (h(i,j,k) + GV%H_subroundoff) * Idt + tendency(i,j,k) = ((uFlx(I-1,j,k)-uFlx(I,j,k)) + (vFlx(i,J-1,k)-vFlx(i,J,k))) * & + G%IareaT(i,j) * Idt endif endif enddo ; enddo ; enddo From 280e98b7a970b233b7a6bdebc84bbed75ca1451a Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Mon, 4 Jan 2021 12:27:27 -0700 Subject: [PATCH 03/22] Added a counter-based PRNG to MOM_random See https://arxiv.org/abs/2004.06278. Not an exact reproduction of "squares" because Fortran doesn't have a uint64 type, and not all compilers provide integers with > 64 bits... --- src/framework/MOM_random.F90 | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index 161236572c..97e4a1ace8 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -17,6 +17,7 @@ module MOM_random public :: random_0d_constructor public :: random_01 +public :: random_01_CB public :: random_norm public :: random_2d_constructor public :: random_2d_01 @@ -46,6 +47,28 @@ real function random_01(CS) end function random_01 +!> Returns a random number between 0 and 1 +!! See https://arxiv.org/abs/2004.06278. Not an exact reproduction of "squares" because Fortran +!! doesn't have a uint64 type, and not all compilers provide integers with > 64 bits... +real function random_01_CB(ctr, key) + integer, parameter :: int64 = selected_int_kind(10) ! Integer with >= 64 bits + integer, intent(in) :: ctr, key ! counter & key inputs, standard integer kind + integer(kind=int64) :: x, y, z ! Follows "Squares" naming convention + + x = (ctr + 1) * (key + 65536) ! 65536 added because keys below that don't work. + y = (ctr + 1) * (key + 65536) + z = y + (key + 65536) + x = x*x + y + x = ior(ishft(x,32),ishft(x,-32)) + x = x*x + z + x = ior(ishft(x,32),ishft(x,-32)) + x = x*x + y + x = ior(ishft(x,32),ishft(x,-32)) + x = x*x + z + random_01_CB = .5*(1. + .5*real(int(ishft(x,-32)))/real(2**30)) + +end function + !> Returns an approximately normally distributed random number with mean 0 and variance 1 real function random_norm(CS) type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators From 62138bebf91c0df1ad535af7dedc3be1e3f56210 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Tue, 5 Jan 2021 16:03:26 -0700 Subject: [PATCH 04/22] Used int64 from F2003 and added doc strings --- src/framework/MOM_random.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index 97e4a1ace8..95775f9050 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -51,8 +51,9 @@ end function random_01 !! See https://arxiv.org/abs/2004.06278. Not an exact reproduction of "squares" because Fortran !! doesn't have a uint64 type, and not all compilers provide integers with > 64 bits... real function random_01_CB(ctr, key) - integer, parameter :: int64 = selected_int_kind(10) ! Integer with >= 64 bits - integer, intent(in) :: ctr, key ! counter & key inputs, standard integer kind + use iso_fortran_env, only : int64 + integer, intent(in) :: ctr !< ctr should be incremented each time you call the function + integer, intent(in) :: key !< key is like a seed: use a different key for each random stream integer(kind=int64) :: x, y, z ! Follows "Squares" naming convention x = (ctr + 1) * (key + 65536) ! 65536 added because keys below that don't work. From dde80bafcb9e097fbf64fdba859c007d8c5162f6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 6 Jan 2021 18:16:39 -0700 Subject: [PATCH 05/22] Avoids calling boundary_k_range at land pts --- src/tracer/MOM_neutral_diffusion.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index e995ca1972..0a413e866b 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -320,7 +320,9 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) call pass_var(hbl,G%Domain) ! get k-indices and zeta do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 - call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j)) + if (G%mask2dT(i,j) > 0.) then + call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j)) + endif enddo; enddo ! TODO: add similar code for BOTTOM boundary layer endif From 0262e930a005ec45ba7938d1b44a54185b6a3563 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 7 Jan 2021 21:10:34 -0700 Subject: [PATCH 06/22] remove omp parallel do directive for a do loop including cpu_clock calls --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index b81cf62631..d24bc890f5 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -391,8 +391,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc) call sfc_bkgnd_mixing(G, US, CS%bkgnd_mixing_csp) - !$OMP parallel do default(shared) private(dRho_int, N2_lay, N2_int, N2_bot, KT_extra, & - !$OMP KS_extra, TKE_to_Kd, maxTKE, dissip, kb) do j=js,je ! Set up variables related to the stratification. From c3ea9d0c3ee1ff982ef75396145020d7d303c6f6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 21 Jan 2021 09:47:44 -0700 Subject: [PATCH 07/22] Fixes latent heat from fprec and frunoff This patch fixes a sign bug, in both MCT and NUOPC, when accounting for the latent heat from fprec and frunnoff. Following MOM6's definition, both fprec and frunoff are > 0 into the ocean. Therefore, the latent heat associated with these terms should be negative. --- config_src/mct_driver/mom_surface_forcing_mct.F90 | 12 ++++++------ .../nuopc_driver/mom_surface_forcing_nuopc.F90 | 10 ++++++---- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/config_src/mct_driver/mom_surface_forcing_mct.F90 b/config_src/mct_driver/mom_surface_forcing_mct.F90 index 82105e040e..ef0527dd1c 100644 --- a/config_src/mct_driver/mom_surface_forcing_mct.F90 +++ b/config_src/mct_driver/mom_surface_forcing_mct.F90 @@ -486,17 +486,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! latent heat flux (W/m^2) fluxes%latent(i,j) = 0.0 - ! contribution from frozen ppt + ! contribution from frozen ppt (notice minus sign since fprec is positive into the ocean) if (associated(IOB%fprec)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif - ! contribution from frozen runoff + ! contribution from frozen runoff (notice minus sign since rofi_flux is positive into the ocean) if (associated(fluxes%frunoff)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif ! contribution from evaporation if (associated(IOB%q_flux)) then diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index 9ecf8bb01a..b219e12971 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -485,15 +485,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0) fluxes%latent(i,j) = 0.0 + ! notice minus sign since fprec is positive into the ocean if (associated(IOB%fprec)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif + ! notice minus sign since frunoff is positive into the ocean if (associated(IOB%frunoff)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * & + fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * & IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion endif if (associated(IOB%q_flux)) then From e209c0ed60ef6a0f7cf00de4d4b94c35632ca3f0 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 27 Jan 2021 12:32:20 -0700 Subject: [PATCH 08/22] incorporated flux correction factors --- config_src/nuopc_driver/mom_cap.F90 | 79 +++++++++++++++--- config_src/nuopc_driver/mom_cap_methods.F90 | 91 ++++++++++++--------- 2 files changed, 119 insertions(+), 51 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index fc6bb5035e..ca44833341 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -35,7 +35,7 @@ module MOM_cap_mod use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh use MOM_cap_time, only: AlarmInit -use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype +use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor, med2mod_areacor #ifdef CESMCOUPLED use shr_file_mod, only: shr_file_setLogUnit, shr_file_getLogUnit #endif @@ -68,6 +68,7 @@ module MOM_cap_mod use ESMF, only: ESMF_COORDSYS_SPH_DEG, ESMF_GridCreate, ESMF_INDEX_DELOCAL use ESMF, only: ESMF_MESHLOC_ELEMENT, ESMF_RC_VAL_OUTOFRANGE, ESMF_StateGet use ESMF, only: ESMF_TimePrint, ESMF_AlarmSet, ESMF_FieldGet, ESMF_Array +use ESMF, only: ESMF_FieldRegridGetArea use ESMF, only: ESMF_ArrayCreate use ESMF, only: ESMF_RC_FILE_OPEN, ESMF_RC_FILE_READ, ESMF_RC_FILE_WRITE use ESMF, only: ESMF_VMBroadcast @@ -888,16 +889,21 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) character(len=256) :: cvalue character(len=256) :: frmt ! format specifier for several error msgs character(len=512) :: err_msg ! error messages + integer :: spatialDim + integer :: numOwnedElements + type(ESMF_Array) :: elemMaskArray + real(ESMF_KIND_R8) , pointer :: ownedElemCoords(:) + real(ESMF_KIND_R8) , pointer :: lat(:), latMesh(:) + real(ESMF_KIND_R8) , pointer :: lon(:), lonMesh(:) + integer(ESMF_KIND_I4) , pointer :: mask(:), maskMesh(:) + real(ESMF_KIND_R8) :: diff_lon, diff_lat + real :: eps_omesh + real(ESMF_KIND_R8) :: L2_to_rad2 + type(ESMF_Field) :: lfield + real(ESMF_KIND_R8), allocatable :: mesh_areas(:) + real(ESMF_KIND_R8), allocatable :: model_areas(:) + real(ESMF_KIND_R8), pointer :: dataPtr_mesh_areas(:) character(len=*), parameter :: subname='(MOM_cap:InitializeRealize)' - integer :: spatialDim - integer :: numOwnedElements - type(ESMF_Array) :: elemMaskArray - real(ESMF_KIND_R8) , pointer :: ownedElemCoords(:) - real(ESMF_KIND_R8) , pointer :: lat(:), latMesh(:) - real(ESMF_KIND_R8) , pointer :: lon(:), lonMesh(:) - integer(ESMF_KIND_I4) , pointer :: mask(:), maskMesh(:) - real(ESMF_KIND_R8) :: diff_lon, diff_lat - real :: eps_omesh !-------------------------------- rc = ESMF_SUCCESS @@ -1426,7 +1432,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) line=__LINE__, & file=__FILE__)) & return - endif !--------------------------------- @@ -1450,6 +1455,58 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) endif + !--------------------------------- + ! determine flux area correction factors - module variables in mom_cap_methods + !--------------------------------- + + ! Area correction factors are ONLY valid for meshes that are read in - so do not need them for + ! grids that are calculated internally + + if (geomtype == ESMF_GEOMTYPE_MESH) then + + ! Determine mesh areas for regridding + call ESMF_MeshGet(Emesh, numOwnedElements=numOwnedElements, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_StateGet(exportState, itemName=trim(fldsFrOcn(2)%stdname), field=lfield, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_FieldRegridGetArea(lfield, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + allocate(mesh_areas(numOwnedElements)) + mesh_areas(:) = dataPtr_mesh_areas(:) + + ! Determine model areas + allocate(model_areas(numOwnedElements)) + k = 0 + L2_to_rad2 = ocean_grid%US%L_to_m**2 / ocean_grid%Rad_Earth**2 + do j = ocean_grid%jsc, ocean_grid%jec + do i = ocean_grid%isc, ocean_grid%iec + k = k + 1 ! Increment position within gindex + model_areas(k) = ocean_grid%AreaT(i,j) * L2_to_rad2 + enddo + enddo + + ! Determine flux correction factors (module variables in mom_) + allocate (mod2med_areacor(numOwnedElements)) + allocate (med2mod_areacor(numOwnedElements)) + do n = 1,numOwnedElements + if (model_areas(n) == mesh_areas(n)) then + mod2med_areacor(n) = 1._ESMF_KIND_R8 + med2mod_areacor(n) = 1._ESMF_KIND_R8 + else + mod2med_areacor(n) = model_areas(n) / mesh_areas(n) + med2mod_areacor(n) = mesh_areas(n) / model_areas(n) + if (abs(mod2med_areacor(n) - 1._ESMF_KIND_R8) > 1.e-13) then + write(6,'(a,i8,2x,d21.14,2x)')' AREACOR mom6: n, abs(mod2med_areacor(n)-1)', & + n, abs(mod2med_areacor(n) - 1._ESMF_KIND_R8) + end if + end if + end do + deallocate(model_areas) + deallocate(mesh_areas) + end if + !--------------------------------- ! Set module variable geomtype in MOM_cap_methods !--------------------------------- diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90 index 70915d0e95..d365268c0c 100644 --- a/config_src/nuopc_driver/mom_cap_methods.F90 +++ b/config_src/nuopc_driver/mom_cap_methods.F90 @@ -43,6 +43,13 @@ module MOM_cap_methods type(ESMF_GeomType_Flag) :: geomtype !< SMF type describing type of !! geometry (mesh or grid) +! area correction factors for fluxes send and received from mediator +! these actors are ONLY valid for meshes that are read in - so do not need them for +! grids that are calculated internally + +real(ESMF_KIND_R8), public, allocatable :: mod2med_areacor(:) ! ratios of model areas to input mesh areas +real(ESMF_KIND_R8), public, allocatable :: med2mod_areacor(:) ! ratios of input mesh areas to model areas + contains !> Sets module variable geometry type @@ -95,7 +102,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! near-IR, direct shortwave (W/m2) !---- call state_getimport(importState, 'mean_net_sw_ir_dir_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_nir_dir, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_nir_dir, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -147,12 +154,14 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, allocate (taux(isc:iec,jsc:jec)) allocate (tauy(isc:iec,jsc:jec)) - call state_getimport(importState, 'mean_zonal_moment_flx', isc, iec, jsc, jec, taux, rc=rc) + call state_getimport(importState, 'mean_zonal_moment_flx', isc, iec, jsc, jec, taux, & + areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - call state_getimport(importState, 'mean_merid_moment_flx', isc, iec, jsc, jec, tauy, rc=rc) + call state_getimport(importState, 'mean_merid_moment_flx', isc, iec, jsc, jec, tauy, & + areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -177,7 +186,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! sensible heat flux (W/m2) !---- call state_getimport(importState, 'mean_sensi_heat_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%t_flux, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%t_flux, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -187,7 +196,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! evaporation flux (W/m2) !---- call state_getimport(importState, 'mean_evap_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%q_flux, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%q_flux, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -197,7 +206,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! liquid precipitation (rain) !---- call state_getimport(importState, 'mean_prec_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%lprec, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lprec, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -207,7 +216,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! frozen precipitation (snow) !---- call state_getimport(importState, 'mean_fprec_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%fprec, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%fprec, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -222,7 +231,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! liquid runoff ice_ocean_boundary%lrunoff (:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Foxx_rofl', & - isc, iec, jsc, jec, ice_ocean_boundary%lrunoff,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -231,7 +240,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! ice runoff ice_ocean_boundary%frunoff (:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Foxx_rofi', & - isc, iec, jsc, jec, ice_ocean_boundary%frunoff,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%frunoff, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -240,7 +249,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! heat content of lrunoff ice_ocean_boundary%lrunoff_hflx(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_runoff_heat_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_hflx, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_hflx, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -249,7 +258,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! heat content of frunoff ice_ocean_boundary%frunoff_hflx(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_calving_heat_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%frunoff_hflx, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%frunoff_hflx, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -260,29 +269,29 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ice_ocean_boundary%salt_flux(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_salt_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%salt_flux,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%salt_flux, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - ! !---- - ! ! snow&ice melt heat flux (W/m^2) - ! !---- + !---- + ! snow&ice melt heat flux (W/m^2) + !---- ice_ocean_boundary%seaice_melt_heat(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'net_heat_flx_to_ocn', & - isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_heat,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_heat, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - ! !---- - ! ! snow&ice melt water flux (W/m^2) - ! !---- + !---- + ! snow&ice melt water flux (W/m^2) + !---- ice_ocean_boundary%seaice_melt(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_fresh_water_to_ocean_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt,rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -296,7 +305,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ice_ocean_boundary%mi(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mass_of_overlying_ice', & - isc, iec, jsc, jec, ice_ocean_boundary%mi, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -477,7 +486,7 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, enddo call State_SetExport(exportState, 'freezing_melting_potential', & - isc, iec, jsc, jec, melt_potential, ocean_grid, rc=rc) + isc, iec, jsc, jec, melt_potential, ocean_grid, areacor=mod2med_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -669,7 +678,7 @@ subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) end subroutine State_GetFldPtr_2d !> Map import state field to output array -subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, rc) +subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, areacor, rc) type(ESMF_State) , intent(in) :: state !< ESMF state character(len=*) , intent(in) :: fldname !< Field name integer , intent(in) :: isc !< The start i-index of cell centers within @@ -682,6 +691,8 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, r !! the computational domain real (ESMF_KIND_R8) , intent(inout) :: output(isc:iec,jsc:jec)!< Output 2D array logical, optional , intent(in) :: do_sum !< If true, sums the data + real (ESMF_KIND_R8), optional, intent(in) :: areacor(:) !< flux area correction factors + !! applicable to meshes integer , intent(out) :: rc !< Return code ! local variables @@ -702,10 +713,7 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, r ! get field pointer call state_getfldptr(state, trim(fldname), dataptr1d, rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! bail out ! determine output array n = 0 @@ -719,14 +727,16 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, r endif enddo enddo + if (present(areacor)) then + do n = 1,size(dataPtr1d) + dataPtr1d(n) = dataPtr1d(n) * areacor(n) + end do + end if else if (geomtype == ESMF_GEOMTYPE_GRID) then call state_getfldptr(state, trim(fldname), dataptr2d, rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! bail out lbnd1 = lbound(dataPtr2d,1) lbnd2 = lbound(dataPtr2d,2) @@ -750,7 +760,7 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, r end subroutine State_GetImport !> Map input array to export state -subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid, rc) +subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid, areacor, rc) type(ESMF_State) , intent(inout) :: state !< ESMF state character(len=*) , intent(in) :: fldname !< Field name integer , intent(in) :: isc !< The start i-index of cell centers within @@ -763,6 +773,8 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid !! the computational domain real (ESMF_KIND_R8) , intent(in) :: input(isc:iec,jsc:jec)!< Input 2D array type(ocean_grid_type) , intent(in) :: ocean_grid !< Ocean horizontal grid + real (ESMF_KIND_R8), optional, intent(in) :: areacor(:) !< flux area correction factors + !! applicable to meshes integer , intent(out) :: rc !< Return code ! local variables @@ -786,10 +798,7 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid if (geomtype == ESMF_GEOMTYPE_MESH) then call state_getfldptr(state, trim(fldname), dataptr1d, rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! bail out n = 0 do j = jsc, jec @@ -800,14 +809,16 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid dataPtr1d(n) = input(i,j) * ocean_grid%mask2dT(ig,jg) enddo enddo + if (present(areacor)) then + do n = 1,(size(dataPtr1d)) + dataPtr1d(n) = dataPtr1d(n) * areacor(n) + enddo + end if else if (geomtype == ESMF_GEOMTYPE_GRID) then call state_getfldptr(state, trim(fldname), dataptr2d, rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! bail out lbnd1 = lbound(dataPtr2d,1) lbnd2 = lbound(dataPtr2d,2) From 387166bced542edb2710e6b0ba3e39ff5470bece Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 24 Feb 2021 16:26:27 -0700 Subject: [PATCH 09/22] reduce the line length in mom_surface_forcing_mct --- config_src/mct_driver/mom_surface_forcing_mct.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/config_src/mct_driver/mom_surface_forcing_mct.F90 b/config_src/mct_driver/mom_surface_forcing_mct.F90 index ef0527dd1c..e5e1309a91 100644 --- a/config_src/mct_driver/mom_surface_forcing_mct.F90 +++ b/config_src/mct_driver/mom_surface_forcing_mct.F90 @@ -496,7 +496,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(fluxes%frunoff)) then fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * & + IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif ! contribution from evaporation if (associated(IOB%q_flux)) then From 1e2fbf463bfbcbc12ef9f404d6ac9bb5a4c12646 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 25 Feb 2021 08:04:26 -0700 Subject: [PATCH 10/22] Add missing areacor --- config_src/nuopc_driver/mom_cap_methods.F90 | 22 ++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90 index d365268c0c..e706aefa6f 100644 --- a/config_src/nuopc_driver/mom_cap_methods.F90 +++ b/config_src/nuopc_driver/mom_cap_methods.F90 @@ -19,6 +19,7 @@ module MOM_cap_methods use MOM_surface_forcing_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type use MOM_domains, only: pass_var +use MOM_error_handler, only: MOM_error, FATAL, is_root_pe use mpp_domains_mod, only: mpp_get_compute_domain ! By default make data private @@ -101,6 +102,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ! near-IR, direct shortwave (W/m2) !---- + call state_getimport(importState, 'mean_net_sw_ir_dir_flx', & isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_nir_dir, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -112,7 +114,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! near-IR, diffuse shortwave (W/m2) !---- call state_getimport(importState, 'mean_net_sw_ir_dif_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_nir_dif, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_nir_dif, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -122,7 +124,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! visible, direct shortwave (W/m2) !---- call state_getimport(importState, 'mean_net_sw_vis_dir_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_vis_dir, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_vis_dir, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -132,7 +134,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! visible, diffuse shortwave (W/m2) !---- call state_getimport(importState, 'mean_net_sw_vis_dif_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_vis_dif, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_vis_dif, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -142,7 +144,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! Net longwave radiation (W/m2) ! ------- call state_getimport(importState, 'mean_net_lw_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%lw_flux, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lw_flux, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -715,6 +717,13 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, a call state_getfldptr(state, trim(fldname), dataptr1d, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! bail out + ! option to apply area correction + if (present(areacor)) then + do n = 1,size(dataPtr1d) + dataPtr1d(n) = dataPtr1d(n) * areacor(n) + end do + end if + ! determine output array n = 0 do j = jsc,jec @@ -727,11 +736,6 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, a endif enddo enddo - if (present(areacor)) then - do n = 1,size(dataPtr1d) - dataPtr1d(n) = dataPtr1d(n) * areacor(n) - end do - end if else if (geomtype == ESMF_GEOMTYPE_GRID) then From 87298d770f5367794a1603dc33255ed9fe1cea5f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 25 Feb 2021 10:19:49 -0700 Subject: [PATCH 11/22] Fix LBD module after merging main_candidate_2021-02-19 --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 0d6dc14bd7..8f022821ea 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -274,7 +274,6 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) endif ! post tendency of tracer content - !### This seems to be dimensionally inconsistent with the calculation of tendency above. if (tracer%id_lbdxy_cont > 0) then call post_data(tracer%id_lbdxy_cont, tendency, CS%diag) endif @@ -293,7 +292,6 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) ! post tendency of tracer concentration; this step must be ! done after posting tracer content tendency, since we alter ! the tendency array and its units. - !### This seems to be dimensionally inconsistent with the calculation of tendency above. if (tracer%id_lbdxy_conc > 0) then do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + CS%H_subroundoff ) From 5667b70f3623a2d9f62ef58449d4da7528d830b8 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 26 Feb 2021 14:00:45 -0700 Subject: [PATCH 12/22] Remove unused imports and delete blank line --- config_src/nuopc_driver/mom_cap_methods.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90 index e706aefa6f..625bde4cfe 100644 --- a/config_src/nuopc_driver/mom_cap_methods.F90 +++ b/config_src/nuopc_driver/mom_cap_methods.F90 @@ -19,7 +19,6 @@ module MOM_cap_methods use MOM_surface_forcing_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type use MOM_domains, only: pass_var -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe use mpp_domains_mod, only: mpp_get_compute_domain ! By default make data private @@ -102,7 +101,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ! near-IR, direct shortwave (W/m2) !---- - call state_getimport(importState, 'mean_net_sw_ir_dir_flx', & isc, iec, jsc, jec, ice_ocean_boundary%sw_flux_nir_dir, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & From 7ab7daaf28104f5237320d9b9555e8531a4b6a3e Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Tue, 23 Mar 2021 19:40:16 -0600 Subject: [PATCH 13/22] updates and bugfixes for area correction factors --- config_src/nuopc_driver/mom_cap.F90 | 64 +++++++++++++-------- config_src/nuopc_driver/mom_cap_methods.F90 | 43 +++++++------- 2 files changed, 63 insertions(+), 44 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index ca44833341..10c73150cf 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -38,6 +38,7 @@ module MOM_cap_mod use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor, med2mod_areacor #ifdef CESMCOUPLED use shr_file_mod, only: shr_file_setLogUnit, shr_file_getLogUnit +use shr_mpi_mod, only : shr_mpi_min, shr_mpi_max #endif use time_utils_mod, only: esmf2fms_time @@ -71,7 +72,7 @@ module MOM_cap_mod use ESMF, only: ESMF_FieldRegridGetArea use ESMF, only: ESMF_ArrayCreate use ESMF, only: ESMF_RC_FILE_OPEN, ESMF_RC_FILE_READ, ESMF_RC_FILE_WRITE -use ESMF, only: ESMF_VMBroadcast +use ESMF, only: ESMF_VMBroadcast, ESMF_VMReduce, ESMF_REDUCE_MAX, ESMF_REDUCE_MIN use ESMF, only: ESMF_AlarmCreate, ESMF_ClockGetAlarmList, ESMF_AlarmList_Flag use ESMF, only: ESMF_AlarmGet, ESMF_AlarmIsCreated, ESMF_ALARMLIST_ALL, ESMF_AlarmIsEnabled use ESMF, only: ESMF_STATEITEM_NOTFOUND, ESMF_FieldWrite @@ -903,6 +904,14 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) real(ESMF_KIND_R8), allocatable :: mesh_areas(:) real(ESMF_KIND_R8), allocatable :: model_areas(:) real(ESMF_KIND_R8), pointer :: dataPtr_mesh_areas(:) + real(ESMF_KIND_R8) :: max_mod2med_areacor + real(ESMF_KIND_R8) :: max_med2mod_areacor + real(ESMF_KIND_R8) :: min_mod2med_areacor + real(ESMF_KIND_R8) :: min_med2mod_areacor + real(ESMF_KIND_R8) :: max_mod2med_areacor_glob + real(ESMF_KIND_R8) :: max_med2mod_areacor_glob + real(ESMF_KIND_R8) :: min_mod2med_areacor_glob + real(ESMF_KIND_R8) :: min_med2mod_areacor_glob character(len=*), parameter :: subname='(MOM_cap:InitializeRealize)' !-------------------------------- @@ -1458,7 +1467,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) !--------------------------------- ! determine flux area correction factors - module variables in mom_cap_methods !--------------------------------- - ! Area correction factors are ONLY valid for meshes that are read in - so do not need them for ! grids that are calculated internally @@ -1467,6 +1475,13 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! Determine mesh areas for regridding call ESMF_MeshGet(Emesh, numOwnedElements=numOwnedElements, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + allocate (mod2med_areacor(numOwnedElements)) + allocate (med2mod_areacor(numOwnedElements)) + mod2med_areacor(:) = 1._ESMF_KIND_R8 + med2mod_areacor(:) = 1._ESMF_KIND_R8 + +#ifdef CESMCOUPLED call ESMF_StateGet(exportState, itemName=trim(fldsFrOcn(2)%stdname), field=lfield, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return call ESMF_FieldRegridGetArea(lfield, rc=rc) @@ -1474,37 +1489,40 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return allocate(mesh_areas(numOwnedElements)) - mesh_areas(:) = dataPtr_mesh_areas(:) - - ! Determine model areas allocate(model_areas(numOwnedElements)) + + ! Determine model areas and flux correction factors (module variables in mom_) k = 0 L2_to_rad2 = ocean_grid%US%L_to_m**2 / ocean_grid%Rad_Earth**2 do j = ocean_grid%jsc, ocean_grid%jec do i = ocean_grid%isc, ocean_grid%iec k = k + 1 ! Increment position within gindex + mesh_areas(k) = dataPtr_mesh_areas(k) model_areas(k) = ocean_grid%AreaT(i,j) * L2_to_rad2 + mod2med_areacor(k) = model_areas(k) / mesh_areas(k) + med2mod_areacor(k) = mesh_areas(k) / model_areas(k) enddo enddo - - ! Determine flux correction factors (module variables in mom_) - allocate (mod2med_areacor(numOwnedElements)) - allocate (med2mod_areacor(numOwnedElements)) - do n = 1,numOwnedElements - if (model_areas(n) == mesh_areas(n)) then - mod2med_areacor(n) = 1._ESMF_KIND_R8 - med2mod_areacor(n) = 1._ESMF_KIND_R8 - else - mod2med_areacor(n) = model_areas(n) / mesh_areas(n) - med2mod_areacor(n) = mesh_areas(n) / model_areas(n) - if (abs(mod2med_areacor(n) - 1._ESMF_KIND_R8) > 1.e-13) then - write(6,'(a,i8,2x,d21.14,2x)')' AREACOR mom6: n, abs(mod2med_areacor(n)-1)', & - n, abs(mod2med_areacor(n) - 1._ESMF_KIND_R8) - end if - end if - end do - deallocate(model_areas) deallocate(mesh_areas) + deallocate(model_areas) + + ! Write diagnostic output for correction factors + min_mod2med_areacor = minval(mod2med_areacor) + max_mod2med_areacor = maxval(mod2med_areacor) + min_med2mod_areacor = minval(med2mod_areacor) + max_med2mod_areacor = maxval(med2mod_areacor) + call shr_mpi_max(max_mod2med_areacor, max_mod2med_areacor_glob, mpicom) + call shr_mpi_min(min_mod2med_areacor, min_mod2med_areacor_glob, mpicom) + call shr_mpi_max(max_med2mod_areacor, max_med2mod_areacor_glob, mpicom) + call shr_mpi_min(min_med2mod_areacor, min_med2mod_areacor_glob, mpicom) + if (localPet == 0) then + write(logunit,'(2A,2g23.15,A )') trim(subname),' : min_mod2med_areacor, max_mod2med_areacor ',& + min_mod2med_areacor_glob, max_mod2med_areacor_glob, 'MOM6' + write(logunit,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& + min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'MOM6' + end if +#endif + end if !--------------------------------- diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90 index 625bde4cfe..2b5b54e55a 100644 --- a/config_src/nuopc_driver/mom_cap_methods.F90 +++ b/config_src/nuopc_driver/mom_cap_methods.F90 @@ -186,7 +186,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! sensible heat flux (W/m2) !---- call state_getimport(importState, 'mean_sensi_heat_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%t_flux, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%t_flux, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -196,7 +196,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! evaporation flux (W/m2) !---- call state_getimport(importState, 'mean_evap_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%q_flux, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%q_flux, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -206,7 +206,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! liquid precipitation (rain) !---- call state_getimport(importState, 'mean_prec_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%lprec, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lprec, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -216,7 +216,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! frozen precipitation (snow) !---- call state_getimport(importState, 'mean_fprec_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%fprec, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%fprec, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -231,7 +231,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! liquid runoff ice_ocean_boundary%lrunoff (:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Foxx_rofl', & - isc, iec, jsc, jec, ice_ocean_boundary%lrunoff, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -240,7 +240,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! ice runoff ice_ocean_boundary%frunoff (:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Foxx_rofi', & - isc, iec, jsc, jec, ice_ocean_boundary%frunoff, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%frunoff, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -249,7 +249,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! heat content of lrunoff ice_ocean_boundary%lrunoff_hflx(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_runoff_heat_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_hflx, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_hflx, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -258,7 +258,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! heat content of frunoff ice_ocean_boundary%frunoff_hflx(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_calving_heat_flx', & - isc, iec, jsc, jec, ice_ocean_boundary%frunoff_hflx, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%frunoff_hflx, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -269,7 +269,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ice_ocean_boundary%salt_flux(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_salt_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%salt_flux, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%salt_flux, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -280,7 +280,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ice_ocean_boundary%seaice_melt_heat(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'net_heat_flx_to_ocn', & - isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_heat, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt_heat, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -291,7 +291,7 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ice_ocean_boundary%seaice_melt(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mean_fresh_water_to_ocean_rate', & - isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt, areacor=mod2med_areacor, rc=rc) + isc, iec, jsc, jec, ice_ocean_boundary%seaice_melt, areacor=med2mod_areacor, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -715,22 +715,23 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, a call state_getfldptr(state, trim(fldname), dataptr1d, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return ! bail out - ! option to apply area correction - if (present(areacor)) then - do n = 1,size(dataPtr1d) - dataPtr1d(n) = dataPtr1d(n) * areacor(n) - end do - end if - - ! determine output array + ! determine output array and apply area correction if present n = 0 do j = jsc,jec do i = isc,iec n = n + 1 if (present(do_sum)) then - output(i,j) = output(i,j) + dataPtr1d(n) + if (present(areacor)) then + output(i,j) = output(i,j) + dataPtr1d(n) * areacor(n) + else + output(i,j) = output(i,j) + dataPtr1d(n) + end if else - output(i,j) = dataPtr1d(n) + if (present(areacor)) then + output(i,j) = dataPtr1d(n) * areacor(n) + else + output(i,j) = dataPtr1d(n) + end if endif enddo enddo From b8aaa8824694bd33add29c8cd4b03bb77ccac62c Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Tue, 23 Mar 2021 20:46:26 -0600 Subject: [PATCH 14/22] changed computation of model_areas --- config_src/nuopc_driver/mom_cap.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 10c73150cf..862147d0d8 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -1498,7 +1498,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) do i = ocean_grid%isc, ocean_grid%iec k = k + 1 ! Increment position within gindex mesh_areas(k) = dataPtr_mesh_areas(k) - model_areas(k) = ocean_grid%AreaT(i,j) * L2_to_rad2 + model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 mod2med_areacor(k) = model_areas(k) / mesh_areas(k) med2mod_areacor(k) = mesh_areas(k) / model_areas(k) enddo From 973632cd6c8e7675f2e5f593975892712981abc3 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 24 Mar 2021 20:55:04 -0600 Subject: [PATCH 15/22] calculate area correction factors only over non-zero mask values --- config_src/nuopc_driver/mom_cap.F90 | 128 ++++++++++++++-------------- 1 file changed, 63 insertions(+), 65 deletions(-) diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index 862147d0d8..8e2acb5195 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -1130,10 +1130,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) end if end do - deallocate(ownedElemCoords) - deallocate(lonMesh , lon ) - deallocate(latMesh , lat ) - deallocate(maskMesh, mask) ! realize the import and export fields using the mesh call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", mesh=Emesh, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -1147,6 +1143,69 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) file=__FILE__)) & return + !--------------------------------- + ! determine flux area correction factors - module variables in mom_cap_methods + !--------------------------------- + ! Area correction factors are ONLY valid for meshes that are read in - so do not need them for + ! grids that are calculated internally + + ! Determine mesh areas for regridding + call ESMF_MeshGet(Emesh, numOwnedElements=numOwnedElements, spatialDim=spatialDim, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + allocate (mod2med_areacor(numOwnedElements)) + allocate (med2mod_areacor(numOwnedElements)) + mod2med_areacor(:) = 1._ESMF_KIND_R8 + med2mod_areacor(:) = 1._ESMF_KIND_R8 + +#ifdef CESMCOUPLED + ! Determine model areas and flux correction factors (module variables in mom_) + call ESMF_StateGet(exportState, itemName=trim(fldsFrOcn(2)%stdname), field=lfield, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_FieldRegridGetArea(lfield, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + + allocate(mesh_areas(numOwnedElements)) + allocate(model_areas(numOwnedElements)) + k = 0 + do j = ocean_grid%jsc, ocean_grid%jec + do i = ocean_grid%isc, ocean_grid%iec + k = k + 1 ! Increment position within gindex + if (mask(k) /= 0) then + mesh_areas(k) = dataPtr_mesh_areas(k) + model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 + mod2med_areacor(k) = model_areas(k) / mesh_areas(k) + med2mod_areacor(k) = mesh_areas(k) / model_areas(k) + end if + end do + end do + deallocate(mesh_areas) + deallocate(model_areas) + + ! Write diagnostic output for correction factors + min_mod2med_areacor = minval(mod2med_areacor) + max_mod2med_areacor = maxval(mod2med_areacor) + min_med2mod_areacor = minval(med2mod_areacor) + max_med2mod_areacor = maxval(med2mod_areacor) + call shr_mpi_max(max_mod2med_areacor, max_mod2med_areacor_glob, mpicom) + call shr_mpi_min(min_mod2med_areacor, min_mod2med_areacor_glob, mpicom) + call shr_mpi_max(max_med2mod_areacor, max_med2mod_areacor_glob, mpicom) + call shr_mpi_min(min_med2mod_areacor, min_med2mod_areacor_glob, mpicom) + if (localPet == 0) then + write(logunit,'(2A,2g23.15,A )') trim(subname),' : min_mod2med_areacor, max_mod2med_areacor ',& + min_mod2med_areacor_glob, max_mod2med_areacor_glob, 'MOM6' + write(logunit,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& + min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'MOM6' + end if +#endif + + deallocate(ownedElemCoords) + deallocate(lonMesh , lon ) + deallocate(latMesh , lat ) + deallocate(maskMesh, mask) + else if (geomtype == ESMF_GEOMTYPE_GRID) then !--------------------------------- @@ -1464,67 +1523,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) endif - !--------------------------------- - ! determine flux area correction factors - module variables in mom_cap_methods - !--------------------------------- - ! Area correction factors are ONLY valid for meshes that are read in - so do not need them for - ! grids that are calculated internally - - if (geomtype == ESMF_GEOMTYPE_MESH) then - - ! Determine mesh areas for regridding - call ESMF_MeshGet(Emesh, numOwnedElements=numOwnedElements, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - - allocate (mod2med_areacor(numOwnedElements)) - allocate (med2mod_areacor(numOwnedElements)) - mod2med_areacor(:) = 1._ESMF_KIND_R8 - med2mod_areacor(:) = 1._ESMF_KIND_R8 - -#ifdef CESMCOUPLED - call ESMF_StateGet(exportState, itemName=trim(fldsFrOcn(2)%stdname), field=lfield, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_FieldRegridGetArea(lfield, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - allocate(mesh_areas(numOwnedElements)) - allocate(model_areas(numOwnedElements)) - - ! Determine model areas and flux correction factors (module variables in mom_) - k = 0 - L2_to_rad2 = ocean_grid%US%L_to_m**2 / ocean_grid%Rad_Earth**2 - do j = ocean_grid%jsc, ocean_grid%jec - do i = ocean_grid%isc, ocean_grid%iec - k = k + 1 ! Increment position within gindex - mesh_areas(k) = dataPtr_mesh_areas(k) - model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 - mod2med_areacor(k) = model_areas(k) / mesh_areas(k) - med2mod_areacor(k) = mesh_areas(k) / model_areas(k) - enddo - enddo - deallocate(mesh_areas) - deallocate(model_areas) - - ! Write diagnostic output for correction factors - min_mod2med_areacor = minval(mod2med_areacor) - max_mod2med_areacor = maxval(mod2med_areacor) - min_med2mod_areacor = minval(med2mod_areacor) - max_med2mod_areacor = maxval(med2mod_areacor) - call shr_mpi_max(max_mod2med_areacor, max_mod2med_areacor_glob, mpicom) - call shr_mpi_min(min_mod2med_areacor, min_mod2med_areacor_glob, mpicom) - call shr_mpi_max(max_med2mod_areacor, max_med2mod_areacor_glob, mpicom) - call shr_mpi_min(min_med2mod_areacor, min_med2mod_areacor_glob, mpicom) - if (localPet == 0) then - write(logunit,'(2A,2g23.15,A )') trim(subname),' : min_mod2med_areacor, max_mod2med_areacor ',& - min_mod2med_areacor_glob, max_mod2med_areacor_glob, 'MOM6' - write(logunit,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& - min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'MOM6' - end if -#endif - - end if - !--------------------------------- ! Set module variable geomtype in MOM_cap_methods !--------------------------------- From a4bf34459c4e9490eb435827d0cf3d0bb6bca08a Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Mon, 29 Mar 2021 10:49:30 -0600 Subject: [PATCH 16/22] add if clause to a set_diffusivity OMP block including clock calls. --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 99dee11b9a..0cb39b2f15 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -417,7 +417,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! parameterization of Kd. !$OMP parallel do default(shared) private(dRho_int,N2_lay,Kd_lay_2d,Kd_int_2d,Kv_bkgnd,N2_int,& - !$OMP N2_bot,KT_extra,KS_extra,TKE_to_Kd,maxTKE,dissip,kb) + !$OMP N2_bot,KT_extra,KS_extra,TKE_to_Kd,maxTKE,dissip,kb)& + !$OMP if(.not. CS%use_CVMix_ddiff) do j=js,je ! Set up variables related to the stratification. From 69c86cd9dc799d995471c34133e175c5da95524c Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Mon, 5 Apr 2021 08:39:30 -0600 Subject: [PATCH 17/22] add support for threading in cmeps --- config_src/drivers/nuopc_cap/mom_cap.F90 | 62 ++++++++++++++++++++++-- 1 file changed, 59 insertions(+), 3 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 15d9415471..ca0b40d509 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -47,7 +47,7 @@ module MOM_cap_mod use, intrinsic :: iso_fortran_env, only: output_unit -use ESMF, only: ESMF_ClockAdvance, ESMF_ClockGet, ESMF_ClockPrint +use ESMF, only: ESMF_ClockAdvance, ESMF_ClockGet, ESMF_ClockPrint, ESMF_VMget use ESMF, only: ESMF_ClockGetAlarm, ESMF_ClockGetNextTime, ESMF_ClockAdvance use ESMF, only: ESMF_ClockSet, ESMF_Clock, ESMF_GeomType_Flag, ESMF_LOGMSG_INFO use ESMF, only: ESMF_Grid, ESMF_GridCreate, ESMF_GridAddCoord @@ -96,6 +96,7 @@ module MOM_cap_mod use NUOPC_Model, only: model_label_SetRunClock => label_SetRunClock use NUOPC_Model, only: model_label_Finalize => label_Finalize use NUOPC_Model, only: SetVM +!$use omp_lib , only : omp_set_num_threads implicit none; private @@ -143,7 +144,8 @@ module MOM_cap_mod integer :: scalar_field_count = 0 integer :: scalar_field_idx_grid_nx = 0 integer :: scalar_field_idx_grid_ny = 0 -character(len=*),parameter :: u_FILE_u = & +integer :: nthrds !< number of openmp threads per task +character(len=*),parameter :: u_file_u = & __FILE__ #ifdef CESMCOUPLED @@ -399,6 +401,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) logical :: existflag integer :: userRc integer :: localPet + integer :: localPeCount integer :: iostat integer :: readunit character(len=512) :: restartfile ! Path/Name of restart file @@ -435,7 +438,31 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) CALL ESMF_TimeIntervalGet(TINT, S=DT_OCEAN, RC=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - !TODO: next two lines not present in NCAR + !--------------------------------- + ! openmp threads + !--------------------------------- + + call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + + if(localPeCount == 1) then + call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + read(cvalue,*) nthrds + else + nthrds = localPeCount + endif + +!$ call omp_set_num_threads(nthrds) + print *,__FILE__,__LINE__,nthrds + call fms_init(mpi_comm_mom) call constants_init call field_manager_init @@ -790,6 +817,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) real(ESMF_KIND_R8), pointer :: dataPtr_ycor(:,:) integer :: mpicom integer :: localPet + integer :: localPeCount integer :: lsize integer :: ig,jg, ni,nj,k integer, allocatable :: gindex(:) ! global index space @@ -847,6 +875,32 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call ESMF_VMGet(vm, petCount=npet, mpiCommunicator=mpicom, localPet=localPet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + !--------------------------------- + ! openmp threads + !--------------------------------- + + call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + + if(localPeCount == 1) then + call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + read(cvalue,*) nthrds + else + nthrds = localPeCount + endif + +!$ call omp_set_num_threads(nthrds) + print *,__FILE__,__LINE__,nthrds + + !--------------------------------- ! global mom grid size !--------------------------------- @@ -1459,6 +1513,8 @@ subroutine ModelAdvance(gcomp, rc) call shr_file_setLogUnit (logunit) +!$ call omp_set_num_threads(nthrds) + ! query the Component for its clock, importState and exportState call ESMF_GridCompGet(gcomp, clock=clock, importState=importState, & exportState=exportState, rc=rc) From 2da6ee08ec78c903a9d7e6ee52b0998a92bdbf80 Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Mon, 5 Apr 2021 10:27:40 -0600 Subject: [PATCH 18/22] initialize localPet --- config_src/drivers/nuopc_cap/mom_cap.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index ca0b40d509..7379bd099a 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -426,7 +426,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ESMF_VMGetCurrent(vm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMGet(VM, mpiCommunicator=mpi_comm_mom, rc=rc) + call ESMF_VMGet(VM, mpiCommunicator=mpi_comm_mom, localPet=localPet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_ClockGet(CLOCK, currTIME=MyTime, TimeStep=TINT, RC=rc) From 8a8d9ea25ac065c6f462cb5e2fc27fed7f05038e Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Fri, 9 Apr 2021 07:06:45 -0600 Subject: [PATCH 19/22] response to review --- config_src/drivers/nuopc_cap/mom_cap.F90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 7379bd099a..2f6dfa86d3 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -145,7 +145,7 @@ module MOM_cap_mod integer :: scalar_field_idx_grid_nx = 0 integer :: scalar_field_idx_grid_ny = 0 integer :: nthrds !< number of openmp threads per task -character(len=*),parameter :: u_file_u = & +character(len=*),parameter :: u_FILE_u = & __FILE__ #ifdef CESMCOUPLED @@ -461,7 +461,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) endif !$ call omp_set_num_threads(nthrds) - print *,__FILE__,__LINE__,nthrds call fms_init(mpi_comm_mom) call constants_init @@ -898,8 +897,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) endif !$ call omp_set_num_threads(nthrds) - print *,__FILE__,__LINE__,nthrds - !--------------------------------- ! global mom grid size From 3b121cf233e073248b1d555debad1a035859169a Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 21 Apr 2021 16:25:34 -0600 Subject: [PATCH 20/22] Replaces ESMF_LogFoundError to ChkErr --- config_src/drivers/nuopc_cap/mom_cap.F90 | 32 +++++++----------------- 1 file changed, 9 insertions(+), 23 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index e81d220d0f..4776183c90 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -53,7 +53,7 @@ module MOM_cap_mod use ESMF, only: ESMF_Grid, ESMF_GridCreate, ESMF_GridAddCoord use ESMF, only: ESMF_GridGetCoord, ESMF_GridAddItem, ESMF_GridGetItem use ESMF, only: ESMF_GridComp, ESMF_GridCompSetEntryPoint, ESMF_GridCompGet -use ESMF, only: ESMF_LogFoundError, ESMF_LogWrite, ESMF_LogSetError +use ESMF, only: ESMF_LogWrite, ESMF_LogSetError use ESMF, only: ESMF_LOGERR_PASSTHRU, ESMF_KIND_R8, ESMF_RC_VAL_WRONG use ESMF, only: ESMF_GEOMTYPE_MESH, ESMF_GEOMTYPE_GRID, ESMF_SUCCESS use ESMF, only: ESMF_METHOD_INITIALIZE, ESMF_MethodRemove, ESMF_State @@ -464,18 +464,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !--------------------------------- call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - + if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) nthrds else nthrds = localPeCount @@ -900,18 +893,11 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) !--------------------------------- call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - + if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) nthrds else nthrds = localPeCount @@ -1075,7 +1061,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! Determine mesh areas for regridding call ESMF_MeshGet(Emesh, numOwnedElements=numOwnedElements, spatialDim=spatialDim, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate (mod2med_areacor(numOwnedElements)) allocate (med2mod_areacor(numOwnedElements)) @@ -1085,11 +1071,11 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) #ifdef CESMCOUPLED ! Determine model areas and flux correction factors (module variables in mom_) call ESMF_StateGet(exportState, itemName=trim(fldsFrOcn(2)%stdname), field=lfield, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldRegridGetArea(lfield, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate(mesh_areas(numOwnedElements)) allocate(model_areas(numOwnedElements)) From 6196c039fb5d40759e77b8cb317f136f0aaa430e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 22 Apr 2021 15:47:18 -0600 Subject: [PATCH 21/22] Add isPresent and isSet when retrieving nthreads --- config_src/drivers/nuopc_cap/mom_cap.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 4776183c90..f5f7985e15 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -467,7 +467,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then - call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) + call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) nthrds else @@ -822,6 +823,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) integer :: lbnd3,ubnd3,lbnd4,ubnd4 integer :: nblocks_tot logical :: found + logical :: isPresent, isSet integer(ESMF_KIND_I4), pointer :: dataPtr_mask(:,:) real(ESMF_KIND_R8), pointer :: dataPtr_area(:,:) real(ESMF_KIND_R8), pointer :: dataPtr_xcen(:,:) @@ -896,7 +898,8 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then - call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) + call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) nthrds else From e6ce6a8865c83d67455402b8e20d5023dd5ad157 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 26 Apr 2021 11:50:36 -0600 Subject: [PATCH 22/22] use isPresent,isSet to conditionally retrieve configuration --- config_src/drivers/nuopc_cap/mom_cap.F90 | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index f5f7985e15..2d79674606 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -418,6 +418,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) character(len=512) :: diro character(len=512) :: logfile character(ESMF_MAXSTR) :: cvalue + character(len=64) :: logmsg logical :: isPresent, isPresentDiro, isPresentLogfile, isSet logical :: existflag integer :: userRc @@ -467,13 +468,19 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then - call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, & + call NUOPC_CompAttributeGet(gcomp, name="nthreads", value=cvalue, & isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - read(cvalue,*) nthrds + if (isPresent .and. isSet) then + read(cvalue,*) nthrds + else + nthrds = localPeCount + endif else nthrds = localPeCount endif + write(logmsg,*) nthrds + call ESMF_LogWrite(trim(subname)//': nthreads = '//trim(logmsg), ESMF_LOGMSG_INFO) !$ call omp_set_num_threads(nthrds) @@ -898,10 +905,14 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then - call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, & + call NUOPC_CompAttributeGet(gcomp, name="nthreads", value=cvalue, & isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - read(cvalue,*) nthrds + if (isPresent .and. isSet) then + read(cvalue,*) nthrds + else + nthrds = localPeCount + endif else nthrds = localPeCount endif