diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 7dc95e13fa..2d79674606 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -35,22 +35,25 @@ 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, state_diagnose +use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, mod2med_areacor +use MOM_cap_methods, only: med2mod_areacor, state_diagnose use MOM_cap_methods, only: ChkErr + #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 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 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 @@ -69,9 +72,10 @@ 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 +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 @@ -93,6 +97,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 @@ -141,6 +146,7 @@ module MOM_cap_mod integer :: scalar_field_count = 0 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 = & __FILE__ @@ -412,10 +418,12 @@ 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 integer :: localPet + integer :: localPeCount integer :: iostat integer :: readunit character(len=512) :: restartfile ! Path/Name of restart file @@ -440,7 +448,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) @@ -452,7 +460,30 @@ 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 (ChkErr(rc,__LINE__,u_FILE_u)) return + + if(localPeCount == 1) then + call NUOPC_CompAttributeGet(gcomp, name="nthreads", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + 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) + call fms_init(mpi_comm_mom) call constants_init call field_manager_init @@ -799,6 +830,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(:,:) @@ -807,6 +839,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 @@ -814,16 +847,29 @@ 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(:) + 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)' - 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 @@ -851,6 +897,28 @@ 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 (ChkErr(rc,__LINE__,u_FILE_u)) return + + if(localPeCount == 1) then + call NUOPC_CompAttributeGet(gcomp, name="nthreads", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) nthrds + else + nthrds = localPeCount + endif + else + nthrds = localPeCount + endif + +!$ call omp_set_num_threads(nthrds) + !--------------------------------- ! global mom grid size !--------------------------------- @@ -992,10 +1060,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 (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1003,6 +1067,69 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call MOM_RealizeFields(exportState, fldsFrOcn_num, fldsFrOcn, "Ocn export", mesh=Emesh, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) 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 (ChkErr(rc,__LINE__,u_FILE_u)) 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 (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldRegridGetArea(lfield, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) 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 !--------------------------------- @@ -1229,7 +1356,6 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call MOM_RealizeFields(exportState, fldsFrOcn_num, fldsFrOcn, "Ocn export", grid=gridOut, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - endif !--------------------------------- @@ -1405,6 +1531,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) diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 78014f1c63..0a4e12c647 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -46,8 +46,13 @@ module MOM_cap_methods type(ESMF_GeomType_Flag) :: geomtype !< SMF type describing type of !! geometry (mesh or grid) -character(len=*),parameter :: u_FILE_u = & - __FILE__ +! 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 +character(len=*),parameter :: u_FILE_u = __FILE__ contains @@ -100,35 +105,35 @@ 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 (ChkErr(rc,__LINE__,u_FILE_u)) return !---- ! 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 (ChkErr(rc,__LINE__,u_FILE_u)) return !---- ! 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 (ChkErr(rc,__LINE__,u_FILE_u)) return !---- ! 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 (ChkErr(rc,__LINE__,u_FILE_u)) return ! ------- ! 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 (ChkErr(rc,__LINE__,u_FILE_u)) return !---- @@ -137,12 +142,13 @@ 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 (ChkErr(rc,__LINE__,u_FILE_u)) return - 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 (ChkErr(rc,__LINE__,u_FILE_u)) return - ! rotate taux and tauy from true zonal/meridional to local coordinates do j = jsc, jec jg = j + ocean_grid%jsc - jsc @@ -161,28 +167,28 @@ 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---- ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---- ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---- ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---- @@ -194,25 +200,25 @@ 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---- @@ -220,23 +226,23 @@ 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! !---- - ! ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! !---- - ! ! 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=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return !---- @@ -247,10 +253,9 @@ 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 (ChkErr(rc,__LINE__,u_FILE_u)) return - !---- ! Partitioned Stokes Drift Components !---- @@ -451,7 +456,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 (ChkErr(rc,__LINE__,u_FILE_u)) return deallocate(melt_potential) @@ -619,7 +624,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 @@ -632,6 +637,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 @@ -654,15 +661,23 @@ subroutine State_GetImport(state, fldname, isc, iec, jsc, jec, output, do_sum, r call state_getfldptr(state, trim(fldname), dataptr1d, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! 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 @@ -694,7 +709,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 @@ -707,6 +722,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 @@ -741,6 +758,11 @@ 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 diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index 09cb23056d..709fd27731 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -12,6 +12,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 @@ -55,6 +56,29 @@ 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) + 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. + 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 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. diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index b51813678b..8f022821ea 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -226,10 +226,10 @@ 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 - !### Probably this needs to be multiplied by (h(i,j,k) + GV%H_subroundoff) for consistency - ! the way it is used later in this routine. - tendency(i,j,k) = (tracer%t(i,j,k)-tracer_old(i,j,k)) * 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 @@ -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 ) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index d8edff3751..1bf7401cbd 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -327,8 +327,10 @@ 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, GV%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j)) - enddo ; enddo + 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