diff --git a/.gitignore b/.gitignore index 25f7524d1c..c57b950fc2 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,15 @@ html +# Build output +*.o +*.mod +MOM6 +build/ +deps/ +pkg/MARBL + + # Autoconf output aclocal.m4 autom4te.cache/ diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 3574943918..83eddf7265 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -293,7 +293,7 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) then read(value,*) dbug - end if + endif write(logmsg,'(i6)') dbug call ESMF_LogWrite('MOM_cap:dbug = '//trim(logmsg), ESMF_LOGMSG_INFO) @@ -370,7 +370,7 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) write(logmsg,*) use_mommesh call ESMF_LogWrite('MOM_cap:use_mommesh = '//trim(logmsg), ESMF_LOGMSG_INFO) - if(use_mommesh)then + if (use_mommesh) then geomtype = ESMF_GEOMTYPE_MESH call NUOPC_CompAttributeGet(gcomp, name='mesh_ocn', isPresent=isPresent, isSet=isSet, rc=rc) if (.not. isPresent .and. .not. isSet) then @@ -443,6 +443,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) character(len=32) :: calendar character(len=:), allocatable :: rpointer_filename integer :: inst_index + logical :: i2o_per_cat real(8) :: MPI_Wtime, timeiads !-------------------------------- @@ -560,6 +561,34 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) time0 = set_date (YEAR,MONTH,DAY,HOUR,MINUTE,SECOND) + !----------------- + ! optional input from cice columns due to ice thickness categories + !----------------- + + Ice_ocean_boundary%ice_ncat = 0 + if (cesm_coupled) then + ! Note that flds_i2o_per_cat is set by the env_run.xml variable CPL_I2O_PER_CAT + ! This xml variable is set by MOM_interface's buildnml script; it has the same + ! value as USE_MARBL in the case + call NUOPC_CompAttributeGet(gcomp, name='flds_i2o_per_cat', value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) i2o_per_cat + if (is_root_pe()) then + write(stdout,*) 'i2o_per_cat = ',i2o_per_cat + endif + + ! Note that ice_ncat is set by the env_run.xml variable ICE_NCAT which is set + ! by the ice component (default is 1) + if (i2o_per_cat) then + call NUOPC_CompAttributeGet(gcomp, name='ice_ncat', value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) Ice_ocean_boundary%ice_ncat + endif + if (is_root_pe()) then + write(stdout,*) 'ice_ncat = ', Ice_ocean_boundary%ice_ncat + endif + end if + if (is_root_pe()) then write(stdout,*) subname//'start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,second endif @@ -663,74 +692,70 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call get_domain_extent(ocean_public%domain, isc, iec, jsc, jec) - allocate ( Ice_ocean_boundary% u_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% v_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% t_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% q_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% salt_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% lw_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_vis_dir (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_vis_dif (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_nir_dir (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_nir_dif (isc:iec,jsc:jec), & - Ice_ocean_boundary% lprec (isc:iec,jsc:jec), & - Ice_ocean_boundary% fprec (isc:iec,jsc:jec), & - Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),& - Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), & - Ice_ocean_boundary% mi (isc:iec,jsc:jec), & - Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec), & - Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec), & - Ice_ocean_boundary% p (isc:iec,jsc:jec), & - Ice_ocean_boundary% lrunoff (isc:iec,jsc:jec), & - Ice_ocean_boundary% frunoff (isc:iec,jsc:jec)) - - Ice_ocean_boundary%u_flux = 0.0 - Ice_ocean_boundary%v_flux = 0.0 - Ice_ocean_boundary%t_flux = 0.0 - Ice_ocean_boundary%q_flux = 0.0 - Ice_ocean_boundary%salt_flux = 0.0 - Ice_ocean_boundary%lw_flux = 0.0 - Ice_ocean_boundary%sw_flux_vis_dir = 0.0 - Ice_ocean_boundary%sw_flux_vis_dif = 0.0 - Ice_ocean_boundary%sw_flux_nir_dir = 0.0 - Ice_ocean_boundary%sw_flux_nir_dif = 0.0 - Ice_ocean_boundary%lprec = 0.0 - Ice_ocean_boundary%fprec = 0.0 - Ice_ocean_boundary%seaice_melt = 0.0 - Ice_ocean_boundary%seaice_melt_heat= 0.0 - Ice_ocean_boundary%mi = 0.0 - Ice_ocean_boundary%ice_fraction = 0.0 - Ice_ocean_boundary%u10_sqr = 0.0 - Ice_ocean_boundary%p = 0.0 - Ice_ocean_boundary%lrunoff = 0.0 - Ice_ocean_boundary%frunoff = 0.0 + allocate(Ice_ocean_boundary% u_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% v_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% t_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% q_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% salt_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% lw_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% sw_flux_vis_dir (isc:iec,jsc:jec), & + Ice_ocean_boundary% sw_flux_vis_dif (isc:iec,jsc:jec), & + Ice_ocean_boundary% sw_flux_nir_dir (isc:iec,jsc:jec), & + Ice_ocean_boundary% sw_flux_nir_dif (isc:iec,jsc:jec), & + Ice_ocean_boundary% lprec (isc:iec,jsc:jec), & + Ice_ocean_boundary% fprec (isc:iec,jsc:jec), & + Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),& + Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), & + Ice_ocean_boundary% mi (isc:iec,jsc:jec), & + Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec), & + Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec), & + Ice_ocean_boundary% p (isc:iec,jsc:jec), & + Ice_ocean_boundary% lrunoff (isc:iec,jsc:jec), & + Ice_ocean_boundary% frunoff (isc:iec,jsc:jec), & + source=0.0) + + ! Allocate memory for fields coming from multiple ice categories + if (Ice_ocean_boundary%ice_ncat > 0) & + allocate(Ice_ocean_boundary% afracr(isc:iec,jsc:jec), & + Ice_ocean_boundary% swnet_afracr(isc:iec,jsc:jec), & + Ice_ocean_boundary% swpen_ifrac_n(isc:iec,jsc:jec,1:Ice_ocean_boundary%ice_ncat), & + Ice_ocean_boundary% ifrac_n(isc:iec,jsc:jec,1:Ice_ocean_boundary%ice_ncat), & + source=0.0) if (cesm_coupled) then - allocate (Ice_ocean_boundary% hrain (isc:iec,jsc:jec), & - Ice_ocean_boundary% hsnow (isc:iec,jsc:jec), & - Ice_ocean_boundary% hrofl (isc:iec,jsc:jec), & - Ice_ocean_boundary% hrofi (isc:iec,jsc:jec), & - Ice_ocean_boundary% hevap (isc:iec,jsc:jec), & - Ice_ocean_boundary% hcond (isc:iec,jsc:jec)) - - Ice_ocean_boundary%hrain = 0.0 - Ice_ocean_boundary%hsnow = 0.0 - Ice_ocean_boundary%hrofl = 0.0 - Ice_ocean_boundary%hrofi = 0.0 - Ice_ocean_boundary%hevap = 0.0 - Ice_ocean_boundary%hcond = 0.0 + allocate(Ice_ocean_boundary% hrain (isc:iec,jsc:jec), & + Ice_ocean_boundary% hsnow (isc:iec,jsc:jec), & + Ice_ocean_boundary% hrofl (isc:iec,jsc:jec), & + Ice_ocean_boundary% hrofi (isc:iec,jsc:jec), & + Ice_ocean_boundary% hevap (isc:iec,jsc:jec), & + Ice_ocean_boundary% hcond (isc:iec,jsc:jec), & + source=0.0) + + ! Needed for MARBL + ! These are allocated separately to make it easier to pull out + ! of the cesm_coupled block if other models want to add BGC + allocate(Ice_ocean_boundary% nhx_dep (isc:iec,jsc:jec), & + Ice_ocean_boundary% noy_dep (isc:iec,jsc:jec), & + Ice_ocean_boundary% atm_fine_dust_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% atm_coarse_dust_flux (isc:iec,jsc:jec),& + Ice_ocean_boundary% seaice_dust_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% atm_bc_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% seaice_bc_flux (isc:iec,jsc:jec), & + Ice_ocean_boundary% atm_co2_prog (isc:iec,jsc:jec), & + Ice_ocean_boundary% atm_co2_diag (isc:iec,jsc:jec), & + source=0.0) endif call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method) if (use_waves) then if (wave_method == "EFACTOR") then - allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec) ) - Ice_ocean_boundary%lamult = 0.0 + allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec), source=0.0) else if (wave_method == "SURFACE_BANDS") then call query_ocean_state(ocean_state, NumWaveBands=Ice_ocean_boundary%num_stk_bands) - allocate(Ice_ocean_boundary%ustkb(isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), source=0.0) - allocate(Ice_ocean_boundary%vstkb(isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), source=0.0) - allocate(Ice_ocean_boundary%stk_wavenumbers(Ice_ocean_boundary%num_stk_bands), source=0.0) + allocate(Ice_ocean_boundary%ustkb(isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & + Ice_ocean_boundary%vstkb(isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & + Ice_ocean_boundary%stk_wavenumbers(Ice_ocean_boundary%num_stk_bands), & + source=0.0) call query_ocean_state(ocean_state, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, unscale=.true.) else call MOM_error(FATAL, "Unsupported WAVE_METHOD encountered in NUOPC cap.") @@ -776,6 +801,32 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofl" , "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_hrofi" , "will provide") + if (Ice_ocean_boundary%ice_ncat > 0) then + call fld_list_add(fldsToOcn_num, fldsToOcn, "Sf_afracr", "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_swnet_afracr", "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_swpen_ifrac_n", "will provide", & + ungridded_lbound=1, ungridded_ubound=Ice_ocean_boundary%ice_ncat) + call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac_n", "will provide", & + ungridded_lbound=1, ungridded_ubound=Ice_ocean_boundary%ice_ncat) + endif + + if (cesm_coupled) then + ! Fields needed for MARBL + call fld_list_add(fldsToOcn_num, fldsToOcn, "Faxa_ndep" , "will provide", & !-> nitrogen deposition + ungridded_lbound=1, ungridded_ubound=2) + call fld_list_add(fldsToOcn_num, fldsToOcn, "Faxa_dstwet" , "will provide", & + ungridded_lbound=1, ungridded_ubound=4) + call fld_list_add(fldsToOcn_num, fldsToOcn, "Faxa_dstdry" , "will provide", & + ungridded_lbound=1, ungridded_ubound=4) + call fld_list_add(fldsToOcn_num, fldsToOcn, "Faxa_bcph" , "will provide", & + ungridded_lbound=1, ungridded_ubound=3) + call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_flxdst" , "will provide") !-> ice runoff + call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_bcphi" , "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_bcpho" , "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "Sa_co2prog" , "will provide") !-> prognostic CO2 from atm + call fld_list_add(fldsToOcn_num, fldsToOcn, "Sa_co2diag" , "will provide") !-> diagnostic CO2 from atm + endif + if (use_waves) then if (wave_method == "EFACTOR") then call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide") @@ -799,6 +850,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_dhdy" , "will provide") call fld_list_add(fldsFrOcn_num, fldsFrOcn, "Fioo_q" , "will provide") call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_bldepth" , "will provide") + if (cesm_coupled) then + call fld_list_add(fldsFrOcn_num, fldsFrOcn, "Faoo_fco2_ocn", "will provide") + endif do n = 1,fldsToOcn_num call NUOPC_Advertise(importState, standardName=fldsToOcn(n)%stdname, name=fldsToOcn(n)%shortname, rc=rc) @@ -1142,7 +1196,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) "EPS_OMESH= ',i8,2(f21.13,3x),2(d21.5))" write(err_msg, frmt)n,lonMesh(n),lon(n), diff_lon, eps_omesh call MOM_error(FATAL, err_msg) - end if + endif diff_lat = abs(latMesh(n) - lat(n)) if (diff_lat > eps_omesh) then frmt = "('ERROR: Difference between ESMF Mesh and MOM6 domain coords is"//& @@ -1150,17 +1204,18 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) "EPS_OMESH= ',i8,2(f21.13,3x),2(d21.5))" write(err_msg, frmt)n,latMesh(n),lat(n), diff_lat, eps_omesh call MOM_error(FATAL, err_msg) - end if + endif if (abs(maskMesh(n) - mask(n)) > 0) then frmt = "('ERROR: ESMF mesh and MOM6 domain masks are inconsistent! - "//& "MOM n, maskMesh(n), mask(n) = ',3(i8,2x))" write(err_msg, frmt)n,maskMesh(n),mask(n) call MOM_error(FATAL, err_msg) - end if + endif end do ! realize the import and export fields using the mesh - call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", mesh=Emesh, rc=rc) + call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", & + ice_ocean_boundary=Ice_ocean_boundary, mesh=Emesh, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call MOM_RealizeFields(exportState, fldsFrOcn_num, fldsFrOcn, "Ocn export", mesh=Emesh, rc=rc) @@ -1176,10 +1231,9 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) 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 + allocate(mod2med_areacor(numOwnedElements), & + med2mod_areacor(numOwnedElements), & + source=1._ESMF_KIND_R8) #ifdef CESMCOUPLED ! Determine model areas and flux correction factors (module variables in mom_) @@ -1201,7 +1255,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth_L**2 mod2med_areacor(k) = model_areas(k) / mesh_areas(k) med2mod_areacor(k) = mesh_areas(k) / model_areas(k) - end if + endif end do end do deallocate(mesh_areas) @@ -1222,7 +1276,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) min_areacor_glob(1), max_areacor_glob(1), 'MOM6' write(stdout,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& min_areacor_glob(2), max_areacor_glob(2), 'MOM6' - end if + endif #endif deallocate(ownedElemCoords) @@ -1409,7 +1463,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) dataPtr_mask(i1,j1) = ocean_grid%mask2dT(ig,jg) dataPtr_xcen(i1,j1) = ocean_grid%geolonT(ig,jg) dataPtr_ycen(i1,j1) = ocean_grid%geolatT(ig,jg) - if(grid_attach_area) then + if (grid_attach_area) then dataPtr_area(i1,j1) = ocean_grid%US%L_to_m**2 * ocean_grid%areaT(ig,jg) endif enddo @@ -1451,7 +1505,8 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) gridOut = gridIn ! for now out same as in - call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", grid=gridIn, rc=rc) + call MOM_RealizeFields(importState, fldsToOcn_num, fldsToOcn, "Ocn import", & + ice_ocean_boundary=Ice_ocean_boundary, grid=gridIn, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call MOM_RealizeFields(exportState, fldsFrOcn_num, fldsFrOcn, "Ocn export", grid=gridOut, rc=rc) @@ -1736,7 +1791,7 @@ subroutine ModelAdvance(gcomp, rc) if (dbug > 0) then call state_diagnose(importState,subname//':IS ',rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + endif !--------------- ! Get ocean grid @@ -1755,10 +1810,10 @@ subroutine ModelAdvance(gcomp, rc) ! Update MOM6 !--------------- - if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM update_ocean_model: ") + if (profile_memory) call ESMF_VMLogMemInfo("Entering MOM update_ocean_model: ") call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled, & cesm_coupled) - if(profile_memory) call ESMF_VMLogMemInfo("Leaving MOM update_ocean_model: ") + if (profile_memory) call ESMF_VMLogMemInfo("Leaving MOM update_ocean_model: ") !--------------- ! Export Data @@ -1770,7 +1825,7 @@ subroutine ModelAdvance(gcomp, rc) if (dbug > 0) then call state_diagnose(exportState,subname//':ES ',rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if + endif endif !--------------- @@ -2025,7 +2080,7 @@ subroutine ModelSetRunClock(gcomp, rc) if (isPresent .and. isSet) then call ESMF_LogWrite(subname//" Restart_n = "//trim(cvalue), ESMF_LOGMSG_INFO) read(cvalue,*) restart_n - if (restart_n /= 0)then + if (restart_n /= 0) then call NUOPC_CompAttributeGet(gcomp, name="restart_option", value=cvalue, & isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -2068,7 +2123,7 @@ subroutine ModelSetRunClock(gcomp, rc) call ESMF_AlarmSet(restart_alarm, clock=mclock, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_LogWrite(subname//" Restart alarm is Created and Set", ESMF_LOGMSG_INFO) - end if + endif ! create a 1-shot alarm at the driver stop time stop_alarm = ESMF_AlarmCreate(mclock, ringtime=dstopTime, name = "stop_alarm", rc=rc) @@ -2176,9 +2231,9 @@ subroutine ocean_model_finalize(gcomp, rc) write_restart = .true. else write_restart = .false. - end if - if (write_restart)call ESMF_LogWrite("No Restart Alarm, writing restart at Finalize ", & - ESMF_LOGMSG_INFO) + endif + if (write_restart) call ESMF_LogWrite("No Restart Alarm, writing restart at Finalize ", & + ESMF_LOGMSG_INFO) call ocean_model_end(ocean_public, ocean_State, Time, write_restart=write_restart) @@ -2227,16 +2282,17 @@ subroutine State_SetScalar(value, scalar_id, State, mytask, scalar_name, scalar_ end subroutine State_SetScalar !> Realize the import and export fields using either a grid or a mesh. -subroutine MOM_RealizeFields(state, nfields, field_defs, tag, grid, mesh, rc) - type(ESMF_State) , intent(inout) :: state !< ESMF_State object for - !! import/export fields. - integer , intent(in) :: nfields !< Number of fields. - type(fld_list_type) , intent(inout) :: field_defs(:) !< Structure with field's - !! information. - character(len=*) , intent(in) :: tag !< Import or export. - type(ESMF_Grid) , intent(in), optional :: grid!< ESMF grid. - type(ESMF_Mesh) , intent(in), optional :: mesh!< ESMF mesh. - integer , intent(inout) :: rc !< Return code. +subroutine MOM_RealizeFields(state, nfields, field_defs, tag, ice_ocean_boundary, grid, mesh, rc) + type(ESMF_State) , intent(inout) :: state !< ESMF_State object for + !! import/export fields. + integer , intent(in) :: nfields !< Number of fields. + type(fld_list_type) , intent(inout) :: field_defs(:) !< Structure with field's + !! information. + type(ice_ocean_boundary_type), intent(inout), optional :: ice_ocean_boundary !< May need to nullify atm_co2 + character(len=*) , intent(in) :: tag !< Import or export. + type(ESMF_Grid) , intent(in) , optional :: grid!< ESMF grid. + type(ESMF_Mesh) , intent(in) , optional :: mesh!< ESMF mesh. + integer , intent(inout) :: rc !< Return code. ! local variables integer :: i @@ -2316,6 +2372,18 @@ subroutine MOM_RealizeFields(state, nfields, field_defs, tag, grid, mesh, rc) call ESMF_LogWrite(subname // tag // " Field "// trim(field_defs(i)%stdname) // " is not connected.", & ESMF_LOGMSG_INFO) + if (present(ice_ocean_boundary)) then + if (trim(field_defs(i)%stdname) == 'Sa_co2prog') then + if (is_root_pe()) write(stdout,*) subname // tag // " Nullifying ice_ocean_boundary%atm_co2_prog" + deallocate(ice_ocean_boundary%atm_co2_prog) + nullify(ice_ocean_boundary%atm_co2_prog) + elseif (trim(field_defs(i)%stdname) == 'Sa_co2diag') then + if (is_root_pe()) write(stdout,*) subname // tag // " Nullifying ice_ocean_boundary%atm_co2_diag" + deallocate(ice_ocean_boundary%atm_co2_diag) + nullify(ice_ocean_boundary%atm_co2_diag) + endif + endif + ! remove a not connected Field from State call ESMF_StateRemove(state, (/field_defs(i)%shortname/), rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -2390,7 +2458,7 @@ subroutine fld_list_add(num, fldlist, stdname, transferOffer, shortname, ungridd if (present(ungridded_lbound) .and. present(ungridded_ubound)) then fldlist(num)%ungridded_lbound = ungridded_lbound fldlist(num)%ungridded_ubound = ungridded_ubound - end if + endif end subroutine fld_list_add diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 125bae5748..d5ec9dc259 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -42,6 +42,7 @@ module MOM_cap_methods !> Get field pointer interface State_GetFldPtr module procedure State_GetFldPtr_1d + module procedure State_GetFldPtr_1d_from_2d module procedure State_GetFldPtr_2d end interface @@ -82,12 +83,14 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! Local Variables integer :: i, j, ib, ig, jg, n integer :: isc, iec, jsc, jec + integer :: esmf_ind integer :: nsc ! number of stokes drift components character(len=128) :: fldname real(ESMF_KIND_R8), allocatable :: taux(:,:) real(ESMF_KIND_R8), allocatable :: tauy(:,:) real(ESMF_KIND_R8), allocatable :: stkx(:,:,:) real(ESMF_KIND_R8), allocatable :: stky(:,:,:) + logical :: med_has_co2 character(len=*) , parameter :: subname = '(mom_import)' rc = ESMF_SUCCESS @@ -271,6 +274,159 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, if (ChkErr(rc,__LINE__,u_FILE_u)) return endif + !--------------! + ! MARBL fields ! + !--------------! + + ! seaice_dust_flux, nhx_dep, and noy_dep are single fields from the coupler + ! atm_fine_dust_flux, atm_coarse_dust_flux, atm_bc_flux, and seaice_bc_flux + ! are all sums of multiple fields and will be treated slightly differently + ! For those fields, we use do_sum = .true. + + !---- + ! nhx deposition + !---- + if (associated(ice_ocean_boundary%nhx_dep)) then + call state_getimport(importState, 'Faxa_ndep', & + isc, iec, jsc, jec, ice_ocean_boundary%nhx_dep(:,:), & + areacor=med2mod_areacor, esmf_ind=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + + !---- + ! noy deposition + !---- + if (associated(ice_ocean_boundary%noy_dep)) then + call state_getimport(importState, 'Faxa_ndep', & + isc, iec, jsc, jec, ice_ocean_boundary%noy_dep(:,:), & + areacor=med2mod_areacor, esmf_ind=2, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + + !---- + ! atmospheric CO2 concentration + ! might not be passed from atmosphere component, + ! in which the pointer(s) will not be associated + !---- + if (associated(ice_ocean_boundary%atm_co2_prog)) then + call state_getimport(importState, 'Sa_co2prog', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_co2_prog(:,:), rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + if (associated(ice_ocean_boundary%atm_co2_diag)) then + call state_getimport(importState, 'Sa_co2diag', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_co2_diag(:,:), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + endif + + !---- + ! fine dust flux from atmosphere + !---- + if (associated(ice_ocean_boundary%atm_fine_dust_flux)) then + ice_ocean_boundary%atm_fine_dust_flux(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Faxa_dstwet', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_fine_dust_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., esmf_ind=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getimport(importState, 'Faxa_dstdry', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_fine_dust_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., esmf_ind=1, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + + !---- + ! coarse dust flux from atmosphere + !---- + if (associated(ice_ocean_boundary%atm_coarse_dust_flux)) then + ice_ocean_boundary%atm_coarse_dust_flux(:,:) = 0._ESMF_KIND_R8 + do esmf_ind=2,4 + call state_getimport(importState, 'Faxa_dstwet', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_coarse_dust_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., esmf_ind=esmf_ind, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getimport(importState, 'Faxa_dstdry', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_coarse_dust_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., esmf_ind=esmf_ind, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + enddo + end if + + !---- + ! dust flux from sea ice + !---- + if (associated(ice_ocean_boundary%seaice_dust_flux)) then + call state_getimport(importState, 'Fioi_flxdst', & + isc, iec, jsc, jec, ice_ocean_boundary%seaice_dust_flux, & + areacor=med2mod_areacor, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + + !---- + ! black carbon flux from atmosphere + !---- + if (associated(ice_ocean_boundary%atm_bc_flux)) then + ice_ocean_boundary%atm_bc_flux(:,:) = 0._ESMF_KIND_R8 + do esmf_ind=1,3 + call state_getimport(importState, 'Faxa_bcph', & + isc, iec, jsc, jec, ice_ocean_boundary%atm_bc_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., esmf_ind=esmf_ind, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + enddo + endif + + !---- + ! black carbon flux from sea ice + !---- + if (associated(ice_ocean_boundary%seaice_bc_flux)) then + ice_ocean_boundary%seaice_bc_flux(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Fioi_bcpho', & + isc, iec, jsc, jec, ice_ocean_boundary%seaice_bc_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getimport(importState, 'Fioi_bcphi', & + isc, iec, jsc, jec, ice_ocean_boundary%seaice_bc_flux(:,:), & + areacor=med2mod_areacor, do_sum=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + + ! Fields coming from coupler per ice category + if (ice_ocean_boundary%ice_ncat > 0) then + call state_getimport(importState, 'Sf_afracr', & + isc, iec, jsc, jec, ice_ocean_boundary%afracr(:,:), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + call state_getimport(importState, 'Foxx_swnet_afracr', & + isc, iec, jsc, jec, ice_ocean_boundary%swnet_afracr(:,:), & + 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, 'Fioi_swpen_ifrac_n', & + isc, iec, jsc, jec, 1, ice_ocean_boundary%ice_ncat, & + ice_ocean_boundary%swpen_ifrac_n(:,:,:), & + 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, 'Si_ifrac_n', & + isc, iec, jsc, jec, 1, ice_ocean_boundary%ice_ncat, & + ice_ocean_boundary%ifrac_n(:,:,:), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + endif ! multiple ice categories + !---- ! salt flux from ice !---- @@ -529,16 +685,13 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ! Sea-surface zonal and meridional slopes !---------------- - allocate(ssh(ocean_grid%isd:ocean_grid%ied,ocean_grid%jsd:ocean_grid%jed)) ! local indices with halos - allocate(dhdx(isc:iec, jsc:jec)) !global indices without halos - allocate(dhdy(isc:iec, jsc:jec)) !global indices without halos + allocate(ssh(ocean_grid%isd:ocean_grid%ied,ocean_grid%jsd:ocean_grid%jed), & ! local indices with halos + dhdx(isc:iec, jsc:jec), & !global indices without halos + dhdy(isc:iec, jsc:jec), & !global indices without halos + source=0.0_ESMF_KIND_R8) allocate(dhdx_rot(isc:iec, jsc:jec)) !global indices without halos allocate(dhdy_rot(isc:iec, jsc:jec)) !global indices without halos - ssh = 0.0_ESMF_KIND_R8 - dhdx = 0.0_ESMF_KIND_R8 - dhdy = 0.0_ESMF_KIND_R8 - ! Make a copy of ssh in order to do a halo update (ssh has local indexing with halos) do j = ocean_grid%jsc, ocean_grid%jec jloc = j + ocean_grid%jdg_offset @@ -629,6 +782,16 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, call State_SetExport(exportState, 'So_dhdy', isc, iec, jsc, jec, dhdy_rot, ocean_grid, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ------- + ! CO2 Flux + ! ------- + call ESMF_StateGet(exportState, 'Faoo_fco2_ocn', itemFlag, rc=rc) + if (itemFlag /= ESMF_STATEITEM_NOTFOUND) then + call State_SetExport(exportState, 'Faoo_fco2_ocn', isc, iec, jsc, jec, & + ocean_public%fco2_ocn, ocean_grid, areacor=mod2med_areacor, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif + deallocate(ssh, dhdx, dhdy, dhdx_rot, dhdy_rot) end subroutine mom_export @@ -654,6 +817,32 @@ subroutine State_GetFldPtr_1d(State, fldname, fldptr, rc) end subroutine State_GetFldPtr_1d +!> Get specific 1D field pointer from 2D field +subroutine State_GetFldPtr_1d_from_2d(State, fldname, esmf_ind, fldptr, rc) + type(ESMF_State) , intent(in) :: State !< ESMF state + character(len=*) , intent(in) :: fldname !< Field name + real(ESMF_KIND_R8), pointer :: fldptr(:)!< Pointer to the 1D field + integer, intent(in) :: esmf_ind !< Index into 2D ESMF array + integer, optional , intent(out) :: rc !< Return code + + ! local variables + real(ESMF_KIND_R8), pointer :: fldptr2d(:,:)!< Pointer to the 1D field + type(ESMF_Field) :: lfield + integer :: lrc + character(len=*),parameter :: subname='(MOM_cap:State_GetFldPtr)' + + call ESMF_StateGet(State, itemName=trim(fldname), field=lfield, rc=lrc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldGet(lfield, farrayPtr=fldptr2d, rc=lrc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (.not. associated(fldptr)) allocate(fldptr(size(fldptr2d,2))) + fldptr = fldptr2d(esmf_ind, :) + + if (present(rc)) rc = lrc + +end subroutine State_GetFldPtr_1d_from_2d + !> Get field pointer 2D subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) type(ESMF_State) , intent(in) :: State !< ESMF state @@ -676,7 +865,7 @@ subroutine State_GetFldPtr_2d(State, fldname, fldptr, rc) end subroutine State_GetFldPtr_2d !> Map 2d import state field to output array -subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum, areacor, rc) +subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum, areacor, esmf_ind, 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 @@ -691,18 +880,25 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum 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, optional, intent(in) :: esmf_ind integer , intent(out) :: rc !< Return code ! local variables type(ESMF_StateItem_Flag) :: itemFlag integer :: n, i, j, i1, j1 integer :: lbnd1,lbnd2 + logical :: do_sum_loc real(ESMF_KIND_R8), pointer :: dataPtr1d(:) real(ESMF_KIND_R8), pointer :: dataPtr2d(:,:) character(len=*) , parameter :: subname='(MOM_cap_methods:state_getimport_2d)' ! ---------------------------------------------- rc = ESMF_SUCCESS + if (present(do_sum)) then + do_sum_loc = do_sum + else + do_sum_loc = .false. + endif call ESMF_StateGet(State, trim(fldname), itemFlag, rc=rc) if (itemFlag /= ESMF_STATEITEM_NOTFOUND) then @@ -710,7 +906,11 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum if (geomtype == ESMF_GEOMTYPE_MESH) then ! get field pointer - call state_getfldptr(state, trim(fldname), dataptr1d, rc) + if (present(esmf_ind)) then + call state_getfldptr(state, trim(fldname), esmf_ind, dataptr1d, rc) + else + call state_getfldptr(state, trim(fldname), dataptr1d, rc) + endif if (ChkErr(rc,__LINE__,u_FILE_u)) return ! determine output array and apply area correction if present @@ -718,23 +918,23 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum do j = jsc,jec do i = isc,iec n = n + 1 - if (present(do_sum)) then + if (do_sum_loc) then 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 + endif else if (present(areacor)) then output(i,j) = dataPtr1d(n) * areacor(n) else output(i,j) = dataPtr1d(n) - end if + endif endif enddo enddo - else if (geomtype == ESMF_GEOMTYPE_GRID) then + elseif (geomtype == ESMF_GEOMTYPE_GRID) then call state_getfldptr(state, trim(fldname), dataptr2d, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -746,7 +946,7 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum j1 = j + lbnd2 - jsc do i = isc, iec i1 = i + lbnd1 - isc - if (present(do_sum)) then + if (do_sum_loc) then output(i,j) = output(i,j) + dataPtr2d(i1,j1) else output(i,j) = dataPtr2d(i1,j1) @@ -784,11 +984,17 @@ subroutine State_GetImport_3d(state, fldname, isc, iec, jsc, jec, lbd, ubd, outp type(ESMF_StateItem_Flag) :: itemFlag integer :: n, i, j, i1, j1, u integer :: lbnd1,lbnd2 + logical :: do_sum_loc real(ESMF_KIND_R8), pointer :: dataPtr2d(:,:) character(len=*) , parameter :: subname='(MOM_cap_methods:state_getimport_3d)' ! ---------------------------------------------- rc = ESMF_SUCCESS + if (present(do_sum)) then + do_sum_loc = do_sum + else + do_sum_loc = .false. + endif call ESMF_StateGet(State, trim(fldname), itemFlag, rc=rc) if (itemFlag /= ESMF_STATEITEM_NOTFOUND) then @@ -805,18 +1011,18 @@ subroutine State_GetImport_3d(state, fldname, isc, iec, jsc, jec, lbd, ubd, outp do j = jsc,jec do i = isc,iec n = n + 1 - if (present(do_sum)) then + if (do_sum_loc) then if (present(areacor)) then output(i,j,u) = output(i,j,u) + dataPtr2d(u,n) * areacor(n) else output(i,j,u) = output(i,j,u) + dataPtr2d(u,n) - end if + endif else if (present(areacor)) then output(i,j,u) = dataPtr2d(u,n) * areacor(n) else output(i,j,u) = dataPtr2d(u,n) - end if + endif endif enddo enddo @@ -887,7 +1093,7 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid do n = 1,(size(dataPtr1d)) dataPtr1d(n) = dataPtr1d(n) * areacor(n) enddo - end if + endif ! if a maskmap is provided, set exports of all eliminated cells to zero. if (associated(ocean_grid%Domain%maskmap)) then diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 9ac40daaa4..329f436e48 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -108,17 +108,18 @@ module MOM_ocean_model_nuopc !! a global max across ocean and non-ocean processors can be !! used to determine its value. real, pointer, dimension(:,:) :: & - t_surf => NULL(), & !< SST on t-cell (degrees Kelvin) - s_surf => NULL(), & !< SSS on t-cell (psu) - u_surf => NULL(), & !< i-velocity at the locations indicated by stagger, m/s. - v_surf => NULL(), & !< j-velocity at the locations indicated by stagger, m/s. + t_surf => NULL(), & !< SST on t-cell (degrees Kelvin) + s_surf => NULL(), & !< SSS on t-cell (psu) + u_surf => NULL(), & !< i-velocity at the locations indicated by stagger, m/s. + v_surf => NULL(), & !< j-velocity at the locations indicated by stagger, m/s. sea_lev => NULL(), & !< Sea level in m after correction for surface pressure, - !! i.e. dzt(1) + eta_t + patm/rho0/grav (m) - frazil =>NULL(), & !< Accumulated heating (in Joules/m^2) from frazil - !! formation in the ocean. + !! i.e. dzt(1) + eta_t + patm/rho0/grav (m) + frazil =>NULL(), & !< Accumulated heating (in Joules/m^2) from frazil + !! formation in the ocean. melt_potential => NULL(), & !< Instantaneous heat used to melt sea ice (in J/m^2) - area => NULL(), & !< cell area of the ocean surface, in m2. - OBLD => NULL() !< Ocean boundary layer depth, in m. + area => NULL(), & !< cell area of the ocean surface, in m2. + OBLD => NULL(), & !< Ocean boundary layer depth, in m. + fco2_ocn => NULL() !< Ocean CO2 flux, in kg CO2/m^2/s type(coupler_2d_bc_type) :: fields !< A structure that may contain named !! arrays of tracer-related surface fields. integer :: avg_kount !< A count of contributions to running @@ -255,6 +256,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! min(HFrz, OBLD), where OBLD is the boundary layer depth. !! If HFrz <= 0 (default), melt potential will not be computed. logical :: use_melt_pot !< If true, allocate melt_potential array + logical :: use_MARBL !< If true, allocate surface co2 array ! This include declares and sets the variable "version". @@ -378,12 +380,14 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, & "If true, enables surface wave modules.", default=.false.) + call get_param(param_file, mdl, "USE_MARBL_TRACERS", use_MARBL, & + default=.false., do_not_log=.true.) ! Consider using a run-time flag to determine whether to do the diagnostic ! vertical integrals, since the related 3-d sums are not negligible in cost. call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, & do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & - use_meltpot=use_melt_pot) + use_meltpot=use_melt_pot, use_marbl_tracers=use_MARBL) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves) @@ -538,6 +542,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%grid, OS%US, OS%forcing_CSp) if (OS%fluxes%fluxes_used) then + + ! enable_averages() is necessary to post forcing fields to diagnostics + call enable_averages(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag) + if (do_thermo) & call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, dt_coupling, & OS%grid, OS%US, OS%forcing_CSp, OS%sfc_state, & @@ -781,7 +789,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time, write_restart) type(time_type), intent(in) :: Time !< The model time, used for writing restarts. logical, intent(in) :: write_restart !< true => write restart file - if(write_restart)call ocean_model_save_restart(Ocean_state, Time) + if (write_restart) call ocean_model_save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%diag, end_diag_manager=.true.) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) @@ -853,25 +861,19 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, endif call mpp_get_compute_domain(Ocean_sfc%Domain, isc, iec, jsc, jec) - allocate ( Ocean_sfc%t_surf (isc:iec,jsc:jec), & - Ocean_sfc%s_surf (isc:iec,jsc:jec), & - Ocean_sfc%u_surf (isc:iec,jsc:jec), & - Ocean_sfc%v_surf (isc:iec,jsc:jec), & - Ocean_sfc%sea_lev(isc:iec,jsc:jec), & - Ocean_sfc%area (isc:iec,jsc:jec), & - Ocean_sfc%OBLD (isc:iec,jsc:jec), & - Ocean_sfc%melt_potential(isc:iec,jsc:jec), & - Ocean_sfc%frazil (isc:iec,jsc:jec)) - - Ocean_sfc%t_surf = 0.0 ! time averaged sst (Kelvin) passed to atmosphere/ice model - Ocean_sfc%s_surf = 0.0 ! time averaged sss (psu) passed to atmosphere/ice models - Ocean_sfc%u_surf = 0.0 ! time averaged u-current (m/sec) passed to atmosphere/ice models - Ocean_sfc%v_surf = 0.0 ! time averaged v-current (m/sec) passed to atmosphere/ice models - Ocean_sfc%sea_lev = 0.0 ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav - Ocean_sfc%frazil = 0.0 ! time accumulated frazil (J/m^2) passed to ice model - Ocean_sfc%melt_potential = 0.0 ! time accumulated melt potential (J/m^2) passed to ice model - Ocean_sfc%OBLD = 0.0 ! ocean boundary layer depth, in m - Ocean_sfc%area = 0.0 + allocate(Ocean_sfc%t_surf (isc:iec,jsc:jec), & ! time averaged sst (Kelvin) passed to atmosphere/ice model + Ocean_sfc%s_surf (isc:iec,jsc:jec), & ! time averaged sss (psu) passed to atmosphere/ice models + Ocean_sfc%u_surf (isc:iec,jsc:jec), & ! time averaged u-current (m/sec) passed to atmosphere/ice models + Ocean_sfc%v_surf (isc:iec,jsc:jec), & ! time averaged v-current (m/sec) passed to atmosphere/ice models + Ocean_sfc%sea_lev(isc:iec,jsc:jec), & ! time averaged thickness of top model grid cell (m) plus + ! patm/rho0/grav + Ocean_sfc%frazil (isc:iec,jsc:jec), & ! time accumulated frazil (J/m^2) passed to ice model + Ocean_sfc%melt_potential(isc:iec,jsc:jec), & ! time accumulated melt potential (J/m^2) passed to ice model + Ocean_sfc%area (isc:iec,jsc:jec), & + Ocean_sfc%OBLD (isc:iec,jsc:jec), & ! ocean boundary layer depth, in m + Ocean_sfc%fco2_ocn(isc:iec,jsc:jec), & ! time averaged co2 flux (kg/m^2/s) passed to atmosphere model + source=0.0) + Ocean_sfc%axes = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics if (present(gas_fields_ocn)) then @@ -968,6 +970,12 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ enddo ; enddo endif + if (allocated(sfc_state%fco2)) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%fco2_ocn(i,j) = US%RZ_T_to_kg_m2s * sfc_state%fco2(i+i0,j+j0) + enddo ; enddo + endif + if (Ocean_sfc%stagger == AGRID) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index d699697140..3e8f80e265 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -38,6 +38,8 @@ module MOM_surface_forcing_nuopc use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init use user_revise_forcing, only : user_revise_forcing_CS use iso_fortran_env, only : int64 +use MARBL_forcing_mod, only : marbl_forcing_CS, MARBL_forcing_init +use MARBL_forcing_mod, only : convert_driver_fields_to_forcings implicit none ; private @@ -79,6 +81,7 @@ module MOM_surface_forcing_nuopc !! pressure limited by max_p_surf instead of the !! full atmospheric pressure. The default is true. logical :: use_CFC !< enables the MOM_CFC_cap tracer package. + logical :: use_marbl_tracers !< enables the MARBL tracer package. logical :: enthalpy_cpl !< Controls if enthalpy terms are provided by the coupler or computed !! internally. real :: gust_const !< constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa] @@ -152,6 +155,8 @@ module MOM_surface_forcing_nuopc type(MOM_restart_CS), pointer :: restart_CSp => NULL() type(user_revise_forcing_CS), pointer :: urf_CS => NULL() + + type(marbl_forcing_CS), pointer :: marbl_forcing_CSp => NULL() !< parameters for getting MARBL forcing end type surface_forcing_CS !> Structure corresponding to forcing, but with the elements, units, and conventions @@ -186,6 +191,19 @@ module MOM_surface_forcing_nuopc !< on ocean surface [Pa] real, pointer, dimension(:,:) :: ice_fraction =>NULL() !< fractional ice area [nondim] real, pointer, dimension(:,:) :: u10_sqr =>NULL() !< wind speed squared at 10m [m2/s2] + real, pointer, dimension(:,:) :: nhx_dep =>NULL() !< Nitrogen deposition [kg/m^2/s] + real, pointer, dimension(:,:) :: noy_dep =>NULL() !< Nitrogen deposition [kg/m^2/s] + real, pointer, dimension(:,:) :: atm_co2_prog =>NULL() !< Prognostic atmospheric co2 concentration [ppm] + real, pointer, dimension(:,:) :: atm_co2_diag =>NULL() !< Diagnostic atmospheric co2 concentration [ppm] + real, pointer, dimension(:,:) :: atm_fine_dust_flux =>NULL() !< Fine dust flux from atmosphere [kg/m^2/s] + real, pointer, dimension(:,:) :: atm_coarse_dust_flux =>NULL() !< Coarse dust flux from atmosphere [kg/m^2/s] + real, pointer, dimension(:,:) :: seaice_dust_flux =>NULL() !< Dust flux from seaice [kg/m^2/s] + real, pointer, dimension(:,:) :: atm_bc_flux =>NULL() !< Black carbon flux from atmosphere [kg/m^2/s] + real, pointer, dimension(:,:) :: seaice_bc_flux =>NULL() !< Black carbon flux from seaice [kg/m^2/s] + real, pointer, dimension(:,:) :: afracr =>NULL() + real, pointer, dimension(:,:) :: swnet_afracr =>NULL() + real, pointer, dimension(:,:,:) :: swpen_ifrac_n =>NULL() + real, pointer, dimension(:,:,:) :: ifrac_n =>NULL() real, pointer, dimension(:,:) :: mi =>NULL() !< mass of ice [kg/m2] real, pointer, dimension(:,:) :: ice_rigidity =>NULL() !< rigidity of the sea ice, sea-ice and !! ice-shelves, expressed as a coefficient @@ -208,6 +226,10 @@ module MOM_surface_forcing_nuopc !! flux-exchange code, based on what the sea-ice !! model is providing. Otherwise, the value from !! the surface_forcing_CS is used. + + ! Forcing when receiving multiple ice categories from CMEPS + integer :: ice_ncat !< Number of ice categories coming from coupler + !! (1 => not using separate categories) end type ice_ocean_boundary_type integer :: id_clock_forcing @@ -297,8 +319,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (fluxes%dt_buoy_accum < 0) then call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., & press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, & - cfc=CS%use_CFC, hevap=CS%enthalpy_cpl, tau_mag=.true.) - !call safe_alloc_ptr(fluxes%omega_w2x,isd,ied,jsd,jed) + cfc=CS%use_CFC, marbl=CS%use_marbl_tracers, hevap=CS%enthalpy_cpl, & + tau_mag=.true., ice_ncat=IOB%ice_ncat) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) @@ -561,6 +583,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo + ! Copy MARBL-specific IOB fields into fluxes; also set some MARBL-specific forcings to other values + ! (constants, values from netCDF, etc) + call convert_driver_fields_to_forcings(IOB%atm_fine_dust_flux, IOB%atm_coarse_dust_flux, & + IOB%seaice_dust_flux, IOB%atm_bc_flux, IOB%seaice_bc_flux, & + IOB%nhx_dep, IOB%noy_dep, IOB%atm_co2_prog, IOB%atm_co2_diag, & + IOB%afracr, IOB%swnet_afracr, IOB%ifrac_n, IOB%swpen_ifrac_n, & + Time, G, US, i0, j0, fluxes, CS%marbl_forcing_CSp) + ! wave to ocean coupling if ( associated(IOB%lamult)) then do j=js,je; do i=is,ie @@ -1209,6 +1239,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, & default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_marbl_tracers, & + default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "ENTHALPY_FROM_COUPLER", CS%enthalpy_cpl, & "If True, the heat (enthalpy) associated with mass entering/leaving the "//& "ocean is provided via coupler.", default=.false.) @@ -1388,6 +1421,10 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call data_override_init(Ocean_domain_in=G%Domain%mpp_domain) endif + ! Set up MARBL forcing control structure + call MARBL_forcing_init(G, US, param_file, diag, Time, CS%inputdir, CS%use_marbl_tracers, & + CS%marbl_forcing_CSp) + if (present(restore_salt)) then ; if (restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) @@ -1496,6 +1533,60 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks endif + ! MARBL forcing + if (associated(iobt%atm_fine_dust_flux)) then + chks = field_chksum(iobt%atm_fine_dust_flux) + if (root) write(outunit,110) 'iobt%atm_fine_dust_flux ', chks + endif + if (associated(iobt%atm_coarse_dust_flux)) then + chks = field_chksum(iobt%atm_coarse_dust_flux) + if (root) write(outunit,110) 'iobt%atm_coarse_dust_flux ', chks + endif + if (associated(iobt%seaice_dust_flux)) then + chks = field_chksum(iobt%seaice_dust_flux) + if (root) write(outunit,110) 'iobt%seaice_dust_flux ', chks + endif + if (associated(iobt%atm_bc_flux)) then + chks = field_chksum(iobt%atm_bc_flux) + if (root) write(outunit,110) 'iobt%atm_bc_flux ', chks + endif + if (associated(iobt%seaice_bc_flux)) then + chks = field_chksum(iobt%seaice_bc_flux) + if (root) write(outunit,110) 'iobt%seaice_bc_flux ', chks + endif + if (associated(iobt%nhx_dep)) then + chks = field_chksum(iobt%nhx_dep) + if (root) write(outunit,100) 'iobt%nhx_dep ', chks + endif + if (associated(iobt%noy_dep)) then + chks = field_chksum(iobt%noy_dep) + if (root) write(outunit,100) 'iobt%noy_dep ', chks + endif + if (associated(iobt%atm_co2_prog)) then + chks = field_chksum(iobt%atm_co2_prog) + if (root) write(outunit,110) 'iobt%atm_co2_prog ', chks + endif + if (associated(iobt%atm_co2_diag)) then + chks = field_chksum(iobt%atm_co2_diag) + if (root) write(outunit,110) 'iobt%atm_co2_diag ', chks + endif + if (associated(iobt%afracr)) then + chks = field_chksum(iobt%afracr) + if (root) write(outunit,100) 'iobt%afracr ', chks + endif + if (associated(iobt%swnet_afracr)) then + chks = field_chksum(iobt%swnet_afracr) + if (root) write(outunit,110) 'iobt%swnet_afracr ', chks + endif + if (associated(iobt%ifrac_n)) then + chks = field_chksum(iobt%ifrac_n) + if (root) write(outunit,100) 'iobt%ifrac_n ', chks + endif + if (associated(iobt%swpen_ifrac_n)) then + chks = field_chksum(iobt%swpen_ifrac_n) + if (root) write(outunit,110) 'iobt%swpen_ifrac_n ', chks + endif + ! enthalpy if (associated(iobt%hrofl)) then chks = field_chksum( iobt%hrofl ) ; if (root) write(outunit,100) 'iobt%hrofl ', chks @@ -1517,6 +1608,7 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) endif 100 FORMAT(" CHECKSUM::",A20," = ",Z20) +110 FORMAT(" CHECKSUM::",A30," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') diff --git a/config_src/drivers/nuopc_cap/time_utils.F90 b/config_src/drivers/nuopc_cap/time_utils.F90 index 81efcd2765..46f922d5bf 100644 --- a/config_src/drivers/nuopc_cap/time_utils.F90 +++ b/config_src/drivers/nuopc_cap/time_utils.F90 @@ -130,7 +130,7 @@ function fms2esmf_time(time, calkind) integer :: rc - if(present(calkind)) then + if (present(calkind)) then l_calkind = calkind else l_calkind = fms2esmf_cal(fms_get_calendar_type()) @@ -154,7 +154,7 @@ function string_to_date(string, rc) ! Local variables integer :: yr,mon,day,hr,min,sec - if(present(rc)) rc = ESMF_SUCCESS + if (present(rc)) rc = ESMF_SUCCESS read(string, '(I4.4,I2.2,I2.2,".",I2.2,I2.2,I2.2)') yr, mon, day, hr, min, sec string_to_date = set_date(yr, mon, day, hr, min, sec) diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 index d17db5a9a1..3a8303e561 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -56,6 +56,8 @@ module MOM_surface_forcing use BFB_surface_forcing, only : BFB_surface_forcing_init, BFB_surface_forcing_CS use dumbbell_surface_forcing, only : dumbbell_surface_forcing_init, dumbbell_surface_forcing_CS use dumbbell_surface_forcing, only : dumbbell_buoyancy_forcing +use MARBL_forcing_mod, only : marbl_forcing_CS, MARBL_forcing_init +use MARBL_forcing_mod, only : convert_driver_fields_to_forcings implicit none ; private @@ -116,6 +118,7 @@ module MOM_surface_forcing !! Dates before 20190101 use original answers. !! Dates after 20190101 use a form of the gyre wind stresses that are !! rotationally invariant and more likely to be the same between compilers. + logical :: use_marbl_tracers !< If true, allocate memory for forcing needed by MARBL logical :: fix_ustar_gustless_bug !< If true correct a bug in the time-averaging of the !! gustless wind friction velocity. ! if WIND_CONFIG=='scurves' then use the following to define a piecewise scurve profile @@ -216,6 +219,7 @@ module MOM_surface_forcing type(MESO_surface_forcing_CS), pointer :: MESO_forcing_CSp => NULL() type(idealized_hurricane_CS), pointer :: idealized_hurricane_CSp => NULL() type(SCM_CVmix_tests_CS), pointer :: SCM_CVmix_tests_CSp => NULL() + type(marbl_forcing_CS), pointer :: marbl_forcing_CSp => NULL() !>@} end type surface_forcing_CS @@ -255,7 +259,7 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US ! Allocate memory for the mechanical and thermodynamic forcing fields. call allocate_mech_forcing(G, forces, stress=.true., ustar=.not.CS%nonBous, press=.true., tau_mag=CS%nonBous) - call allocate_forcing_type(G, fluxes, ustar=.not.CS%nonBous, tau_mag=CS%nonBous, & + call allocate_forcing_type(G, fluxes, ustar=.not.CS%nonBous, marbl=CS%use_marbl_tracers, tau_mag=CS%nonBous, & fix_accum_bug=CS%fix_ustar_gustless_bug) if (trim(CS%buoy_config) /= "NONE") then if ( CS%use_temperature ) then @@ -351,6 +355,10 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US endif endif + if (CS%use_marbl_tracers) then + call MARBL_forcing_from_data_override(fluxes, day_center, G, US, CS) + endif + if (associated(CS%tracer_flow_CSp)) then call call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, CS%Rho0, & CS%tracer_flow_CSp) @@ -1542,6 +1550,94 @@ subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, US, CS) call callTree_leave("buoyancy_forcing_linear") end subroutine buoyancy_forcing_linear + +! Sets the necessary MARBL forcings via the data override facility. +subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS) + type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields + type(time_type), intent(in) :: day !< The time of the fluxes + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(surface_forcing_CS), pointer :: CS !< pointer to control structure returned by + !! a previous surface_forcing_init call + ! Local variables + real, pointer, dimension(:,:) :: atm_co2_prog =>NULL() !< Prognostic atmospheric CO2 concentration [ppm] + real, pointer, dimension(:,:) :: atm_co2_diag =>NULL() !< Diagnostic atmospheric CO2 concentration [ppm] + real, pointer, dimension(:,:) :: atm_fine_dust_flux =>NULL() !< Fine dust flux from atmosphere [kg/m^2/s ~> RZ/T] + real, pointer, dimension(:,:) :: atm_coarse_dust_flux =>NULL() !< Coarse dust flux from atmosphere [kg/m^2/s ~> RZ/T] + real, pointer, dimension(:,:) :: seaice_dust_flux =>NULL() !< Dust flux from seaice [kg/m^2/s ~> RZ/T] + real, pointer, dimension(:,:) :: atm_bc_flux =>NULL() !< Black carbon flux from atmosphere [kg/m^2/s ~> RZ/T] + real, pointer, dimension(:,:) :: seaice_bc_flux =>NULL() !< Black carbon flux from seaice [kg/m^2/s ~> RZ/T] + real, pointer, dimension(:,:) :: nhx_dep =>NULL() !< Nitrogen deposition [kg/m^2/s ~> RZ/T] + real, pointer, dimension(:,:) :: noy_dep =>NULL() !< Nitrogen deposition [kg/m^2/s ~> RZ/T] + integer :: isc, iec, jsc, jec + + ! Necessary null pointers for arguments to convert_driver_fields_to_forcings() + ! Since they are null, MARBL will not use multiple ice categories + real, pointer, dimension(:,:) :: afracr =>NULL() + real, pointer, dimension(:,:) :: swnet_afracr =>NULL() + real, pointer, dimension(:,:,:) :: swpen_ifrac_n =>NULL() + real, pointer, dimension(:,:,:) :: ifrac_n =>NULL() + + call callTree_enter("MARBL_forcing_from_data_override, MOM_surface_forcing.F90") + + if (.not.CS%dataOverrideIsInitialized) then + call data_override_init(G%Domain) + CS%dataOverrideIsInitialized = .True. + endif + + ! Allocate memory for pointers + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec + allocate ( atm_co2_prog (isc:iec,jsc:jec), & + atm_co2_diag (isc:iec,jsc:jec), & + atm_fine_dust_flux (isc:iec,jsc:jec), & + atm_coarse_dust_flux (isc:iec,jsc:jec), & + seaice_dust_flux (isc:iec,jsc:jec), & + atm_bc_flux (isc:iec,jsc:jec), & + seaice_bc_flux (isc:iec,jsc:jec), & + nhx_dep (isc:iec,jsc:jec), & + noy_dep (isc:iec,jsc:jec), & + source=0.0) + + + ! fluxes used directly as MARBL inputs + ! (should be scaled) + call data_override(G%Domain, 'ice_fraction', fluxes%ice_fraction, day) + call data_override(G%Domain, 'u10_sqr', fluxes%u10_sqr, day, scale=US%m_s_to_L_T**2) + + ! fluxes used to compute MARBL inputs + ! These are kept in physical units, and will be scaled appropriately in + ! convert_driver_fields_to_forcings() + call data_override(G%Domain, 'atm_co2_prog', atm_co2_prog, day) + call data_override(G%Domain, 'atm_co2_diag', atm_co2_diag, day) + call data_override(G%Domain, 'atm_fine_dust_flux', atm_fine_dust_flux, day) + call data_override(G%Domain, 'atm_coarse_dust_flux', atm_coarse_dust_flux, day) + call data_override(G%Domain, 'atm_bc_flux', atm_bc_flux, day) + call data_override(G%Domain, 'seaice_dust_flux', seaice_dust_flux, day) + call data_override(G%Domain, 'seaice_bc_flux', seaice_bc_flux, day) + call data_override(G%Domain, 'nhx_dep', nhx_dep, day) + call data_override(G%Domain, 'noy_dep', noy_dep, day) + + call convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, & + seaice_dust_flux, atm_bc_flux, seaice_bc_flux, & + nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, & + afracr, swnet_afracr, ifrac_n, swpen_ifrac_n, & + day, G, US, 0, 0, fluxes, CS%marbl_forcing_CSp) + + deallocate ( atm_co2_prog, & + atm_co2_diag, & + atm_fine_dust_flux, & + atm_coarse_dust_flux, & + seaice_dust_flux, & + atm_bc_flux, & + seaice_bc_flux, & + nhx_dep, & + noy_dep) + + call callTree_leave("MARBL_forcing_from_data_override") + +end subroutine MARBL_forcing_from_data_override + + !> Save a restart file for the forcing fields subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & filename_suffix) @@ -1739,7 +1835,6 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C "The file with the surface salinity toward which to "//& "restore in the variable given by SSS_RESTORE_VAR.", & fail_if_missing=.true.) - if (CS%archaic_OMIP_file) then CS%SST_restore_var = "TEMP" ; CS%SSS_restore_var = "SALT" else @@ -1952,6 +2047,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C call read_netCDF_data(filename, 'gustiness', CS%gust, G%Domain, & rescale=US%Pa_to_RLZ_T2) ! units in file should be [Pa] endif + call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_marbl_tracers, & + default=.false., do_not_log=.true.) ! All parameter settings are now known. @@ -1978,6 +2075,10 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C call SCM_CVmix_tests_surface_forcing_init(Time, G, param_file, CS%SCM_CVmix_tests_CSp) endif + ! Set up MARBL forcing control structure + call MARBL_forcing_init(G, US, param_file, diag, Time, CS%inputdir, CS%use_marbl_tracers, & + CS%marbl_forcing_CSp) + call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles) ! Set up any restart fields associated with the forcing. @@ -2037,6 +2138,7 @@ subroutine surface_forcing_end(CS, fluxes) if (associated(CS)) deallocate(CS) CS => NULL() + call callTree_leave("MARBL_forcing_from_data_override, MOM_surface_forcing.F90") end subroutine surface_forcing_end end module MOM_surface_forcing diff --git a/config_src/external/MARBL/README.md b/config_src/external/MARBL/README.md new file mode 100644 index 0000000000..f19f76dec8 --- /dev/null +++ b/config_src/external/MARBL/README.md @@ -0,0 +1,6 @@ +MARBL +===== + +These APIs reflect those for the MARBL library available at https://github.com/marbl-ecosys/MARBL + +The modules in this directory do not do any computations. They simply reflect the APIs of the above package. diff --git a/config_src/external/MARBL/marbl_constants_mod.F90 b/config_src/external/MARBL/marbl_constants_mod.F90 new file mode 100644 index 0000000000..7a1d44ba97 --- /dev/null +++ b/config_src/external/MARBL/marbl_constants_mod.F90 @@ -0,0 +1,11 @@ +!> A non-functioning template of the MARBL constants module +module marbl_constants_mod + + implicit none + private + + !> Molecular weight of iron + real, public, parameter :: molw_Fe = 55.845 + +end module marbl_constants_mod + diff --git a/config_src/external/MARBL/marbl_interface.F90 b/config_src/external/MARBL/marbl_interface.F90 new file mode 100644 index 0000000000..c31684597c --- /dev/null +++ b/config_src/external/MARBL/marbl_interface.F90 @@ -0,0 +1,134 @@ +!> A non-functioning template of the MARBL interface +module marbl_interface + + use MOM_error_handler, only : MOM_error, FATAL + use marbl_logging, only : marbl_log_type + use marbl_interface_public_types, only : marbl_forcing_fields_type + use marbl_interface_public_types, only : marbl_tracer_metadata_type + use marbl_interface_public_types, only : marbl_saved_state_type + use marbl_interface_public_types, only : marbl_diagnostics_type + use marbl_interface_public_types, only : marbl_domain_type + use marbl_interface_public_types, only : marbl_output_for_GCM_type + implicit none + private ! Only want marbl_interface_class to be public, not supporting functions + + !> A non-functioning template of the MARBL_interface class + !! + !> All variables are dummy representations of actual members of the real marbl_interface_class + !! that are used in the MARBL tracer routines. + type, public :: marbl_interface_class + type(marbl_log_type) :: StatusLog !< dummy log + type(marbl_forcing_fields_type), allocatable :: surface_flux_forcings(:) !< dummy forcing array + type(marbl_forcing_fields_type), allocatable :: interior_tendency_forcings(:) !< dummy forcing array + type(marbl_tracer_metadata_type), allocatable :: tracer_metadata(:) !< dummy metadata array + type(marbl_domain_type) :: domain !< dummy domain + type(marbl_saved_state_type) :: surface_flux_saved_state !< dummy saved state + type(marbl_saved_state_type) :: interior_tendency_saved_state !< dummy saved state + type(marbl_diagnostics_type) :: surface_flux_diags !< dummy diagnostics + type(marbl_diagnostics_type) :: interior_tendency_diags !< dummy diagnostics + type(marbl_output_for_GCM_type) :: surface_flux_output !< dummy output + type(marbl_output_for_GCM_type) :: interior_tendency_output !< dummy output + real, allocatable :: tracers(:,:) !< dummy tracer array + real, allocatable :: tracers_at_surface(:,:) !< dummy tracer surface array + real, allocatable :: bot_flux_to_tend(:) !< dummy array for bot flux to tendency wgts + real, allocatable :: surface_fluxes(:,:) !< dummy fluxes + real, allocatable :: interior_tendencies(:,:) !< dummy tendencies + contains + procedure, public :: put_setting !< dummy put_setting routine + procedure, public :: get_setting !< dummy get_setting routine + procedure, public :: init !< dummy routine + procedure, public :: surface_flux_compute !< dummy surface flux routine + procedure, public :: interior_tendency_compute !< dummy interior tendency routine + procedure, public :: add_output_for_GCM !< dummy add_output_for_GCM routine + procedure, public :: shutdown !< dummy shutdown routine + end type marbl_interface_class + + !> Error message that appears if the dummy interface is called + character(len=*), parameter :: error_msg = "MOM6 built the MARBL stubs rather than the full library" + +contains + + !> Dummy version of MARBL's put_setting() function + subroutine put_setting(self, str_in) + class(marbl_interface_class), intent(in) :: self + character(len=*), intent(in) :: str_in + + call MOM_error(FATAL, error_msg) + end subroutine put_setting + + !> Dummy version of MARBL's get_setting() function + subroutine get_setting(self, str_in, log_out) + class(marbl_interface_class), intent(in) :: self + character(len=*), intent(in) :: str_in + logical, intent(out) :: log_out + + call MOM_error(FATAL, error_msg) + end subroutine get_setting + + !> Dummy version of MARBL's init() function + subroutine init(self, & + gcm_num_levels, & + gcm_num_PAR_subcols, & + gcm_num_elements_surface_flux, & + gcm_delta_z, & + gcm_zw, & + gcm_zt, & + unit_system_opt, & + lgcm_has_global_ops) + + class(marbl_interface_class), intent(inout) :: self + integer, intent(in) :: gcm_num_levels + integer, intent(in) :: gcm_num_PAR_subcols + integer, intent(in) :: gcm_num_elements_surface_flux + real, intent(in) :: gcm_delta_z(gcm_num_levels) + real, intent(in) :: gcm_zw(gcm_num_levels) + real, intent(in) :: gcm_zt(gcm_num_levels) + character(len=*), intent(in) :: unit_system_opt + logical, intent(in) :: lgcm_has_global_ops + + call MOM_error(FATAL, error_msg) + end subroutine init + + !> Dummy version of MARBL's surface_flux_compute() function + subroutine surface_flux_compute(self) + + class(marbl_interface_class), intent(inout) :: self + + call MOM_error(FATAL, error_msg) + + end subroutine surface_flux_compute + + !> Dummy version of MARBL's interior_tendency_compute() function + subroutine interior_tendency_compute(self) + + class(marbl_interface_class), intent(inout) :: self + + call MOM_error(FATAL, error_msg) + + end subroutine interior_tendency_compute + + !> Dummy version of MARBL's add_output_for_GCM() function + subroutine add_output_for_GCM(self, num_elements, field_name, output_id, field_source, num_levels) + + class (marbl_interface_class), intent(inout) :: self + integer, intent(in) :: num_elements + character(len=*), intent(in) :: field_name + integer, intent(out) :: output_id + character(len=*), intent(out) :: field_source + integer, optional, intent(in) :: num_levels + + output_id = 0 + field_source = "" + + end subroutine add_output_for_GCM + + !> Dummy version of MARBL's shutdown() function + subroutine shutdown(self) + + class(marbl_interface_class), intent(inout) :: self + + call MOM_error(FATAL, error_msg) + + end subroutine shutdown + +end module marbl_interface diff --git a/config_src/external/MARBL/marbl_interface_public_types.F90 b/config_src/external/MARBL/marbl_interface_public_types.F90 new file mode 100644 index 0000000000..5c49ea1985 --- /dev/null +++ b/config_src/external/MARBL/marbl_interface_public_types.F90 @@ -0,0 +1,89 @@ +!> A non-functioning template of the public structures provided through MARBL interface +module marbl_interface_public_types + + use marbl_logging, only : marbl_log_type + + implicit none + private ! Only want a few types to be public + + !> A non-functioning template of MARBL diagnostic type + type :: marbl_single_diagnostic_type + character(len=0) :: long_name !< dummy name + character(len=0) :: short_name !< dummy name + character(len=0) :: units !< dummy units + character(len=0) :: vertical_grid !< dummy grid + logical :: compute_now !< dummy flag + logical :: ltruncated_vertical_extent !< dummy flag + integer :: ref_depth !< dummy depth + real, allocatable, dimension(:) :: field_2d !< dummy field + real, allocatable, dimension(:,:) :: field_3d !< dummy field + end type marbl_single_diagnostic_type + + !> A non-functioning template of MARBL diagnostic type + type, public :: marbl_diagnostics_type + type(marbl_single_diagnostic_type), dimension(:), pointer :: diags => NULL() !< dummy point + end type marbl_diagnostics_type + + !> A non-functioning template of MARBL saved state type + type :: marbl_single_saved_state_type + integer :: rank !< dummy rank + character(len=0) :: short_name !< dummy name + character(len=0) :: units !< dummy units + character(len=0) :: vertical_grid !< dummy grid + real, allocatable :: field_2d(:) !< dummy field + real, allocatable :: field_3d(:,:) !< dummy field + end type marbl_single_saved_state_type + + !> A non-functioning template of MARBL saved state type + type, public :: marbl_saved_state_type + integer :: saved_state_cnt !< dummy counter + type(marbl_single_saved_state_type), dimension(:), pointer :: state => NULL() !< dummy pointer + end type marbl_saved_state_type + + !> A non-functioning template of MARBL forcing metadata type + type :: marbl_forcing_fields_metadata_type + character(len=0) :: varname !< dummy name + end type marbl_forcing_fields_metadata_type + + !> A non-functioning template of MARBL forcing type + type, public :: marbl_forcing_fields_type + type(marbl_forcing_fields_metadata_type) :: metadata !< dummy metadata + real, pointer :: field_0d(:) => NULL() !< dummy pointer + real, pointer :: field_1d(:,:) => NULL() !< dummy pointer + end type marbl_forcing_fields_type + + !> A non-functioning template of MARBL tracer metadata type + type, public :: marbl_tracer_metadata_type + character(len=0) :: short_name !< dummy name + character(len=0) :: long_name !< dummy name + character(len=0) :: units !< dummy units + end type marbl_tracer_metadata_type + + !> A non-functioning template of MARBL domain type + type, public :: marbl_domain_type + integer :: kmt !< dummy index + integer :: km !< dummy index + real, allocatable :: zt(:) !< dummy depths + real, allocatable :: zw(:) !< dummy depths + real, allocatable :: delta_z(:) !< dummy thicknesses + end type marbl_domain_type + + !> A non-functioning template of MARBL single output type + type, public :: marbl_single_output_type + ! marbl_single_output : + ! a private type, this contains both the metadata and + ! the actual data for a single field computed in either + ! surface_flux_compute() or interior_tendency_compute() + ! that needs to be passed to the GCM / flux coupler. + ! Data must be accessed via the marbl_output_for_GCM_type + ! data structure. + real, allocatable, dimension(:) :: forcing_field_0d !< dummy forcing_field_0d + real, allocatable, dimension(:,:) :: forcing_field_1d !< forcing_field_1d + end type marbl_single_output_type + + !> A non-functioning template of MARBL output for GCM type + type, public :: marbl_output_for_GCM_type + type(marbl_single_output_type), dimension(:), pointer :: outputs_for_GCM => NULL() !< dummy outputs_for_GCM + end type marbl_output_for_GCM_type + +end module marbl_interface_public_types \ No newline at end of file diff --git a/config_src/external/MARBL/marbl_logging.F90 b/config_src/external/MARBL/marbl_logging.F90 new file mode 100644 index 0000000000..906d881f0e --- /dev/null +++ b/config_src/external/MARBL/marbl_logging.F90 @@ -0,0 +1,38 @@ +!> A non-functioning template of the MARBL logging module +module marbl_logging + + implicit none + private + + !> A non-functioning template of the marbl status log type + type, public :: marbl_status_log_entry_type + integer :: ElementInd !< dummy index + logical :: lonly_master_writes !< dummy flag + character(len=0) :: LogMessage !< dummy message + type(marbl_status_log_entry_type), pointer :: next !< dummy pointer + end type marbl_status_log_entry_type + + !> A non-functioning template of the marbl status log type + type, public :: marbl_log_type + logical, public :: labort_marbl !< dummy flag + type(marbl_status_log_entry_type), pointer :: FullLog !< dummy pointer + contains + procedure, public :: log_error_trace !< dummy trace routine + procedure, public :: erase !< dummy erase routine + end type marbl_log_type + +contains + + !> dummy trace routine + subroutine log_error_trace(self, RoutineName, CodeLoc, ElemInd) + class(marbl_log_type), intent(inout) :: self + character(len=*), intent(in) :: RoutineName, CodeLoc + integer, optional, intent(in) :: ElemInd + end subroutine log_error_trace + + !> dummy erase routine + subroutine erase(self) + class(marbl_log_type), intent(inout) :: self + end subroutine erase + +end module marbl_logging \ No newline at end of file diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 9423197f89..52aac958e0 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 9423197f894112edfcb1502245f7d7b873d551f9 +Subproject commit 52aac958e05cdb2471dc73f9ef7fb4e816c550f2 diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index b8b3174b4a..6e4969142e 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -53,6 +53,12 @@ module MOM_forcing_type module procedure allocate_mech_forcing_from_ref end interface allocate_mech_forcing +!> Allocate arrays if optional flag is present and true (works for 2D and 3D) +interface myAlloc + module procedure myAlloc_2d + module procedure myAlloc_3d +end interface myAlloc + !> Determine the friction velocity from a forcing type or a mechanical forcing type. interface find_ustar module procedure find_ustar_fluxes @@ -212,6 +218,19 @@ module MOM_forcing_type ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] + ! Forcing fields required for MARBL + real, pointer, dimension(:,:) :: & + noy_dep => NULL(), & !< NOy Deposition [conc Z T-1 ~> conc m s-1] + nhx_dep => NULL(), & !< NHx Deposition [conc Z T-1 ~> conc m s-1] + atm_co2 => NULL(), & !< Atmospheric CO2 Concentration [ppm] + atm_alt_co2 => NULL(), & !< Alternate atmospheric CO2 Concentration [ppm] + dust_flux => NULL(), & !< Flux of dust into the ocean [R Z T-1 ~> kgN m-2 s-1] + iron_flux => NULL() !< Flux of dust into the ocean [conc Z T-1 ~> conc m s-1] + + real, pointer, dimension(:,:,:) :: & + fracr_cat => NULL(), & !< per-category ice fraction + qsw_cat => NULL() !< per-category shortwave + real, pointer, dimension(:,:) :: & lamult => NULL() !< Langmuir enhancement factor [nondim] @@ -3202,8 +3221,9 @@ end subroutine forcing_diagnostics !> Conditionally allocate fields within the forcing type subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & - shelf, iceberg, salt, fix_accum_bug, cfc, waves, & - shelf_sfc_accumulation, lamult, hevap, tau_mag) + shelf, iceberg, salt, fix_accum_bug, cfc, marbl, & + waves, shelf_sfc_accumulation, lamult, hevap, & + ice_ncat, tau_mag) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields logical, optional, intent(in) :: water !< If present and true, allocate water fluxes @@ -3217,6 +3237,8 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & !! accumulation of ustar_gustless logical, optional, intent(in) :: cfc !< If present and true, allocate fields needed !! for cfc surface fluxes + logical, optional, intent(in) :: marbl !< If present and true, allocate fields needed + !! for MARBL surface fluxes logical, optional, intent(in) :: waves !< If present and true, allocate wave fields logical, optional, intent(in) :: shelf_sfc_accumulation !< If present and true, and shelf is true, !! then allocate surface flux deposition from the atmosphere @@ -3225,6 +3247,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: hevap !< If present and true, allocate heat content evap. !! This field must be allocated when enthalpy is provided !! via coupler. + integer, optional, intent(in) :: ice_ncat !< number of ice categories logical, optional, intent(in) :: tau_mag !< If present and true, allocate tau_mag and related fields ! Local variables @@ -3291,20 +3314,37 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & if (shelf_sfc_acc) call myAlloc(fluxes%shelf_sfc_mass_flux,isd,ied,jsd,jed, shelf_sfc_acc) endif; endif - !These fields should only on allocated when iceberg area is being passed through the coupler. + !These fields should only be allocated when iceberg area is being passed through the coupler. call myAlloc(fluxes%ustar_berg,isd,ied,jsd,jed, iceberg) call myAlloc(fluxes%area_berg,isd,ied,jsd,jed, iceberg) call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) - !These fields should only on allocated when USE_CFC_CAP is activated. + !These fields should only be allocated when USE_CFC_CAP is activated. call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, cfc) call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, cfc) - !These fields should only on allocated when wave coupling is activated. + !These fields should only be allocated when wave coupling is activated. call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, waves) call myAlloc(fluxes%lamult,isd,ied,jsd,jed, lamult) if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug + + !These fields should only be allocated when USE_MARBL is activated. + call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%noy_dep,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%nhx_dep,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%atm_co2,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%atm_alt_co2,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%dust_flux,isd,ied,jsd,jed, marbl) + call myAlloc(fluxes%iron_flux,isd,ied,jsd,jed, marbl) + + ! These fields should only be allocated when receiving multiple ice categories + if (present(ice_ncat)) then + call myAlloc(fluxes%fracr_cat,isd,ied,jsd,jed,1,ice_ncat+1, ice_ncat > 0) + call myAlloc(fluxes%qsw_cat,isd,ied,jsd,jed,1,ice_ncat+1, ice_ncat > 0) + endif + end subroutine allocate_forcing_by_group !> Allocate elements of a new forcing type based on their status in an existing type. @@ -3495,7 +3535,7 @@ end subroutine get_mech_forcing_groups !> Allocates and zeroes-out array. -subroutine myAlloc(array, is, ie, js, je, flag) +subroutine myAlloc_2d(array, is, ie, js, je, flag) real, dimension(:,:), pointer :: array !< Array to be allocated integer, intent(in) :: is !< Start i-index integer, intent(in) :: ie !< End i-index @@ -3506,7 +3546,22 @@ subroutine myAlloc(array, is, ie, js, je, flag) if (present(flag)) then ; if (flag) then ; if (.not.associated(array)) then allocate(array(is:ie,js:je), source=0.0) endif ; endif ; endif -end subroutine myAlloc +end subroutine myAlloc_2d + +subroutine myAlloc_3d(array, is, ie, js, je, ks, ke, flag) + real, dimension(:,:,:), pointer :: array !< Array to be allocated + integer, intent(in) :: is !< Start i-index + integer, intent(in) :: ie !< End i-index + integer, intent(in) :: js !< Start j-index + integer, intent(in) :: je !< End j-index + integer, intent(in) :: ks !< Start k-index + integer, intent(in) :: ke !< End k-index + logical, optional, intent(in) :: flag !< Flag to indicate to allocate + + if (present(flag)) then ; if (flag) then ; if (.not.associated(array)) then + allocate(array(is:ie,js:je,ks:ke), source=0.0) + endif ; endif ; endif +end subroutine myAlloc_3d !> Deallocate the forcing type subroutine deallocate_forcing_type(fluxes) @@ -3562,6 +3617,14 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) if (associated(fluxes%ice_fraction)) deallocate(fluxes%ice_fraction) if (associated(fluxes%u10_sqr)) deallocate(fluxes%u10_sqr) + if (associated(fluxes%noy_dep)) deallocate(fluxes%noy_dep) + if (associated(fluxes%nhx_dep)) deallocate(fluxes%nhx_dep) + if (associated(fluxes%atm_co2)) deallocate(fluxes%atm_co2) + if (associated(fluxes%atm_alt_co2)) deallocate(fluxes%atm_alt_co2) + if (associated(fluxes%dust_flux)) deallocate(fluxes%dust_flux) + if (associated(fluxes%iron_flux)) deallocate(fluxes%iron_flux) + if (associated(fluxes%fracr_cat)) deallocate(fluxes%fracr_cat) + if (associated(fluxes%qsw_cat)) deallocate(fluxes%qsw_cat) call coupler_type_destructor(fluxes%tr_fluxes) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 0eab1a5b17..c432e73223 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -57,7 +57,8 @@ module MOM_variables ocean_heat, & !< The total heat content of the ocean in [C R Z ~> degC kg m-2]. ocean_salt, & !< The total salt content of the ocean in [1e-3 S R Z ~> kgSalt m-2]. taux_shelf, & !< The zonal stresses on the ocean under shelves [R L Z T-2 ~> Pa]. - tauy_shelf !< The meridional stresses on the ocean under shelves [R L Z T-2 ~> Pa]. + tauy_shelf, & !< The meridional stresses on the ocean under shelves [R L Z T-2 ~> Pa]. + fco2 !< CO2 flux from the ocean to the atmosphere [R Z T-1 ~> kgCO2 m-2 s-1] logical :: T_is_conT = .false. !< If true, the temperature variable SST is actually the !! conservative temperature in [C ~> degC]. logical :: S_is_absS = .false. !< If true, the salinity variable SSS is actually the @@ -337,7 +338,7 @@ module MOM_variables !! the ocean model. Unused fields are unallocated. subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & gas_fields_ocn, use_meltpot, use_iceshelves, & - omit_frazil) + omit_frazil, use_marbl_tracers) type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(surface), intent(inout) :: sfc_state !< ocean surface state type to be allocated. logical, optional, intent(in) :: use_temperature !< If true, allocate the space for thermodynamic variables. @@ -354,9 +355,10 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & !! under ice shelves. logical, optional, intent(in) :: omit_frazil !< If present and false, do not allocate the space to !! pass frazil fluxes to the coupler + logical, optional, intent(in) :: use_marbl_tracers !< If true, allocate the space for CO2 flux from MARBL ! local variables - logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil + logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil, alloc_fco2 integer :: is, ie, js, je, isd, ied, jsd, jed integer :: isdB, iedB, jsdB, jedB @@ -369,6 +371,7 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot alloc_iceshelves = .false. ; if (present(use_iceshelves)) alloc_iceshelves = use_iceshelves alloc_frazil = .true. ; if (present(omit_frazil)) alloc_frazil = .not.omit_frazil + alloc_fco2 = .false. ; if (present(use_marbl_tracers)) alloc_fco2 = use_marbl_tracers if (sfc_state%arrays_allocated) return @@ -408,6 +411,10 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & call coupler_type_spawn(gas_fields_ocn, sfc_state%tr_fields, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) + if (alloc_fco2) then + allocate(sfc_state%fco2(isd:ied,jsd:jed), source=0.0) + endif + sfc_state%arrays_allocated = .true. end subroutine allocate_surface_state @@ -429,6 +436,7 @@ subroutine deallocate_surface_state(sfc_state) if (allocated(sfc_state%ocean_mass)) deallocate(sfc_state%ocean_mass) if (allocated(sfc_state%ocean_heat)) deallocate(sfc_state%ocean_heat) if (allocated(sfc_state%ocean_salt)) deallocate(sfc_state%ocean_salt) + if (allocated(sfc_state%fco2)) deallocate(sfc_state%fco2) call coupler_type_destructor(sfc_state%tr_fields) sfc_state%arrays_allocated = .false. diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 index e131e8db9d..9c725f7f1a 100644 --- a/src/framework/MOM_interpolate.F90 +++ b/src/framework/MOM_interpolate.F90 @@ -10,13 +10,27 @@ module MOM_interpolate use MOM_interp_infra, only : horiz_interp_type, get_external_field_info use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights use MOM_interp_infra, only : external_field -use MOM_time_manager, only : time_type +use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(<), operator(>) implicit none ; private +!> Data type used to store information about forcing datasets that are time series +!! E.g. how do we align the data in the model with the time axis in the file? +type, public :: forcing_timeseries_dataset + character(len=200) :: file_name !< name of file containing river flux forcing + logical :: l_time_varying !< .true. => forcing is dependent on model time, .false. => static forcing + ! logical :: l_FMS_modulo !< .true. => let FMS handle determining time level to read (e.g. for climatologies) + type(time_type) :: data_forcing !< convert data_forcing_year to time type + type(time_type) :: data_start !< convert data_start_year to time type + type(time_type) :: data_end !< convert data_end_year to time type + type(time_type) :: m2d_offset !< add to model time to get data time +end type forcing_timeseries_dataset + public :: time_interp_external, init_external_field, time_interp_external_init, get_external_field_info public :: horiz_interp_type, run_horiz_interp, build_horiz_interp_weights public :: external_field +public :: forcing_timeseries_set_time_type_vars +public :: map_model_time_to_forcing_time !> Read a field based on model time, and rotate to the model domain. interface time_interp_external @@ -210,4 +224,65 @@ subroutine time_interp_external_3d(field, time, data_in, interp, & end subroutine time_interp_external_3d +!> Set time_type variables in forcing_timeseries_dataset type based on integer input +!! TODO: make this part of forcing_timeseries_dataset class if OO is okay in MOM6? +subroutine forcing_timeseries_set_time_type_vars(data_start_year, data_end_year, data_ref_year, & + model_ref_year, data_forcing_year, forcing_dataset) + + integer, intent(in) :: data_start_year !< first year of data to read + !! (this is ignored for static forcing) + integer, intent(in) :: data_end_year !< last year of data to read + !! (this is ignored for static forcing) + integer, intent(in) :: data_ref_year !< for time-varying forcing, align + !! data_ref_year in file with + !! model_ref_year in model + integer, intent(in) :: model_ref_year !< for time-varying forcing, align + !! data_ref_year in file with + !! model_ref_year in model + integer, intent(in) :: data_forcing_year !< for static forcing, read file at this + !! date (this is ignored for time-varying + !! forcing) + type(forcing_timeseries_dataset), intent(inout) :: forcing_dataset !< information about forcing file + + if (forcing_dataset%l_time_varying) then + forcing_dataset%data_start = real_to_time(year_to_sec(data_start_year)) + forcing_dataset%data_end = real_to_time(year_to_sec(data_end_year)) + forcing_dataset%m2d_offset = real_to_time(year_to_sec(data_ref_year - model_ref_year)) + else + forcing_dataset%data_forcing = real_to_time(year_to_sec(data_forcing_year)) + endif + +end subroutine forcing_timeseries_set_time_type_vars + +!> If necessary, apply an offset to convert from model time to forcing time and then +!! ensure result is within acceptable bounds +function map_model_time_to_forcing_time(Time, forcing_dataset) + + type(time_type), intent(in) :: Time !< Model time + type(forcing_timeseries_dataset), intent(in) :: forcing_dataset !< information about forcing file + type(time_type) :: map_model_time_to_forcing_time !< time to read forcing file + + if (forcing_dataset%l_time_varying) then + map_model_time_to_forcing_time = Time + forcing_dataset%m2d_offset + ! If Time + offset is not between data_start and data_end, use whichever of those values is closer + if (map_model_time_to_forcing_time < forcing_dataset%data_start) & + map_model_time_to_forcing_time = forcing_dataset%data_start + if (map_model_time_to_forcing_time > forcing_dataset%data_end) & + map_model_time_to_forcing_time = forcing_dataset%data_end + else + map_model_time_to_forcing_time = forcing_dataset%data_forcing + endif + +end function map_model_time_to_forcing_time + +!> real_to_time converts from seconds since 0001-01-01 to time_type so we need to convert from years -> seconds +function year_to_sec(year) + + integer, intent(in) :: year + real :: year_to_sec + + year_to_sec = 86400. * 365. * real(year-1) + +end function year_to_sec + end module MOM_interpolate diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 808430df2c..fa39971d70 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -36,7 +36,7 @@ module MOM_tracer_initialization_from_Z !> Initializes a tracer from a z-space data file, including any lateral regridding that is needed. subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_nam, & src_var_unit_conversion, src_var_record, homogenize, & - useALEremapping, remappingScheme, src_var_gridspec ) + useALEremapping, remappingScheme, src_var_gridspec, ongrid ) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -54,6 +54,9 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ character(len=*), optional, intent(in) :: remappingScheme !< remapping scheme to use. character(len=*), optional, intent(in) :: src_var_gridspec !< Source variable name in a gridspec file. !! This is not implemented yet. + logical, optional, intent(in) :: ongrid !< If true, then data are assumed to have been interpolated to + !! the model horizontal grid. In this case, only extrapolation + !! is performed by horiz_interp_and_extrap_tracer() ! Local variables real :: land_fill = 0.0 ! A value to use to replace missing values [CU ~> conc] real :: convert ! A conversion factor into the model's internal units [CU conc-1 ~> 1] @@ -111,10 +114,10 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ "initial conditions.", default=.false.) call get_param(PF, mdl, "Z_INIT_ALE_REMAPPING", useALE, & "If True, then remap straight to model coordinate from file.",& - default=.true.) + default=.false.) call get_param(PF, mdl, "Z_INIT_REMAPPING_SCHEME", remapScheme, & "The remapping scheme to use if using Z_INIT_ALE_REMAPPING is True.", & - default="PLM") + default="PPM_IH4") call get_param(PF, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) @@ -145,7 +148,8 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ call horiz_interp_and_extrap_tracer(src_file, src_var_nam, recnum, & G, tr_z, mask_z, z_in, z_edges_in, missing_value, & - scale=convert, homogenize=homog, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date) + scale=convert, homogenize=homog, m_to_Z=US%m_to_Z, & + answer_date=hor_regrid_answer_date, ongrid=ongrid) kd = size(z_edges_in,1)-1 call pass_var(tr_z,G%Domain) diff --git a/src/parameterizations/MARBL b/src/parameterizations/MARBL new file mode 120000 index 0000000000..c78d57b86a --- /dev/null +++ b/src/parameterizations/MARBL @@ -0,0 +1 @@ +../../pkg/MARBL/src/ \ No newline at end of file diff --git a/src/tracer/MARBL_forcing_mod.F90 b/src/tracer/MARBL_forcing_mod.F90 new file mode 100644 index 0000000000..9375f9ab08 --- /dev/null +++ b/src/tracer/MARBL_forcing_mod.F90 @@ -0,0 +1,378 @@ +!> This module provides a common datatype to provide forcing for MARBL tracers +!! regardless of driver +module MARBL_forcing_mod + +!! This module exists to house code used by multiple drivers in config_src/ +!! for passing forcing fields to MARBL +!! (This comment can go in the wiki on the NCAR fork?) + +use MOM_diag_mediator, only : safe_alloc_ptr, diag_ctrl, register_diag_field, post_data +use MOM_time_manager, only : time_type +use MOM_error_handler, only : MOM_error, WARNING, FATAL +use MOM_file_parser, only : get_param, log_param, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_interpolate, only : external_field, init_external_field, time_interp_external +use MOM_io, only : slasher +use marbl_constants_mod, only : molw_Fe +use MOM_forcing_type, only : forcing + +implicit none ; private + +#include + +public :: MARBL_forcing_init +public :: convert_driver_fields_to_forcings + +!> Data type used to store diagnostic index returned from register_diag_field() +!! For the forcing fields that can be written via post_data() +type, private :: marbl_forcing_diag_ids + integer :: atm_fine_dust !< Atmospheric fine dust component of dust_flux + integer :: atm_coarse_dust !< Atmospheric coarse dust component of dust_flux + integer :: atm_bc !< Atmospheric black carbon component of iron_flux + integer :: ice_dust !< Sea-ice dust component of dust_flux + integer :: ice_bc !< Sea-ice black carbon component of iron_flux +end type marbl_forcing_diag_ids + +!> Control structure for this module +type, public :: marbl_forcing_CS + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. + + real :: dust_ratio_thres !< coarse/fine dust ratio threshold + real :: dust_ratio_to_fe_bioavail_frac !< ratio of dust to iron bioavailability fraction + real :: fe_bioavail_frac_offset !< offset for iron bioavailability fraction + real :: atm_fe_to_bc_ratio !< atmospheric iron to black carbon ratio + real :: atm_bc_fe_bioavail_frac !< atmospheric black carbon to iron bioavailablity fraction ratio + real :: seaice_fe_to_bc_ratio !< sea-ice iron to black carbon ratio + real :: seaice_bc_fe_bioavail_frac !< sea-ice black carbon to iron bioavailablity fraction ratio + real :: iron_frac_in_atm_fine_dust !< Fraction of fine dust from the atmosphere that is iron + real :: iron_frac_in_atm_coarse_dust !< Fraction of coarse dust from the atmosphere that is iron + real :: iron_frac_in_seaice_dust !< Fraction of dust from the sea ice that is iron + real :: atm_co2_const !< atmospheric CO2 (if specifying a constant value) [ppm] + real :: atm_alt_co2_const !< alternate atmospheric CO2 for _ALT_CO2 tracers + !! (if specifying a constant value) [ppm] + + type(marbl_forcing_diag_ids) :: diag_ids !< used for registering and posting some MARBL forcing fields as diagnostics + + logical :: use_marbl_tracers !< most functions can return immediately + !! MARBL tracers are turned off + integer :: atm_co2_iopt !< Integer version of atm_co2_opt, which determines source of atm_co2 + integer :: atm_alt_co2_iopt !< Integer version of atm_alt_co2_opt, which determines source of atm_alt_co2 + +end type marbl_forcing_CS + +! Module parameters +integer, parameter :: atm_co2_constant_iopt = 0 !< module parameter denoting atm_co2_opt = 'constant' +integer, parameter :: atm_co2_prognostic_iopt = 1 !< module parameter denoting atm_co2_opt = 'diagnostic' +integer, parameter :: atm_co2_diagnostic_iopt = 2 !< module parameter denoting atm_co2_opt = 'prognostic' + +contains + + subroutine MARBL_forcing_init(G, US, param_file, diag, day, inputdir, use_marbl, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + character(len=*), intent(in) :: inputdir !< Directory containing input files + logical, intent(in) :: use_marbl !< Is MARBL tracer package active? + type(marbl_forcing_CS), pointer, intent(inout) :: CS !< A pointer that is set to point to control + !! structure for MARBL forcing + + character(len=40) :: mdl = "MARBL_forcing_mod" ! This module's name. + character(len=15) :: atm_co2_opt + character(len=200) :: err_message + + if (associated(CS)) then + call MOM_error(WARNING, "marbl_forcing_init called with an associated control structure.") + return + endif + + allocate(CS) + CS%diag => diag + + CS%use_marbl_tracers = .true. + if (.not. use_marbl) then + CS%use_marbl_tracers = .false. + return + endif + + call get_param(param_file, mdl, "DUST_RATIO_THRES", CS%dust_ratio_thres, & + "TODO: Add description", units="add_units", default=69.00594) + call get_param(param_file, mdl, "DUST_RATIO_TO_FE_BIOAVAIL_FRAC", & + CS%dust_ratio_to_fe_bioavail_frac, & + "TODO: Add description", units="add_units", default=1./366.314) + call get_param(param_file, mdl, "FE_BIOAVAIL_FRAC_OFFSET", CS%fe_bioavail_frac_offset, & + "TODO: Add description", units="add_units", default=0.0146756) + call get_param(param_file, mdl, "ATM_FE_TO_BC_RATIO", CS%atm_fe_to_bc_ratio, & + "TODO: Add description", units="add_units", default=1.) + call get_param(param_file, mdl, "ATM_BC_FE_BIOAVAIL_FRAC", CS%atm_bc_fe_bioavail_frac, & + "TODO: Add description", units="add_units", default=0.06) + call get_param(param_file, mdl, "SEAICE_FE_TO_BC_RATIO", CS%seaice_fe_to_bc_ratio, & + "TODO: Add description", units="add_units", default=1.) + call get_param(param_file, mdl, "SEAICE_BC_FE_BIOAVAIL_FRAC", CS%seaice_bc_fe_bioavail_frac, & + "TODO: Add description", units="add_units", default=0.06) + call get_param(param_file, mdl, "IRON_FRAC_IN_ATM_FINE_DUST", CS%iron_frac_in_atm_fine_dust, & + "Fraction of fine dust from the atmosphere that is iron", units="add_units", default=0.035) + call get_param(param_file, mdl, "IRON_FRAC_IN_ATM_COARSE_DUST", & + CS%iron_frac_in_atm_coarse_dust, & + "Fraction of coarse dust from the atmosphere that is iron", units="add_units", & + default=0.035) + call get_param(param_file, mdl, "IRON_FRAC_IN_SEAICE_DUST", CS%iron_frac_in_seaice_dust, & + "Fraction of dust from sea ice that is iron", units="add_units", default=0.035) + call get_param(param_file, mdl, "ATM_CO2_OPT", atm_co2_opt, & + "Source of atmospheric CO2 [constant, diagnostic, or prognostic]", & + default="constant") + select case (trim(atm_co2_opt)) + case("prognostic") + CS%atm_co2_iopt = atm_co2_prognostic_iopt + case("diagnostic") + CS%atm_co2_iopt = atm_co2_diagnostic_iopt + case("constant") + CS%atm_co2_iopt = atm_co2_constant_iopt + case DEFAULT + write(err_message, "(3A)") "'", trim(atm_co2_opt), "' is not a valid ATM_CO2_OPT value" + call MOM_error(FATAL, err_message) + end select + if (CS%atm_co2_iopt == atm_co2_constant_iopt) then + call get_param(param_file, mdl, "ATM_CO2_CONST", CS%atm_co2_const, & + "Value to send to MARBL as xco2", & + default=284.317, units="ppm") + endif + call get_param(param_file, mdl, "ATM_ALT_CO2_OPT", atm_co2_opt, & + "Source of alternate atmospheric CO2 [constant, diagnostic, or prognostic]", & + default="constant") + select case (trim(atm_co2_opt)) + case("prognostic") + CS%atm_alt_co2_iopt = atm_co2_prognostic_iopt + case("diagnostic") + CS%atm_alt_co2_iopt = atm_co2_diagnostic_iopt + case("constant") + CS%atm_alt_co2_iopt = atm_co2_constant_iopt + case DEFAULT + write(err_message, "(3A)") "'", trim(atm_co2_opt), "' is not a valid ATM_ALT_CO2_OPT value" + call MOM_error(FATAL, err_message) + end select + if (CS%atm_alt_co2_iopt == atm_co2_constant_iopt) then + call get_param(param_file, mdl, "ATM_ALT_CO2_CONST", CS%atm_alt_co2_const, & + "Value to send to MARBL as xco2_alt_co2", & + default=284.317, units="ppm") + endif + + ! Register diagnostic fields for outputing forcing values + ! These fields are posted from convert_driver_fields_to_forcings(), and they are received + ! in physical units so no conversion is necessary here. + CS%diag_ids%atm_fine_dust = register_diag_field("ocean_model", "ATM_FINE_DUST_FLUX_CPL", & + CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "ATM_FINE_DUST_FLUX from cpl", "kg/m^2/s") + CS%diag_ids%atm_coarse_dust = register_diag_field("ocean_model", "ATM_COARSE_DUST_FLUX_CPL", & + CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "ATM_COARSE_DUST_FLUX from cpl", "kg/m^2/s") + CS%diag_ids%atm_bc = register_diag_field("ocean_model", "ATM_BLACK_CARBON_FLUX_CPL", & + CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "ATM_BLACK_CARBON_FLUX from cpl", "kg/m^2/s") + + CS%diag_ids%ice_dust = register_diag_field("ocean_model", "SEAICE_DUST_FLUX_CPL", & + CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "SEAICE_DUST_FLUX from cpl", "kg/m^2/s") + CS%diag_ids%ice_bc = register_diag_field("ocean_model", "SEAICE_BLACK_CARBON_FLUX_CPL", & + CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "SEAICE_BLACK_CARBON_FLUX from cpl", "kg/m^2/s") + + end subroutine MARBL_forcing_init + + ! Note: ice fraction and u10_sqr are handled in mom_surface_forcing because of CFCs + subroutine convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, & + seaice_dust_flux, atm_bc_flux, seaice_bc_flux, & + nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, & + afracr, swnet_afracr, ifrac_n, & + swpen_ifrac_n, Time, G, US, i0, j0, fluxes, CS) + + real, dimension(:,:), pointer, intent(in) :: atm_fine_dust_flux !< atmosphere fine dust flux from IOB + !! [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: atm_coarse_dust_flux !< atmosphere coarse dust flux from IOB + !! [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: seaice_dust_flux !< sea ice dust flux from IOB [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: atm_bc_flux !< atmosphere black carbon flux from IOB + !! [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: seaice_bc_flux !< sea ice black carbon flux from IOB + !! [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: afracr !< open ocean fraction + real, dimension(:,:), pointer, intent(in) :: nhx_dep !< NHx flux from atmosphere [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: noy_dep !< NOy flux from atmosphere [kg m-2 s-1] + real, dimension(:,:), pointer, intent(in) :: atm_co2_prog !< Prognostic atmospheric CO2 concentration + real, dimension(:,:), pointer, intent(in) :: atm_co2_diag !< Diagnostic atmospheric CO2 concentration + real, dimension(:,:), pointer, intent(in) :: swnet_afracr !< shortwave flux * open ocean fraction + real, dimension(:,:,:), pointer, intent(in) :: ifrac_n !< per-category ice fraction + real, dimension(:,:,:), pointer, intent(in) :: swpen_ifrac_n !< per-category shortwave flux * ice fraction + type(time_type), intent(in) :: Time !< The time of the fluxes, used for + !! interpolating the salinity to the + !! right time, when it is being + !! restored. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, intent(in) :: i0 !< i index offset + integer, intent(in) :: j0 !< j index offset + type(forcing), intent(inout) :: fluxes !< MARBL-specific forcing fields + type(marbl_forcing_CS), pointer, intent(inout) :: CS !< A pointer that is set to point to + !! control structure for MARBL forcing + + integer :: i, j, is, ie, js, je, m + real :: atm_fe_bioavail_frac !< TODO: define this (local) term + real :: seaice_fe_bioavail_frac !< TODO: define this (local) term + real :: iron_flux_conversion !< TODO: define this (local) term + real :: ndep_conversion !< Combination of unit conversion factors for rescaling + !! nitrogen deposition [kg(N) m-2 s-1 ~> mol m-3 Z T-1] + + if (.not. CS%use_marbl_tracers) return + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + ndep_conversion = (1.e6/14.) * (US%m_to_Z * US%T_to_s) ! kg / m^2 / s -> conc Z / T + iron_flux_conversion = (1.e6 / molw_Fe) * (US%m_to_Z * US%T_to_s) ! kg / m^2 / s -> conc Z / T + + ! Post fields from coupler to diagnostics + ! TODO: units from diag register are incorrect; we should be converting these in the cap, I think + if (CS%diag_ids%atm_fine_dust > 0) & + call post_data(CS%diag_ids%atm_fine_dust, atm_fine_dust_flux(is-i0:ie-i0,js-j0:je-j0), & + CS%diag, mask=G%mask2dT(is:ie,js:je)) + if (CS%diag_ids%atm_coarse_dust > 0) & + call post_data(CS%diag_ids%atm_coarse_dust, atm_coarse_dust_flux(is-i0:ie-i0,js-j0:je-j0), & + CS%diag, mask=G%mask2dT(is:ie,js:je)) + if (CS%diag_ids%atm_bc > 0) & + call post_data(CS%diag_ids%atm_bc, atm_bc_flux(is-i0:ie-i0,js-j0:je-j0), CS%diag, & + mask=G%mask2dT(is:ie,js:je)) + if (CS%diag_ids%ice_dust > 0) & + call post_data(CS%diag_ids%ice_dust, seaice_dust_flux(is-i0:ie-i0,js-j0:je-j0), CS%diag, & + mask=G%mask2dT(is:ie,js:je)) + if (CS%diag_ids%ice_bc > 0) & + call post_data(CS%diag_ids%ice_bc, seaice_bc_flux(is-i0:ie-i0,js-j0:je-j0), CS%diag, & + mask=G%mask2dT(is:ie,js:je)) + + do j=js,je ; do i=is,ie + ! Nitrogen Deposition + fluxes%nhx_dep(i,j) = (G%mask2dT(i,j) * ndep_conversion) * nhx_dep(i-i0,j-j0) + fluxes%noy_dep(i,j) = (G%mask2dT(i,j) * ndep_conversion) * noy_dep(i-i0,j-j0) + enddo ; enddo + + ! Atmospheric CO2 + select case (CS%atm_co2_iopt) + case (atm_co2_prognostic_iopt) + if (associated(atm_co2_prog)) then + do j=js,je ; do i=is,ie + fluxes%atm_co2(i,j) = G%mask2dT(i,j) * atm_co2_prog(i-i0,j-j0) + enddo ; enddo + else + call MOM_error(FATAL, & + "ATM_CO2_OPT = 'prognostic' but atmosphere is not providing this field") + endif + case (atm_co2_diagnostic_iopt) + if (associated(atm_co2_diag)) then + do j=js,je ; do i=is,ie + fluxes%atm_co2(i,j) = G%mask2dT(i,j) * atm_co2_diag(i-i0,j-j0) + enddo ; enddo + else + call MOM_error(FATAL, & + "ATM_CO2_OPT = 'diagnostic' but atmosphere is not providing this field") + endif + case (atm_co2_constant_iopt) + do j=js,je ; do i=is,ie + fluxes%atm_co2(i,j) = G%mask2dT(i,j) * CS%atm_co2_const + enddo ; enddo + end select + + ! Alternate Atmospheric CO2 + select case (CS%atm_alt_co2_iopt) + case (atm_co2_prognostic_iopt) + if (associated(atm_co2_prog)) then + do j=js,je ; do i=is,ie + fluxes%atm_alt_co2(i,j) = G%mask2dT(i,j) * atm_co2_prog(i-i0,j-j0) + enddo ; enddo + else + call MOM_error(FATAL, & + "ATM_ALT_CO2_OPT = 'prognostic' but atmosphere is not providing this field") + endif + case (atm_co2_diagnostic_iopt) + if (associated(atm_co2_diag)) then + do j=js,je ; do i=is,ie + fluxes%atm_alt_co2(i,j) = G%mask2dT(i,j) * atm_co2_diag(i-i0,j-j0) + enddo ; enddo + else + call MOM_error(FATAL, & + "ATM_ALT_CO2_OPT = 'diagnostic' but atmosphere is not providing this field") + endif + case (atm_co2_constant_iopt) + do j=js,je ; do i=is,ie + fluxes%atm_alt_co2(i,j) = G%mask2dT(i,j) * CS%atm_co2_const + enddo ; enddo + end select + + ! Dust flux + if (associated(atm_fine_dust_flux)) then + do j=js,je ; do i=is,ie + fluxes%dust_flux(i,j) = US%kg_m2s_to_RZ_T * G%mask2dT(i,j) * & + (atm_fine_dust_flux(i-i0,j-j0) + atm_coarse_dust_flux(i-i0,j-j0) + & + seaice_dust_flux(i-i0,j-j0)) + enddo ; enddo + endif + + if (associated(atm_bc_flux)) then + do j=js,je ; do i=is,ie + ! TODO: abort if atm_fine_dust_flux and atm_coarse_dust_flux are not associated? + ! Contribution of atmospheric dust to iron flux + if (atm_coarse_dust_flux(i-i0,j-j0) < & + CS%dust_ratio_thres * atm_fine_dust_flux(i-i0,j-j0)) then + atm_fe_bioavail_frac = CS%fe_bioavail_frac_offset + CS%dust_ratio_to_fe_bioavail_frac * & + (CS%dust_ratio_thres - atm_coarse_dust_flux(i-i0,j-j0) / atm_fine_dust_flux(i-i0,j-j0)) + else + atm_fe_bioavail_frac = CS%fe_bioavail_frac_offset + endif + + ! Contribution of atmospheric dust to iron flux + fluxes%iron_flux(i,j) = (atm_fe_bioavail_frac * & + (CS%iron_frac_in_atm_fine_dust * atm_fine_dust_flux(i-i0,j-j0) + & + CS%iron_frac_in_atm_coarse_dust * atm_coarse_dust_flux(i-i0,j-j0))) + + ! Contribution of atmospheric black carbon to iron flux + fluxes%iron_flux(i,j) = fluxes%iron_flux(i,j) + (CS%atm_bc_fe_bioavail_frac * & + (CS%atm_fe_to_bc_ratio * atm_bc_flux(i-i0,j-j0))) + + seaice_fe_bioavail_frac = atm_fe_bioavail_frac + ! Contribution of seaice dust to iron flux + fluxes%iron_flux(i,j) = fluxes%iron_flux(i,j) + (seaice_fe_bioavail_frac * & + (CS%iron_frac_in_seaice_dust * seaice_dust_flux(i-i0,j-j0))) + + ! Contribution of seaice black carbon to iron flux + fluxes%iron_flux(i,j) = fluxes%iron_flux(i,j) + (CS%seaice_bc_fe_bioavail_frac * & + (CS%seaice_fe_to_bc_ratio * seaice_bc_flux(i-i0,j-j0))) + + ! Unit conversion (kg / m^2 / s -> conc Z/T) + fluxes%iron_flux(i,j) = (G%mask2dT(i,j) * iron_flux_conversion) * fluxes%iron_flux(i,j) + + enddo ; enddo + endif + + ! Per ice-category forcings + ! If the cap receives per-category fields, memory should be allocated in fluxes + if (associated(ifrac_n)) then + do j=js,je ; do i=is,ie + fluxes%fracr_cat(i,j,1) = min(1., afracr(i-i0,j-j0)) + fluxes%qsw_cat(i,j,1) = swnet_afracr(i-i0,j-j0) + do m=1,size(ifrac_n, 3) + fluxes%fracr_cat(i,j,m+1) = min(1., ifrac_n(i-i0,j-j0,m)) + fluxes%qsw_cat(i,j,m+1) = swpen_ifrac_n(i-i0,j-j0,m) + enddo + where (fluxes%fracr_cat(i,j,:) > 0.) + fluxes%qsw_cat(i,j,:) = fluxes%qsw_cat(i,j,:) / fluxes%fracr_cat(i,j,:) + elsewhere + fluxes%fracr_cat(i,j,:) = 0. + fluxes%qsw_cat(i,j,:) = 0. + endwhere + fluxes%fracr_cat(i,j,:) = G%mask2dT(i,j) * fluxes%fracr_cat(i,j,:) + fluxes%qsw_cat(i,j,:) = G%mask2dT(i,j) * fluxes%qsw_cat(i,j,:) + enddo; enddo + endif + + end subroutine convert_driver_fields_to_forcings + +end module MARBL_forcing_mod diff --git a/src/tracer/MARBL_tracers.F90 b/src/tracer/MARBL_tracers.F90 new file mode 100644 index 0000000000..9c856fef85 --- /dev/null +++ b/src/tracer/MARBL_tracers.F90 @@ -0,0 +1,2206 @@ +!> A tracer package for tracers computed in the MARBL library +!! +!! Currently configured for use with marbl0.36.0 +!! https://github.com/marbl-ecosys/MARBL/releases/tag/marbl0.36.0 +!! (clone entire repo into pkg/MARBL) +module MARBL_tracers + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_coms, only : EFP_type, root_PE, broadcast +use MOM_debugging, only : hchksum +use MOM_diag_mediator, only : diag_ctrl +use MOM_error_handler, only : is_root_PE, MOM_error, FATAL, WARNING, NOTE +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : external_field, init_external_field, time_interp_external +use MOM_CVMix_KPP, only : KPP_NonLocalTransport, KPP_CS +use MOM_hor_index, only : hor_index_type +use MOM_interpolate, only : forcing_timeseries_dataset +use MOM_interpolate, only : forcing_timeseries_set_time_type_vars +use MOM_interpolate, only : map_model_time_to_forcing_time +use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type +use MOM_remapping, only : reintegrate_column +use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h +use MOM_restart, only : query_initialized, MOM_restart_CS, register_restart_field +use MOM_spatial_means, only : global_mass_int_EFP +use MOM_sponge, only : set_up_sponge_field, sponge_CS +use MOM_time_manager, only : time_type +use MOM_tracer_registry, only : register_tracer +use MOM_tracer_types, only : tracer_type, tracer_registry_type +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z +use MOM_tracer_Z_init, only : read_Z_edges +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface, thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type +use MOM_diag_mediator, only : register_diag_field, post_data!, safe_alloc_ptr + +use MARBL_interface, only : MARBL_interface_class +use MARBL_interface_public_types, only : marbl_diagnostics_type, marbl_saved_state_type + +use coupler_types_mod, only : coupler_type_set_data, ind_csurf +use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux + +implicit none ; private + +#include + +public register_MARBL_tracers, initialize_MARBL_tracers +public MARBL_tracers_column_physics, MARBL_tracers_surface_state +public MARBL_tracers_set_forcing +public MARBL_tracers_stock, MARBL_tracers_get, MARBL_tracers_end + +! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional +! consistency testing. These are noted in comments with units like Z, H, L, and T, along with +! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units +! vary with the Boussinesq approximation, the Boussinesq variant is given first. + +!> Temporary type for diagnostic variables coming from MARBL +!! Allocate exactly one of field_[23]d +type :: temp_MARBL_diag + integer :: id !< index into MOM diagnostic structure + real, allocatable :: field_2d(:,:) !< memory for 2D field + real, allocatable :: field_3d(:,:,:) !< memory for 3D field +end type temp_MARBL_diag + +!> MOM6 needs to know the index of some MARBL tracers to properly apply river fluxes +type :: tracer_ind_type + integer :: no3_ind !< NO3 index + integer :: po4_ind !< PO4 index + integer :: don_ind !< DON index + integer :: donr_ind !< DONr index + integer :: dop_ind !< DOP index + integer :: dopr_ind !< DOPr index + integer :: sio3_ind !< SiO3 index + integer :: fe_ind !< Fe index + integer :: doc_ind !< DOC index + integer :: docr_ind !< DOCr index + integer :: alk_ind !< ALK index + integer :: alk_alt_co2_ind !< ALK_ALT_CO2 index + integer :: dic_ind !< DIC index + integer :: dic_alt_co2_ind !< DIC_ALT_CO2 index +end type tracer_ind_type + +!> MOM needs to store some information about saved_state; besides providing these +!! fields to MARBL, they are also written to restart files +type :: saved_state_for_MARBL_type + character(len=200) :: short_name !< name of variable being saved + character(len=200) :: file_varname !< name of variable in restart file + character(len=200) :: units !< variable units + real, pointer :: field_2d(:,:) !< memory for 2D field + real, pointer :: field_3d(:,:,:) !< memory for 3D field +end type saved_state_for_MARBL_type + +!> All calls to MARBL are done via the interface class +type(MARBL_interface_class) :: MARBL_instances + +!> Pointer to tracer concentration and to tracer_type in tracer registry +type, private :: MARBL_tracer_data + real, pointer :: tr(:,:,:) => NULL() !< The array of tracers used in this subroutine, in g m-3? + type(tracer_type), pointer :: tr_ptr => NULL() !< pointer to tracer inside Tr_reg +end type MARBL_tracer_data + +!> The control structure for the MARBL tracer package +type, public :: MARBL_tracers_CS ; private + integer :: ntr !< The number of tracers that are actually used. + logical :: debug !< If true, write verbose checksums for debugging purposes. + logical :: base_bio_on !< Will MARBL use base biotic tracers? + logical :: abio_dic_on !< Will MARBL use abiotic DIC / DI14C tracers? + logical :: ciso_on !< Will MARBL use isotopic tracers? + + integer :: restore_count !< The number of tracers MARBL is configured to restore + logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler. + logical :: use_ice_category_fields !< Forcing will include multiple ice categories for ice_frac and shortwave + logical :: request_Chl_from_MARBL !< MARBL can provide Chl to use in set_pen_shortwave() + integer :: ice_ncat !< Number of ice categories when use_ice_category_fields = True + real :: IC_min !< Minimum value for tracer initial conditions + character(len=200) :: IC_file !< The file in which the age-tracer initial values cam be found. + logical :: ongrid !< True if IC_file is already interpolated to MOM grid + type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry + type(MARBL_tracer_data), dimension(:), allocatable :: tracer_data !< type containing tracer data and pointer + !! into tracer registry + + integer, allocatable, dimension(:) :: ind_tr !< Indices returned by aof_set_coupler_flux if it is used and the + !! surface tracer concentrations are to be provided to the coupler. + + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control structure + + type(vardesc), allocatable :: tr_desc(:) !< Descriptions and metadata for the tracers + logical :: tracers_may_reinit !< If true the tracers may be initialized if not found in a restart file + + character(len=200) :: fesedflux_file !< name of [netCDF] file containing iron sediment flux + character(len=200) :: feventflux_file !< name of [netCDF] file containing iron vent flux + type(forcing_timeseries_dataset) :: d14c_dataset(3) !< File and time axis information for d14c forcing + real, dimension(3) :: d14c_bands !< forcing is organized into bands: [30 N, 90 N]; [30 S, 30 N]; [90 S, 30 S] + integer :: d14c_id !< id for diagnostic field with d14c forcing + logical :: read_riv_fluxes !< If true, use river fluxes supplied from an input file. + !! This is temporary, we will always read river fluxes + type(forcing_timeseries_dataset) :: riv_flux_dataset !< File and time axis information for river fluxes + character(len=4) :: restoring_source !< location of tracer restoring data + !! valid values: file, none + integer :: restoring_nz !< number of levels in tracer restoring file + real, allocatable, dimension(:) :: & + restoring_z_edges !< The depths of the cell interfaces in the tracer restoring file [Z ~> m] + real, allocatable, dimension(:) :: & + restoring_dz !< The thickness of the cell layers in the tracer restoring file [H ~> m] + integer :: restoring_timescale_nz !< number of levels in tracer restoring timescale file + real, allocatable, dimension(:) :: & + restoring_timescale_z_edges !< The depths of the cell interfaces in the tracer restoring timescale file [Z ~> m] + real, allocatable, dimension(:) :: & + restoring_timescale_dz !< The thickness of the cell layers in the tracer restoring timescale file [H ~> m] + character(len=14) :: restoring_I_tau_source !< location of inverse restoring timescale data + !! valid values: file, grid_dependent + character(len=200) :: restoring_file !< name of [netCDF] file containing tracer restoring data + type(remapping_CS) :: restoring_remapCS !< Remapping parameters and work arrays for tracer restoring / timescale + character(len=200) :: restoring_I_tau_file !< name of [netCDF] file containing inverse restoring timescale + character(len=200) :: restoring_I_tau_var_name !< name of field containing inverse restoring timescale + character(len=35) :: marbl_settings_file !< name of [text] file containing MARBL settings + + real :: bot_flux_mix_thickness !< for bottom flux -> tendency conversion, assume uniform mixing over + !! bottom layer of prescribed thickness [Z ~> m] + real :: Ibfmt !< Reciprocal of bot_flux_mix_thickness [Z-1 ~> m-1] + + type(temp_MARBL_diag), allocatable :: surface_flux_diags(:) !< collect surface flux diagnostics from all columns + !! before posting + type(temp_MARBL_diag), allocatable :: interior_tendency_diags(:) !< collect tendency diagnostics from all columns + !! before posting + type(saved_state_for_MARBL_type), allocatable :: surface_flux_saved_state(:) !< surface_flux saved state + type(saved_state_for_MARBL_type), allocatable :: interior_tendency_saved_state(:) !< interior_tendency saved state + + ! TODO: If we can post data column by column, all we need are integer arrays for ids + ! integer, allocatable :: id_surface_flux_diags(:) !< array of indices for surface_flux diagnostics + ! integer, allocatable :: id_interior_tendency_diags(:) !< array of indices for interior_tendency diagnostics + + type(tracer_ind_type) :: tracer_inds !< Indices to tracers that will have river fluxes added to STF + + !> Need to store global output from both marbl_instance%surface_flux_compute() and + !! marbl_instance%interior_tendency_compute(). For the former, just need id to register + !! because we already copy data into CS%STF; latter requires copying data and indices + !! so currently using temp_MARBL_diag for that. + integer, allocatable :: id_surface_flux_out(:) !< register_diag indices for surface_flux output + type(temp_MARBL_diag), allocatable :: interior_tendency_out(:) !< collect interior tendencies for diagnostic output + type(temp_MARBL_diag), allocatable :: interior_tendency_out_zint(:) !< vertical integral of interior tendencies + !! (full column) + type(temp_MARBL_diag), allocatable :: interior_tendency_out_zint_100m(:) !< vertical integral of interior tendencies + !! (top 100m) + integer :: bot_flux_to_tend_id !< register_diag index for BOT_FLUX_TO_TEND + integer, allocatable :: fracr_cat_id(:) !< register_diag index for per-category ice fraction + integer, allocatable :: qsw_cat_id(:) !< register_diag index for per-category shortwave + + real, allocatable :: STF(:,:,:) !< surface fluxes returned from MARBL to use in tracer_vertdiff + !! (dims: i, j, tracer) [conc Z T-1 ~> conc m s-1] + real, allocatable :: SFO(:,:,:) !< surface flux output returned from MARBL for use in GCM + !! e.g. CO2 flux to pass to atmosphere (dims: i, j, num_sfo) + real, allocatable :: ITO(:,:,:,:) !< interior tendency output returned from MARBL for use in GCM + !! e.g. total chlorophyll to use in shortwave penetration (dims: i, j, k, num_ito) + + integer :: u10_sqr_ind !< index of MARBL forcing field array to copy 10-m wind (squared) into + integer :: sss_ind !< index of MARBL forcing field array to copy sea surface salinity into + integer :: sst_ind !< index of MARBL forcing field array to copy sea surface temperature into + integer :: ifrac_ind !< index of MARBL forcing field array to copy ice fraction into + integer :: dust_dep_ind !< index of MARBL forcing field array to copy dust flux into + integer :: fe_dep_ind !< index of MARBL forcing field array to copy iron flux into + integer :: nox_flux_ind !< index of MARBL forcing field array to copy NOx flux into + integer :: nhy_flux_ind !< index of MARBL forcing field array to copy NHy flux into + integer :: atmpress_ind !< index of MARBL forcing field array to copy atmospheric pressure into + integer :: xco2_ind !< index of MARBL forcing field array to copy CO2 flux into + integer :: xco2_alt_ind !< index of MARBL forcing field array to copy CO2 flux (alternate CO2) into + integer :: d14c_ind !< index of MARBL forcing field array to copy d14C into + + !> external_field types for river fluxes (added to surface fluxes) + type(external_field) :: id_din_riv !< id for time_interp_external. + type(external_field) :: id_don_riv !< id for time_interp_external. + type(external_field) :: id_dip_riv !< id for time_interp_external. + type(external_field) :: id_dop_riv !< id for time_interp_external. + type(external_field) :: id_dsi_riv !< id for time_interp_external. + type(external_field) :: id_dfe_riv !< id for time_interp_external. + type(external_field) :: id_dic_riv !< id for time_interp_external. + type(external_field) :: id_alk_riv !< id for time_interp_external. + type(external_field) :: id_doc_riv !< id for time_interp_external. + + !> external_field type for d14c (needed if abio_dic_on is True) + type(external_field) :: id_d14c(3) !< id for time_interp_external. + + !> Indices for river fluxes (diagnostics) + integer :: no3_riv_flux !< NO3 riverine flux + integer :: po4_riv_flux !< PO4 riverine flux + integer :: don_riv_flux !< DON riverine flux + integer :: donr_riv_flux !< DONr riverine flux + integer :: dop_riv_flux !< DOP riverine flux + integer :: dopr_riv_flux !< DOPr riverine flux + integer :: sio3_riv_flux !< SiO3 riverine flux + integer :: fe_riv_flux !< Fe riverine flux + integer :: doc_riv_flux !< DOC riverine flux + integer :: docr_riv_flux !< DOCr riverine flux + integer :: alk_riv_flux !< ALK riverine flux + integer :: alk_alt_co2_riv_flux !< ALK (alternate CO2) riverine flux + integer :: dic_riv_flux !< DIC riverine flux + integer :: dic_alt_co2_riv_flux !< DIC (alternate CO2) riverine flux + + !> Indices for forcing fields required to compute interior tendencies + integer :: dustflux_ind !< index of MARBL forcing field array to copy dust flux into + integer :: PAR_col_frac_ind !< index of MARBL forcing field array to copy PAR column fraction into + integer :: surf_shortwave_ind !< index of MARBL forcing field array to copy surface shortwave into + integer :: potemp_ind !< index of MARBL forcing field array to copy potential temperature into + integer :: salinity_ind !< index of MARBL forcing field array to copy salinity into + integer :: pressure_ind !< index of MARBL forcing field array to copy pressure into + integer :: fesedflux_ind !< index of MARBL forcing field array to copy iron sediment flux into + integer :: o2_scalef_ind !< index of MARBL forcing field array to copy O2 scale length into + integer :: remin_scalef_ind !< index of MARBL forcing field array to copy remin scale length into + type(external_field), allocatable :: id_tracer_restoring(:) !< id number for time_interp_external + integer, allocatable :: tracer_restoring_ind(:) !< index of MARBL forcing field to copy + !! per-tracer restoring field into + integer, allocatable :: tracer_I_tau_ind(:) !< index of MARBL forcing field to copy per-tracer + !! inverse restoring timescale into + + !> Memory for storing river fluxes, tracer restoring fields, and abiotic forcing + real, allocatable :: d14c(:,:) !< d14c forcing for abiotic DIC and carbon isotope tracer modules + real, allocatable :: RIV_FLUXES(:,:,:) !< river flux forcing for applyTracerBoundaryFluxesInOut + !! (needs to be time-integrated when passed to function!) + !! (dims: i, j, tracer) [conc m s-1] + character(len=15), allocatable :: tracer_restoring_varname(:) !< name of variable being restored + real, allocatable :: I_tau(:,:,:) !< inverse restoring timescale for marbl tracers (dims: i, j, k) [1/s] + real, allocatable, dimension(:,:,:,:) :: restoring_in !< Restoring fields read from file + !! (dims: i, j, restoring_nz, restoring_cnt) [tracer units] + + !> Number of surface flux outputs as well as specific indices for each one + integer :: sfo_cnt !< number of surface flux outputs from MARBL + integer :: ito_cnt !< number of interior tendency outputs from MARBL + integer :: flux_co2_ind !< index to co2 flux surface flux output + integer :: total_Chl_ind !< index to total chlorophyll interior tendency output + + ! TODO: create generic 3D forcing input type to read z coordinate + values + real :: fesedflux_scale_factor !< scale factor for iron sediment flux + integer :: fesedflux_nz !< number of levels in iron sediment flux file + real, allocatable, dimension(:,:,:) :: fesedflux_in !< Field to read iron sediment flux into [conc m s-1] + real, allocatable, dimension(:,:,:) :: feventflux_in !< Field to read iron vent flux into [conc m s-1] + real, allocatable, dimension(:) :: & + fesedflux_z_edges !< The depths of the cell interfaces in the input data [Z ~> m] + ! TODO: this thickness does not need to be 3D, but that's a problem for future Mike + real, allocatable, dimension(:,:,:) :: & + fesedflux_dz !< The thickness of the cell layers in the input data [H ~> m] +end type MARBL_tracers_CS + +! Module parameters +real, parameter :: atm_per_Pa = 1./101325. !< convert from Pa -> atm + +contains + +!> This subroutine is used to read marbl_in, configure MARBL accordingly, and then +!! call MARBL's initialization routine +subroutine configure_MARBL_tracers(GV, US, param_file, CS) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(MARBL_tracers_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module + +# include "version_variable.h" + character(len=40) :: mdl = "MARBL_tracers" ! This module's name. + character(len=256) :: log_message + character(len=256) :: marbl_in_line(1) + character(len=256) :: forcing_sname, field_source + integer :: m, n, nz, marbl_settings_in, read_error, I_tau_count, fi + logical :: chl_from_file, forcing_processed + nz = GV%ke + marbl_settings_in = 615 + + ! (1) Read parameters necessary for general setup of MARBL + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "DEBUG", CS%debug, "If true, write out verbose debugging data.", & + default=.false., debuggingParam=.true.) + call get_param(param_file, mdl, "MARBL_IC_MIN_VAL", CS%IC_min, & + "Minimum value of tracer initial conditions (set to 1e-100 for dim scaling tests)", & + default=0., units="tracer units") + call get_param(param_file, mdl, "MARBL_SETTINGS_FILE", CS%marbl_settings_file, & + "The name of a file from which to read the run-time settings for MARBL.", default="marbl_in") + call get_param(param_file, mdl, "BOT_FLUX_MIX_THICKNESS", CS%bot_flux_mix_thickness, & + "Bottom fluxes are uniformly mixed over layer of this thickness", default=1., units="m", & + scale=US%m_to_Z) + call get_param(param_file, mdl, "USE_ICE_CATEGORIES", CS%use_ice_category_fields, & + "If true, allocate memory for shortwave and ice fraction split by ice thickness category.", & + default=.false.) + call get_param(param_file, mdl, "ICE_NCAT", CS%ice_ncat, & + "Number of ice thickness categories in shortwave and ice fraction forcings.", default=0) + CS%Ibfmt = 1. / CS%bot_flux_mix_thickness + + if (CS%use_ice_category_fields .and. (CS%ice_ncat == 0)) & + call MOM_error(FATAL, & + "Can not configure MARBL to use multiple ice categories without ice_ncat present") + + ! (2) Read marbl settings file and call put_setting() + ! (2a) only master task opens file + if (is_root_PE()) then + ! read the marbl_in into buffer + open(unit=marbl_settings_in, file=CS%marbl_settings_file, iostat=read_error) + if (read_error .ne. 0) then + write(log_message, '(A, I0, 2A)') "IO ERROR ", read_error, " opening namelist file : ", & + trim(CS%marbl_settings_file) + call MOM_error(FATAL, log_message) + endif + endif + + ! (2b) master task reads file and broadcasts line-by-line + marbl_in_line = '' + do + ! i. Read next line on master, iostat value out + ! (Exit loop if read is not successful; either read error or end of file) + if (is_root_PE()) read(marbl_settings_in, "(A)", iostat=read_error) marbl_in_line(1) + call broadcast(read_error, root_PE()) + if (read_error .ne. 0) exit + + ! ii. Broadcast line just read in on root PE to all tasks + call broadcast(marbl_in_line, 256, root_PE()) + + ! iii. All tasks call put_setting (TODO: openMP blocks?) + call MARBL_instances%put_setting(marbl_in_line(1)) + enddo + + ! (2c) we should always reach the EOF to capture the entire file... + if (.not. is_iostat_end(read_error)) then + write(log_message, '(3A, I0)') "IO ERROR reading ", trim(CS%marbl_settings_file), ": ", & + read_error + call MOM_error(FATAL, log_message) + else + if (is_root_PE()) then + write(log_message, '(3A)') "Read '", trim(CS%marbl_settings_file), "' until EOF." + call MOM_error(NOTE, log_message) + endif + endif + if (is_root_PE()) close(marbl_settings_in) + + ! (3) Initialize MARBL and configure MOM6 accordingly + + ! (3a) call marbl%init() + ! TODO: We want to strip gcm_delta_z, gcm_zw, and gcm_zt values out of + ! init because MOM updates them every time step / every column + call MARBL_instances%init(gcm_num_levels = nz, gcm_num_PAR_subcols = CS%ice_ncat + 1, & + gcm_num_elements_surface_flux = 1, & ! FIXME: change to number of grid cells on MPI task + gcm_delta_z = GV%sInterface(2:nz+1) - GV%sInterface(1:nz), gcm_zw = GV%sInterface(2:nz+1), & + gcm_zt = GV%sLayer, unit_system_opt = "mks", lgcm_has_global_ops = .false.) ! FIXME: add global ops + ! Regardless of vertical grid, MOM6 will always use GV%ke levels in all columns + MARBL_instances%domain%kmt = GV%ke + if (MARBL_instances%StatusLog%labort_marbl) & + call MARBL_instances%StatusLog%log_error_trace("MARBL_instances%init", & + "configure_MARBL_tracers") + call print_marbl_log(MARBL_instances%StatusLog) + call MARBL_instances%StatusLog%erase() + CS%ntr = size(MARBL_instances%tracer_metadata) + call marbl_instances%get_setting('base_bio_on', CS%base_bio_on) + call marbl_instances%get_setting('abio_dic_on', CS%abio_dic_on) + call marbl_instances%get_setting('ciso_on', CS%ciso_on) + + ! (3b) Read parameters that depend on how MARBL is configured + if (CS%base_bio_on) then + call get_param(param_file, mdl, "CHL_FROM_FILE", chl_from_file, & + "If true, chl_a is read from a file.", default=.true.) + CS%request_Chl_from_MARBL = (.not. chl_from_file) + else + CS%request_Chl_from_MARBL = .false. + endif + + ! (4) Request fields needed by MOM6 + CS%sfo_cnt = 0 + CS%ito_cnt = 0 + + if (CS%base_bio_on) then + ! CO2 Flux to the atmosphere + call MARBL_instances%add_output_for_GCM(num_elements=1, field_name="flux_co2", & + output_id=CS%flux_co2_ind, field_source=field_source) + if (trim(field_source) == "surface_flux") then + CS%sfo_cnt = CS%sfo_cnt + 1 + else if (trim(field_source) == "interior_tendency") then + CS%ito_cnt = CS%ito_cnt + 1 + end if + + ! Total 3D Chlorophyll + call MARBL_instances%add_output_for_GCM(num_elements=1, num_levels=nz, field_name="total_Chl", & + output_id=CS%total_Chl_ind, field_source=field_source) + if (trim(field_source) == "surface_flux") then + CS%sfo_cnt = CS%sfo_cnt + 1 + else if (trim(field_source) == "interior_tendency") then + CS%ito_cnt = CS%ito_cnt + 1 + end if + end if + + ! (5) Initialize forcing fields + ! i. store all surface forcing indices + CS%u10_sqr_ind = -1 + CS%sss_ind = -1 + CS%sst_ind = -1 + CS%ifrac_ind = -1 + CS%dust_dep_ind = -1 + CS%fe_dep_ind = -1 + CS%nox_flux_ind = -1 + CS%nhy_flux_ind = -1 + CS%atmpress_ind = -1 + CS%xco2_ind = -1 + CS%xco2_alt_ind = -1 + do m=1,size(MARBL_instances%surface_flux_forcings) + select case (trim(MARBL_instances%surface_flux_forcings(m)%metadata%varname)) + case('u10_sqr') + CS%u10_sqr_ind = m + case('sss') + CS%sss_ind = m + case('sst') + CS%sst_ind = m + case('Ice Fraction') + CS%ifrac_ind = m + case('Dust Flux') + CS%dust_dep_ind = m + case('Iron Flux') + CS%fe_dep_ind = m + case('NOx Flux') + CS%nox_flux_ind = m + case('NHy Flux') + CS%nhy_flux_ind = m + case('Atmospheric Pressure') + CS%atmpress_ind = m + case('xco2') + CS%xco2_ind = m + case('xco2_alt_co2') + CS%xco2_alt_ind = m + case('d14c') + CS%d14c_ind = m + case DEFAULT + write(log_message, "(A,1X,A)") & + trim(MARBL_instances%surface_flux_forcings(m)%metadata%varname), & + 'is not a valid surface flux forcing field name.' + call MOM_error(FATAL, log_message) + end select + enddo + + ! ii. store all interior forcing indices + CS%dustflux_ind = -1 + CS%PAR_col_frac_ind = -1 + CS%surf_shortwave_ind = -1 + CS%potemp_ind = -1 + CS%salinity_ind = -1 + CS%pressure_ind = -1 + CS%fesedflux_ind = -1 + CS%o2_scalef_ind = -1 + CS%remin_scalef_ind = -1 + CS%d14c_ind = -1 + allocate(CS%id_tracer_restoring(CS%ntr)) + allocate(CS%tracer_restoring_varname(CS%ntr), source=' ') ! gfortran 13.2 bug? + ! source = '' does not blank out strings + allocate(CS%tracer_restoring_ind(CS%ntr), source=-1) + allocate(CS%tracer_I_tau_ind(CS%ntr), source=-1) + CS%restore_count = 0 + I_tau_count = 0 + do m=1,size(MARBL_instances%interior_tendency_forcings) + select case (trim(MARBL_instances%interior_tendency_forcings(m)%metadata%varname)) + case('Dust Flux') + CS%dustflux_ind = m + case('PAR Column Fraction') + CS%PAR_col_frac_ind = m + case('Surface Shortwave') + CS%surf_shortwave_ind = m + case('Potential Temperature') + CS%potemp_ind = m + case('Salinity') + CS%salinity_ind = m + case('Pressure') + CS%pressure_ind = m + case('Iron Sediment Flux') + CS%fesedflux_ind = m + case('O2 Consumption Scale Factor') + CS%o2_scalef_ind = m + case('Particulate Remin Scale Factor') + CS%remin_scalef_ind = m + case DEFAULT + ! fi stands for forcing_index + fi = index(MARBL_instances%interior_tendency_forcings(m)%metadata%varname, & + 'Restoring Field') + if (fi > 0) then + CS%restore_count = CS%restore_count + 1 + CS%tracer_restoring_ind(CS%restore_count) = m + CS%tracer_restoring_varname(CS%restore_count) = & + MARBL_instances%interior_tendency_forcings(m)%metadata%varname(1:fi-2) + else + fi = index(MARBL_instances%interior_tendency_forcings(m)%metadata%varname, & + 'Restoring Inverse Timescale') + if (fi > 0) then + I_tau_count = I_tau_count + 1 + CS%tracer_I_tau_ind(I_tau_count) = m + else + write(log_message, "(A,1X,A)") & + trim(MARBL_instances%interior_tendency_forcings(m)%metadata%varname), & + 'is not a valid interior tendency forcing field name.' + call MOM_error(FATAL, log_message) + endif + endif + end select + enddo +end subroutine configure_MARBL_tracers + +!> This subroutine is used to register tracer fields and subroutines +!! to be used with MOM. +function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS, MARBL_computes_chl) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(MARBL_tracers_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module + type(tracer_registry_type), pointer :: tr_Reg !< A pointer that is set to point to the control + !! structure for the tracer advection and diffusion module. + type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct + logical, intent(out) :: MARBL_computes_chl !< If MARBL is computing chlorophyll, MOM + !! may use it to compute SW penetration + +! Local variables +! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=40) :: mdl = "MARBL_tracers" ! This module's name. + character(len=256) :: log_message + character(len=200) :: inputdir ! The directory where the input files are. + character(len=48) :: var_name ! The variable's name. + character(len=128) :: desc_name ! The variable's descriptor. + character(len=48) :: units ! The variable's units. + character(len=96) :: file_name ! file name for d14c (looped over three bands) + real, pointer :: tr_ptr(:,:,:) => NULL() + integer :: forcing_file_start_year + integer :: forcing_file_end_year + integer :: forcing_file_data_ref_year + integer :: forcing_file_model_ref_year + integer :: forcing_file_forcing_year + logical :: register_MARBL_tracers + logical :: restoring_has_edges, restoring_use_missing + logical :: restoring_timescale_has_edges, restoring_timescale_use_missing + real :: restoring_missing, restoring_timescale_missing + integer :: isd, ied, jsd, jed, nz, m, k, kbot + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + + if (associated(CS)) then + call MOM_error(WARNING, "register_MARBL_tracers called with an associated control structure.") + return + endif + allocate(CS) + + call configure_MARBL_tracers(GV, US, param_file, CS) + MARBL_computes_chl = CS%base_bio_on + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + ! ** Input directory + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + ! ** Tracer initial conditions + call get_param(param_file, mdl, "MARBL_TRACERS_IC_FILE", CS%IC_file, & + "The file in which the MARBL tracers initial values can be found.", & + default="ecosys_jan_IC_omip_latlon_1x1_180W_c230331.nc") + if (scan(CS%IC_file,'/') == 0) then + ! Add the directory if CS%IC_file is not already a complete path. + CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) + call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_IC_FILE", CS%IC_file) + endif + call get_param(param_file, mdl, "MARBL_TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code if they are not found in the "//& + "restart files. Otherwise it is a fatal error if tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + call get_param(param_file, mdl, "MARBL_TRACERS_INIT_VERTICAL_REMAP_ONLY", CS%ongrid, & + "If true, initial conditions are on the model horizontal grid. Extrapolation over " //& + "missing ocean values is done using an ICE-9 procedure with vertical ALE remapping .", & + default=.false.) + if (CS%base_bio_on) then + ! ** FESEDFLUX + call get_param(param_file, mdl, "MARBL_FESEDFLUX_FILE", CS%fesedflux_file, & + "The file in which the iron sediment flux forcing field can be found.", & + default="fesedflux_total_reduce_oxic_tx0.66v1.c230817.nc") + if (scan(CS%fesedflux_file,'/') == 0) then + ! Add the directory if CS%fesedflux_file is not already a complete path. + CS%fesedflux_file = trim(slasher(inputdir))//trim(CS%fesedflux_file) + call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FESEDFLUX_FILE", CS%fesedflux_file) + endif + ! ** FEVENTFLUX + call get_param(param_file, mdl, "MARBL_FEVENTFLUX_FILE", CS%feventflux_file, & + "The file in which the iron vent flux forcing field can be found.", & + default="feventflux_5gmol_tx0.66v1.c230817.nc") + if (scan(CS%feventflux_file,'/') == 0) then + ! Add the directory if CS%feventflux_file is not already a complete path. + CS%feventflux_file = trim(slasher(inputdir))//trim(CS%feventflux_file) + call log_param(param_file, mdl, "INPUTDIR/MARBL_TRACERS_FEVENTFLUX_FILE", CS%feventflux_file) + endif + ! ** Scale factor for FESEDFLUX + call get_param(param_file, mdl, "MARBL_FESEDFLUX_SCALE_FACTOR", CS%fesedflux_scale_factor, & + "Conversion factor between FESEDFLUX file units and MARBL units", & + units="umol m-1 d-1 -> mmol m-2 s-1", default=0.001/86400.) + + ! ** River fluxes + call get_param(param_file, mdl, "READ_RIV_FLUXES", CS%read_riv_fluxes, & + "If true, use river fluxes supplied from an input file", default=.true.) + if (CS%read_riv_fluxes) then + call get_param(param_file, mdl, "RIV_FLUX_FILE", CS%riv_flux_dataset%file_name, & + "The file in which the river fluxes can be found", & + default="riv_nut.gnews_gnm.JRA025m_to_tx0.66v1_nnsm_e333r100_190910.20210405.nc") + ! call get_param(param_file, mdl, "RIV_FLUX_OFFSET_YEAR", CS%riv) + if (scan(CS%riv_flux_dataset%file_name,'/') == 0) then + ! CS%riv_flux_dataset%file_name = trim(inputdir) // trim(CS%riv_flux_dataset%file_name) + CS%riv_flux_dataset%file_name = trim(slasher(inputdir)) //& + trim(CS%riv_flux_dataset%file_name) + call log_param(param_file, mdl, "INPUTDIR/RIV_FLUX_FILE", CS%riv_flux_dataset%file_name) + endif + call get_param(param_file, mdl, "RIV_FLUX_L_TIME_VARYING", & + CS%riv_flux_dataset%l_time_varying, & + ".true. for time-varying forcing, .false. for static forcing", default=.false.) + if (CS%riv_flux_dataset%l_time_varying) then + call get_param(param_file, mdl, "RIV_FLUX_FILE_START_YEAR", forcing_file_start_year, & + "First year of data to read in RIV_FLUX_FILE", default=1900) + call get_param(param_file, mdl, "RIV_FLUX_FILE_END_YEAR", forcing_file_end_year, & + "Last year of data to read in RIV_FLUX_FILE", default=2000) + call get_param(param_file, mdl, "RIV_FLUX_FILE_DATA_REF_YEAR", forcing_file_data_ref_year, & + "Align this year in RIV_FLUX_FILE with RIV_FLUX_FILE_MODEL_REF_YEAR in model", & + default=1900) + call get_param(param_file, mdl, "RIV_FLUX_FILE_MODEL_REF_YEAR", & + forcing_file_model_ref_year, & + "Align this year in model with RIV_FLUX_FILE_DATA_REF_YEAR in RIV_FLUX_FILE", & + default=1) + else + call get_param(param_file, mdl, "RIV_FLUX_FORCING_YEAR", forcing_file_forcing_year, & + "Year from RIV_FLUX_FILE to use for forcing", default=1900) + endif + call forcing_timeseries_set_time_type_vars(forcing_file_start_year, forcing_file_end_year, & + forcing_file_data_ref_year, forcing_file_model_ref_year, forcing_file_forcing_year, & + CS%riv_flux_dataset) + endif + endif + + if (CS%abio_dic_on) then + call get_param(param_file, mdl, "D14C_L_TIME_VARYING", CS%d14c_dataset(1)%l_time_varying, & + ".true. for time-varying forcing, .false. for static forcing", default=.false.) + CS%d14c_dataset(2)%l_time_varying = CS%d14c_dataset(1)%l_time_varying + CS%d14c_dataset(3)%l_time_varying = CS%d14c_dataset(1)%l_time_varying + if (CS%d14c_dataset(1)%l_time_varying) then + call get_param(param_file, mdl, "D14C_FILE_START_YEAR", forcing_file_start_year, & + "First year of data to read in D14C_FILE", default=1850) + call get_param(param_file, mdl, "D14C_FILE_END_YEAR", forcing_file_end_year, & + "Last year of data to read in D14C_FILE", default=2015) + call get_param(param_file, mdl, "D14C_FILE_DATA_REF_YEAR", forcing_file_data_ref_year, & + "Align this year in D14C_FILE with D14C_FILE_MODEL_REF_YEAR in model", default=1850) + call get_param(param_file, mdl, "D14C_FILE_MODEL_REF_YEAR", forcing_file_model_ref_year, & + "Align this year in model with D14C_FILE_DATA_REF_YEAR in D14C_FILE", default=1) + else + call get_param(param_file, mdl, "D14C_FORCING_YEAR", forcing_file_forcing_year, & + "Year from D14C_FILE to use for forcing", default=1850) + endif + do m=1,3 + write(var_name, "(A,I0)") "MARBL_D14C_FILE_", m + write(file_name, "(A,I0,A)") "atm_delta_C14_CMIP6_sector", m, & + "_global_1850-2015_yearly_v2.0_c240202.nc" + call get_param(param_file, mdl, var_name, CS%d14c_dataset(m)%file_name, & + "The file in which the d14c forcing field can be found.", default=file_name) + call forcing_timeseries_set_time_type_vars(forcing_file_start_year, forcing_file_end_year, & + forcing_file_data_ref_year, forcing_file_model_ref_year, forcing_file_forcing_year, & + CS%d14c_dataset(m)) + if (scan(CS%d14c_dataset(m)%file_name,'/') == 0) then + ! Add the directory if CS%d14c_dataset%file_name is not already a complete path. + CS%d14c_dataset(m)%file_name = trim(slasher(inputdir))//trim(CS%d14c_dataset(m)%file_name) + call log_param(param_file, mdl, "INPUTDIR/D14C_FILE", CS%d14c_dataset(m)%file_name) + endif + enddo +endif + + ! ** Tracer Restoring + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_SOURCE", CS%restoring_source, & + "Source of data for restoring MARBL tracers", default="none") + select case(CS%restoring_source) + case("none") + case("file") + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_FILE", CS%restoring_file, & + "File containing fields to restore MARBL tracers towards") + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_SOURCE", & + CS%restoring_I_tau_source, "Source of data for inverse timescale for restoring MARBL tracers") + + ! Initialize remapping type + call initialize_remapping(CS%restoring_remapCS, 'PCM', boundary_extrapolation=.false., answer_date=99991231) + + ! Set up array for thicknesses in restoring file + call read_Z_edges(CS%restoring_file, "PO4", CS%restoring_z_edges, CS%restoring_nz, & + restoring_has_edges, restoring_use_missing, restoring_missing, scale=US%m_to_Z) + allocate(CS%restoring_dz(CS%restoring_nz)) + do k=CS%restoring_nz,1,-1 + kbot = k + 1 ! level k is between z(k) and z(k+1) + CS%restoring_dz(k) = (CS%restoring_z_edges(k) - CS%restoring_z_edges(kbot)) * GV%Z_to_H + enddo + + select case(CS%restoring_I_tau_source) + case("file") + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_FILE", & + CS%restoring_I_tau_file, & + "File containing the inverse timescale for restoring MARBL tracers") + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_VAR_NAME", & + CS%restoring_I_tau_var_name, & + "Field containing the inverse timescale for restoring MARBL tracers", & + default="I_TAU") + ! Set up array for thicknesses in restoring timescale file + call read_Z_edges(CS%restoring_I_tau_file, CS%restoring_I_tau_var_name, CS%restoring_timescale_z_edges, & + CS%restoring_timescale_nz, restoring_timescale_has_edges, & + restoring_timescale_use_missing, restoring_timescale_missing, scale=US%m_to_Z) + allocate(CS%restoring_timescale_dz(CS%restoring_timescale_nz)) + do k=CS%restoring_timescale_nz,1,-1 + kbot = k + 1 ! level k is between z(k) and z(k+1) + CS%restoring_timescale_dz(k) = (CS%restoring_timescale_z_edges(k) - & + CS%restoring_timescale_z_edges(kbot)) * GV%Z_to_H + enddo + case DEFAULT + write(log_message, "(3A)") "'", trim(CS%restoring_I_tau_source), & + "' is not a valid option for MARBL_TRACER_RESTORING_I_TAU_SOURCE" + call MOM_error(FATAL, log_message) + end select + case DEFAULT + write(log_message, "(3A)") "'", trim(CS%restoring_source), & + "' is not a valid option for MARBL_TRACER_RESTORING_SOURCE" + call MOM_error(FATAL, log_message) + end select + + allocate(CS%ind_tr(CS%ntr)) + allocate(CS%tr_desc(CS%ntr)) + allocate(CS%tracer_data(CS%ntr)) + + do m=1,CS%ntr + allocate(CS%tracer_data(m)%tr(isd:ied,jsd:jed,nz), source=0.0) + write(var_name(:),'(A)') trim(MARBL_instances%tracer_metadata(m)%short_name) + write(desc_name(:),'(A)') trim(MARBL_instances%tracer_metadata(m)%long_name) + write(units(:),'(A)') trim(MARBL_instances%tracer_metadata(m)%units) + CS%tr_desc(m) = var_desc(trim(var_name), trim(units), trim(desc_name), caller=mdl) + + ! This is needed to force the compiler not to do a copy in the registration + ! calls. Curses on the designers and implementers of Fortran90. + tr_ptr => CS%tracer_data(m)%tr(:,:,:) + call query_vardesc(CS%tr_desc(m), name=var_name, & + caller="register_MARBL_tracers") + ! Register the tracer for horizontal advection, diffusion, and restarts. + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, units = units, & + tr_desc=CS%tr_desc(m), registry_diags=.true., & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit, & + Tr_out=CS%tracer_data(m)%tr_ptr) + + ! Set coupled_tracers to be true (hard-coded above) to provide the surface + ! values to the coupler (if any). This is meta-code and its arguments will + ! currently (deliberately) give fatal errors if it is used. + if (CS%coupled_tracers) & + CS%ind_tr(m) = aof_set_coupler_flux(trim(var_name)//'_flux', & + flux_type=' ', implementation=' ', caller="register_MARBL_tracers") + enddo + + ! Set up memory for saved state + call setup_saved_state(MARBL_instances%surface_flux_saved_state, HI, GV, restart_CS, & + CS%tracers_may_reinit, CS%surface_flux_saved_state) + call setup_saved_state(MARBL_instances%interior_tendency_saved_state, HI, GV, restart_CS, & + CS%tracers_may_reinit, CS%interior_tendency_saved_state) + + CS%tr_Reg => tr_Reg + CS%restart_CSp => restart_CS + + call set_riv_flux_tracer_inds(CS) + register_MARBL_tracers = .true. + +end function register_MARBL_tracers + +!> This subroutine initializes the CS%ntr tracer fields in tr(:,:,:,:) +!! and it sets up the tracer output. +subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag, OBC, CS, sponge_CSp) + logical, intent(in) :: restart !< .true. if the fields have already been + !! read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies + !! whether, where, and what open boundary + !! conditions are used. + type(MARBL_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_MARBL_tracers. + type(sponge_CS), pointer :: sponge_CSp !< A pointer to the control structure + !! for the sponges, if they are in use. + + ! Local variables + character(len=200) :: log_message + character(len=48) :: name ! A variable's name in a NetCDF file. + character(len=100) :: longname ! The long name of that variable. + character(len=48) :: units ! The units of the variable. + character(len=48) :: flux_units ! The units for age tracer fluxes, either + ! years m3 s-1 or years kg s-1. + character(len=48) :: tracer_name + logical :: fesedflux_has_edges, fesedflux_use_missing + real :: fesedflux_missing + integer :: i, j, k, kbot, m, diag_size + + if (.not.associated(CS)) return + if (CS%ntr < 1) return + + CS%diag => diag + + ! Allocate memory for surface tracer fluxes + allocate(CS%STF(SZI_(G), SZJ_(G), CS%ntr), & + CS%RIV_FLUXES(SZI_(G), SZJ_(G), CS%ntr), & + CS%SFO(SZI_(G), SZJ_(G), CS%sfo_cnt), & + CS%ITO(SZI_(G), SZJ_(G), SZK_(G), CS%ito_cnt), & + source=0.0) + + ! Allocate memory for d14c forcing + if (CS%abio_dic_on) allocate(CS%d14c(SZI_(G), SZJ_(G))) + + ! Register diagnostics returned from MARBL (surface flux first, then interior tendency) + call register_MARBL_diags(MARBL_instances%surface_flux_diags, diag, day, G, CS%surface_flux_diags) + call register_MARBL_diags(MARBL_instances%interior_tendency_diags, diag, day, G, & + CS%interior_tendency_diags) + + ! Register per-tracer diagnostics computed from MARBL surface flux / interior tendency values + allocate(CS%id_surface_flux_out(CS%ntr)) + allocate(CS%interior_tendency_out(CS%ntr)) + allocate(CS%interior_tendency_out_zint(CS%ntr)) + allocate(CS%interior_tendency_out_zint_100m(CS%ntr)) + do m=1,CS%ntr + write(name, "(2A)") "STF_", trim(MARBL_instances%tracer_metadata(m)%short_name) + write(longname, "(2A)") trim(MARBL_instances%tracer_metadata(m)%long_name), " Surface Flux" + write(units, "(2A)") trim(MARBL_instances%tracer_metadata(m)%units), " m/s" + CS%id_surface_flux_out(m) = register_diag_field("ocean_model", trim(name), & + diag%axesT1, & ! T => tracer grid? 1 => no vertical grid + day, trim(longname), trim(units), conversion=US%Z_to_m*US%s_to_T) + + write(name, "(2A)") "J_", trim(MARBL_instances%tracer_metadata(m)%short_name) + write(longname, "(2A)") trim(MARBL_instances%tracer_metadata(m)%long_name), " Source Sink Term" + write(units, "(2A)") trim(MARBL_instances%tracer_metadata(m)%units), "/s" + CS%interior_tendency_out(m)%id = register_diag_field("ocean_model", trim(name), & + diag%axesTL, & ! T=> tracer grid? L => layer center + day, trim(longname), trim(units)) + if (CS%interior_tendency_out(m)%id > 0) & + allocate(CS%interior_tendency_out(m)%field_3d(SZI_(G),SZJ_(G), SZK_(G)), source=0.0) + + write(name, "(2A)") "Jint_", trim(MARBL_instances%tracer_metadata(m)%short_name) + write(longname, "(2A)") trim(MARBL_instances%tracer_metadata(m)%long_name), & + " Source Sink Term Vertical Integral" + write(units, "(2A)") trim(MARBL_instances%tracer_metadata(m)%units), " m/s" + CS%interior_tendency_out_zint(m)%id = register_diag_field("ocean_model", trim(name), & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, trim(longname), trim(units)) + if (CS%interior_tendency_out_zint(m)%id > 0) & + allocate(CS%interior_tendency_out_zint(m)%field_2d(SZI_(G),SZJ_(G)), source=0.0) + + write(name, "(2A)") "Jint_100m_", trim(MARBL_instances%tracer_metadata(m)%short_name) + write(longname, "(2A)") trim(MARBL_instances%tracer_metadata(m)%long_name), & + " Source Sink Term Vertical Integral, 0-100m" + write(units, "(2A)") trim(MARBL_instances%tracer_metadata(m)%units), " m/s" + CS%interior_tendency_out_zint_100m(m)%id = register_diag_field("ocean_model", trim(name), & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, trim(longname), trim(units)) + if (CS%interior_tendency_out_zint_100m(m)%id > 0) & + allocate(CS%interior_tendency_out_zint_100m(m)%field_2d(SZI_(G),SZJ_(G)), source=0.0) + + enddo + + ! Register diagnostics for MOM to report that are not tracer specific + CS%bot_flux_to_tend_id = register_diag_field("ocean_model", "BOT_FLUX_TO_TEND", & + diag%axesTL, & ! T=> tracer grid? L => layer center + day, "Conversion Factor for Bottom Flux -> Tend", "1/m") + + do m=1,CS%ntr + call query_vardesc(CS%tr_desc(m), name=name, caller="initialize_MARBL_tracers") + if ((.not. restart) .or. & + (CS%tracers_may_reinit .and. & + .not. query_initialized(CS%tracer_data(m)%tr(:,:,:), name, CS%restart_CSp))) then + ! TODO: added the ongrid optional argument, but is there a good way to detect if the file is on grid? + call MOM_initialize_tracer_from_Z(h, CS%tracer_data(m)%tr, G, GV, US, param_file, & + CS%IC_file, name, ongrid=CS%ongrid) + do k=1,GV%ke + do j=G%jsc, G%jec + do i=G%isc, G%iec + ! Ensure tracer concentrations are at / above minimum value + if (CS%tracer_data(m)%tr(i,j,k) < CS%IC_min) CS%tracer_data(m)%tr(i,j,k) = CS%IC_min + enddo + enddo + enddo + endif + enddo + + ! Register diagnostics for river fluxes + CS%no3_riv_flux = register_diag_field("ocean_model", "NO3_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Inorganic Nitrate Riverine Flux", "mmol/m^3 m/s") + CS%po4_riv_flux = register_diag_field("ocean_model", "PO4_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Inorganic Phosphate Riverine Flux", "mmol/m^3 m/s") + CS%don_riv_flux = register_diag_field("ocean_model", "DON_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Organic Nitrogen Riverine Flux", "mmol/m^3 m/s") + CS%donr_riv_flux = register_diag_field("ocean_model", "DONR_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Refractory DON Riverine Flux", "mmol/m^3 m/s") + CS%dop_riv_flux = register_diag_field("ocean_model", "DOP_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Organic Phosphorus Riverine Flux", "mmol/m^3 m/s") + CS%dopr_riv_flux = register_diag_field("ocean_model", "DOPR_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Refractory DOP Riverine Flux", "mmol/m^3 m/s") + CS%sio3_riv_flux = register_diag_field("ocean_model", "SiO3_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Inorganic Silicate Riverine Flux", "mmol/m^3 m/s") + CS%fe_riv_flux = register_diag_field("ocean_model", "Fe_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Inorganic Iron Riverine Flux", "mmol/m^3 m/s") + CS%doc_riv_flux = register_diag_field("ocean_model", "DOC_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Organic Carbon Riverine Flux", "mmol/m^3 m/s") + CS%docr_riv_flux = register_diag_field("ocean_model", "DOCR_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Refractory DOC Riverine Flux", "mmol/m^3 m/s") + CS%alk_riv_flux = register_diag_field("ocean_model", "ALK_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Alkalinity Riverine Flux", "meq/m^3 m/s") + CS%alk_alt_co2_riv_flux = register_diag_field("ocean_model", "ALK_ALT_CO2_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Alkalinity Riverine Flux, Alternative CO2", "meq/m^3 m/s") + CS%dic_riv_flux = register_diag_field("ocean_model", "DIC_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Inorganic Carbon Riverine Flux", "mmol/m^3 m/s") + CS%dic_alt_co2_riv_flux = register_diag_field("ocean_model", "DIC_ALT_CO2_RIV_FLUX", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Dissolved Inorganic Carbon Riverine Flux, Alternative CO2", "mmol/m^3 m/s") + + ! Register diagnostics for d14c forcing + if (CS%abio_dic_on) then + CS%d14c_id = register_diag_field("ocean_model", "D14C_FORCING", & + diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid + day, "Delta-14C in atmospheric CO2", "per mil, relative to Modern") + endif + + ! Register diagnostics for per-category forcing fields + if (CS%ice_ncat > 0) then + allocate(CS%fracr_cat_id(CS%ice_ncat+1)) + allocate(CS%qsw_cat_id(CS%ice_ncat+1)) + do m=1,CS%ice_ncat+1 + write(name, "(A,I0)") "FRACR_CAT_", m + write(longname, "(A,I0)") "Fraction of area in ice category ", m + units = "fraction" + CS%fracr_cat_id(m) = register_diag_field("ocean_model", trim(name), & + diag%axesT1, & ! T => tracer grid? 1 => no vertical grid + day, trim(longname), trim(units)) + write(name, "(A,I0)") "QSW_CAT_", m + write(longname, "(A,I0)") "Shortwave penetrating through ice category ", m + units = "TODO: set units" + CS%qsw_cat_id(m) = register_diag_field("ocean_model", trim(name), & + diag%axesT1, & ! T => tracer grid? 1 => no vertical grid + day, trim(longname), trim(units)) + enddo + endif + + if (CS%base_bio_on) then + ! Read initial fesedflux and feventflux fields + ! (1) get vertical dimension + ! -- comes from fesedflux_file, assume same dimension in feventflux + ! (maybe these fields should be combined?) + ! -- note: read_Z_edges treats depth as positive UP => 0 at surface, negative at depth + fesedflux_use_missing = .false. + call read_Z_edges(CS%fesedflux_file, "FESEDFLUXIN", CS%fesedflux_z_edges, CS%fesedflux_nz, & + fesedflux_has_edges, fesedflux_use_missing, fesedflux_missing, scale=US%m_to_Z) + + ! (2) Allocate memory for fesedflux and feventflux + allocate(CS%fesedflux_in(SZI_(G), SZJ_(G), CS%fesedflux_nz)) + allocate(CS%feventflux_in(SZI_(G), SZJ_(G), CS%fesedflux_nz)) + allocate(CS%fesedflux_dz(SZI_(G), SZJ_(G), CS%fesedflux_nz)) + + ! (3) Read data + ! TODO: Add US term to scale + call MOM_read_data(CS%fesedflux_file, "FESEDFLUXIN", CS%fesedflux_in(:,:,:), G%Domain, & + scale=CS%fesedflux_scale_factor) + call MOM_read_data(CS%feventflux_file, "FESEDFLUXIN", CS%feventflux_in(:,:,:), G%Domain, & + scale=CS%fesedflux_scale_factor) + + ! (4) Relocate values that are below ocean bottom to layer that intersects bathymetry + ! Remember, fesedflux_z_edges = 0 at surface and is < 0 below surface + + do k=CS%fesedflux_nz, 1, -1 + kbot = k + 1 ! level k is between z(k) and z(k+1) + do j=G%jsc, G%jec + do i=G%isc, G%iec + if (G%mask2dT(i,j) == 0) cycle + if (G%bathyT(i,j) + CS%fesedflux_z_edges(1) < 1e-8 * US%m_to_Z) then + write(log_message, *) "Current implementation of fesedflux assumes G%bathyT >=", & + " first edge;first edge = ", -CS%fesedflux_z_edges(1), "bathyT = ", G%bathyT(i,j) + call MOM_error(FATAL, log_message) + endif + ! Also figure out layer thickness while we're here + CS%fesedflux_dz(i,j,k) = (CS%fesedflux_z_edges(k) - CS%fesedflux_z_edges(kbot)) * GV%Z_to_H + ! If top interface is at or below ocean bottom, move flux in current layer up one + ! and set thickness of current level to 0 + if (G%bathyT(i,j) + CS%fesedflux_z_edges(k) < 1e-8 * US%m_to_Z) then + CS%fesedflux_in(i,j,k-1) = CS%fesedflux_in(i,j,k-1) + CS%fesedflux_in(i,j,k) + CS%fesedflux_in(i,j,k) = 0. + CS%feventflux_in(i,j,k-1) = CS%feventflux_in(i,j,k-1) + CS%feventflux_in(i,j,k) + CS%feventflux_in(i,j,k) = 0. + CS%fesedflux_dz(i,j,k) = 0. + elseif (G%bathyT(i,j) + CS%fesedflux_z_edges(kbot) < 1e-8 * US%m_to_Z) then + ! Otherwise, if lower interface is below bathymetry move interface to ocean bottom + CS%fesedflux_dz(i,j,k) = (G%bathyT(i,j) + CS%fesedflux_z_edges(k)) * GV%Z_to_H + endif + enddo + enddo + enddo + + ! Initialize external field for river fluxes + if (CS%read_riv_fluxes) then + CS%id_din_riv = init_external_field(CS%riv_flux_dataset%file_name, 'din_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_don_riv = init_external_field(CS%riv_flux_dataset%file_name, 'don_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_dip_riv = init_external_field(CS%riv_flux_dataset%file_name, 'dip_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_dop_riv = init_external_field(CS%riv_flux_dataset%file_name, 'dop_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_dsi_riv = init_external_field(CS%riv_flux_dataset%file_name, 'dsi_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_dfe_riv = init_external_field(CS%riv_flux_dataset%file_name, 'dfe_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_dic_riv = init_external_field(CS%riv_flux_dataset%file_name, 'dic_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_alk_riv = init_external_field(CS%riv_flux_dataset%file_name, 'alk_riv_flux', & + domain=G%Domain%mpp_domain) + CS%id_doc_riv = init_external_field(CS%riv_flux_dataset%file_name, 'doc_riv_flux', & + domain=G%Domain%mpp_domain) + endif + endif + + if (CS%abio_dic_on) then + ! Initialize external field for d14c forcing + do m=1,3 + CS%id_d14c(m) = init_external_field(CS%d14c_dataset(m)%file_name, "Delta14co2_in_air", & + ignore_axis_atts=.true.) + enddo + endif + + ! Initialize external field for restoring + if (CS%restoring_I_tau_source == "file") then + select case(CS%restoring_source) + case("file") + ! Set up array for reading in raw restoring data + allocate(CS%restoring_in(SZI_(G), SZJ_(G), CS%restoring_nz, CS%restore_count), source=0.) + do m=1,CS%restore_count + CS%id_tracer_restoring(m) = init_external_field(CS%restoring_file, & + trim(CS%tracer_restoring_varname(m)), domain=G%Domain%mpp_domain) + enddo + end select + select case(CS%restoring_I_tau_source) + case("file") + allocate(CS%I_tau(SZI_(G), SZJ_(G), CS%restoring_timescale_nz), source=0.) + call MOM_read_data(CS%restoring_I_tau_file, "RTAU", CS%I_tau(:,:,:), G%Domain) + end select + endif + +end subroutine initialize_MARBL_tracers + +!> This subroutine is used to register tracer fields and subroutines +!! to be used with MOM. +subroutine register_MARBL_diags(MARBL_diags, diag, day, G, id_diags) + + type(marbl_diagnostics_type), intent(in) :: MARBL_diags !< MARBL diagnostics from MARBL_instances + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(diag_ctrl), target, intent(in) :: diag !< Structure used to regulate diagnostic output. + !integer, allocatable, intent(inout) :: id_diags(:) !< allocatable array storing diagnostic index number + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(temp_marbl_diag), allocatable, intent(inout) :: id_diags(:) !< allocatable array storing diagnostic index + !! number and buffer space for collecting diags + !! from all columns + + integer :: m, diag_size + + diag_size = size(MARBL_diags%diags) + allocate(id_diags(diag_size)) + do m = 1, diag_size + id_diags(m)%id = -1 + if (trim(MARBL_diags%diags(m)%vertical_grid) .eq. "none") then ! 2D field + id_diags(m)%id = register_diag_field("ocean_model", & + trim(MARBL_diags%diags(m)%short_name), & + diag%axesT1, & ! T => tracer grid? 1 => no vertical grid + day, & + trim(MARBL_diags%diags(m)%long_name), & + trim(MARBL_diags%diags(m)%units)) + if (id_diags(m)%id > 0) allocate(id_diags(m)%field_2d(SZI_(G),SZJ_(G)), source=0.0) + else ! 3D field + ! TODO: MARBL should provide v_extensive through MARBL_diags + ! (for now, FESEDFLUX is the only one that should be true) + ! Also, known issue where passing v_extensive=.false. isn't + ! treated the same as not passing v_extensive + if (trim(MARBL_diags%diags(m)%short_name).eq."FESEDFLUX") then + id_diags(m)%id = register_diag_field("ocean_model", & + trim(MARBL_diags%diags(m)%short_name), & + diag%axesTL, & ! T=> tracer grid? L => layer center + day, & + trim(MARBL_diags%diags(m)%long_name), & + trim(MARBL_diags%diags(m)%units), & + v_extensive=.true.) + else + id_diags(m)%id = register_diag_field("ocean_model", & + trim(MARBL_diags%diags(m)%short_name), & + diag%axesTL, & ! T=> tracer grid? L => layer center + day, & + trim(MARBL_diags%diags(m)%long_name), & + trim(MARBL_diags%diags(m)%units)) + endif + if (id_diags(m)%id > 0) allocate(id_diags(m)%field_3d(SZI_(G),SZJ_(G), SZK_(G)), source=0.0) + endif + enddo + +end subroutine register_MARBL_diags + +!> This subroutine allocates memory for saved state fields and registers them in the restart files +subroutine setup_saved_state(MARBL_saved_state, HI, GV, restart_CS, tracers_may_reinit, & + local_saved_state) + + type(marbl_saved_state_type), intent(in) :: MARBL_saved_state !< MARBL saved state from + !! MARBL_instances + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(MOM_restart_CS), pointer, intent(in) :: restart_CS !< control structure to add saved state + !! to restarts + logical, intent(in) :: tracers_may_reinit !< used to determine mandatory + !! flag in restart + type(saved_state_for_MARBL_type), allocatable, intent(inout) :: local_saved_state(:) !< allocatable array for local + !! saved state + + integer :: num_fields, m + character(len=200) :: log_message, varname + + num_fields = MARBL_saved_state%saved_state_cnt + allocate(local_saved_state(num_fields)) + + do m=1,num_fields + write(varname, "(2A)") "MARBL_", trim(MARBL_saved_state%state(m)%short_name) + select case (MARBL_saved_state%state(m)%rank) + case (2) + allocate(local_saved_state(m)%field_2d(SZI_(HI),SZJ_(HI)), source=0.0) + call register_restart_field(local_saved_state(m)%field_2d, varname, & + .not.tracers_may_reinit, restart_CS) + case (3) + if (trim(MARBL_saved_state%state(m)%vertical_grid).eq."layer_avg") then + allocate(local_saved_state(m)%field_3d(SZI_(HI),SZJ_(HI), SZK_(GV)), source=0.0) + call register_restart_field(local_saved_state(m)%field_3d, varname, & + .not.tracers_may_reinit, restart_CS) + else + write(log_message, "(3A, I0, A)") "'", trim(MARBL_saved_state%state(m)%vertical_grid), & + "' is an invalid vertical grid for saved state (ind = ", m, ")" + call MOM_error(FATAL, log_message) + endif + case DEFAULT + write(log_message, "(I0, A, I0, A)") MARBL_saved_state%state(m)%rank, & + " is an invalid rank for saved state (ind = ", m, ")" + call MOM_error(FATAL, log_message) + end select + local_saved_state(m)%short_name = trim(MARBL_saved_state%state(m)%short_name) + write(local_saved_state(m)%file_varname, "(2A)") "MARBL_", trim(local_saved_state(m)%short_name) + local_saved_state(m)%units = trim(MARBL_saved_state%state(m)%units) + enddo + +end subroutine setup_saved_state + +!> This subroutine applies diapycnal diffusion and any other column +!! tracer physics or chemistry to the tracers from this file. +subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, tv, & + KPP_CSp, nonLocalTrans, evap_CFL_limit, minimum_forcing_depth) + + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: ea !< an array to which the amount of fluid entrained + !! from the layer above during this call will be + !! added [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: eb !< an array to which the amount of fluid entrained + !! from the layer below during this call will be + !! added [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic + !! and tracer forcing fields. Unused fields have NULL ptrs. + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(MARBL_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_MARBL_tracers. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + type(KPP_CS), optional, pointer :: KPP_CSp !< KPP control structure + real, optional, intent(in) :: nonLocalTrans(:,:,:) !< Non-local transport [nondim] + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which + !! fluxes can be applied [m] + +! Local variables + character(len=256) :: log_message + real, dimension(SZI_(G),SZJ_(G)) :: ref_mask ! Mask for 2D MARBL diags using ref_depth + real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_loc ! Local copy of CS%RIV_FLUXES*dt + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: bot_flux_to_tend + real :: cum_bftt_dz ! sum of bot_flux_to_tend * dz from the bottom layer to current layer + real :: sfc_val ! The surface value for the tracers. + real :: Isecs_per_year ! The number of seconds in a year. + real :: year ! The time in years. + integer :: secs, days ! Integer components of the time type. + real, dimension(0:GV%ke) :: zi ! z-coordinate interface depth [Z ~> m] + real, dimension(GV%ke) :: zc ! z-coordinate layer center depth [Z ~> m] + real, dimension(GV%ke) :: dz ! z-coordinate cell thickness [H ~> m] + integer :: i, j, k, is, ie, js, je, nz, m + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (.not.associated(CS)) return + + ! (1) Compute surface fluxes + ! FIXME: MARBL can handle computing surface fluxes for all columns simultaneously + ! I was just thinking going column-by-column at first might be easier + do j=js,je + do i=is,ie + ! i. only want ocean points in this loop + if (G%mask2dT(i,j) == 0) cycle + + ! ii. Load proper column data + ! * surface flux forcings + ! These fields are getting the correct data + ! TODO: if top layer is vanishly thin, do we actually want (e.g.) top 5m average temp / salinity? + ! How does MOM pass SST and SSS to GFDL coupler? (look in core.F90?) + if (CS%sss_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%sss_ind)%field_0d(1) = tv%S(i,j,1) * US%S_to_ppt + if (CS%sst_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%sst_ind)%field_0d(1) = tv%T(i,j,1) * US%C_to_degC + if (CS%ifrac_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%ifrac_ind)%field_0d(1) = fluxes%ice_fraction(i,j) + + ! MARBL wants u10_sqr in (m/s)^2 + if (CS%u10_sqr_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%u10_sqr_ind)%field_0d(1) = fluxes%u10_sqr(i,j) * & + ((US%L_T_to_m_s)**2) + + ! mct_driver/ocn_cap_methods:93 -- ice_ocean_boundary%p(i,j) comes from coupler + ! We may need a new ice_ocean_boundary%p_atm because %p includes ice in GFDL driver + if (CS%atmpress_ind > 0) then + if (associated(fluxes%p_surf_full)) then + MARBL_instances%surface_flux_forcings(CS%atmpress_ind)%field_0d(1) = & + fluxes%p_surf_full(i,j) * ((US%R_to_kg_m3 * (US%L_T_to_m_s**2)) * atm_per_Pa) + else + ! hardcode value of 1 atm (can't figure out how to get this from solo_driver) + MARBL_instances%surface_flux_forcings(CS%atmpress_ind)%field_0d(1) = 1. + endif + endif + + ! These are okay, but need option to come in from coupler + if (CS%xco2_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%xco2_ind)%field_0d(1) = fluxes%atm_co2(i,j) + if (CS%xco2_alt_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%xco2_alt_ind)%field_0d(1) = fluxes%atm_alt_co2(i,j) + + ! These are okay, but need option to read in from file + if (CS%dust_dep_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%dust_dep_ind)%field_0d(1) = & + fluxes%dust_flux(i,j) * US%RZ_T_to_kg_m2s + + if (CS%fe_dep_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%fe_dep_ind)%field_0d(1) = & + fluxes%iron_flux(i,j) * (US%Z_to_m * US%s_to_T) + + ! MARBL wants ndep in (mmol/m^2/s) + if (CS%nox_flux_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%nox_flux_ind)%field_0d(1) = fluxes%noy_dep(i,j) * & + (US%Z_to_m * US%s_to_T) + if (CS%nhy_flux_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%nhy_flux_ind)%field_0d(1) = fluxes%nhx_dep(i,j) * & + (US%Z_to_m * US%s_to_T) + + if (CS%d14c_ind > 0) & + MARBL_instances%surface_flux_forcings(CS%d14c_ind)%field_0d(1) = CS%d14c(i,j) + + ! * tracers at surface + ! TODO: average over some shallow depth (e.g. 5m) + do m=1,CS%ntr + MARBL_instances%tracers_at_surface(1,m) = CS%tracer_data(m)%tr(i,j,1) + enddo + + ! * surface flux saved state + do m=1,size(MARBL_instances%surface_flux_saved_state%state) + ! (currently only 2D fields are saved from surface_flux_compute()) + MARBL_instances%surface_flux_saved_state%state(m)%field_2d(1) = & + CS%surface_flux_saved_state(m)%field_2d(i,j) + enddo + + ! iii. Compute surface fluxes in MARBL + call MARBL_instances%surface_flux_compute() + if (MARBL_instances%StatusLog%labort_marbl) then + call MARBL_instances%StatusLog%log_error_trace("MARBL_instances%surface_flux_compute()", & + "MARBL_tracers_column_physics") + endif + call print_marbl_log(MARBL_instances%StatusLog) + call MARBL_instances%StatusLog%erase() + + ! iv. Copy output that MOM6 needs to hold on to + ! * saved state + do m=1,size(MARBL_instances%surface_flux_saved_state%state) + CS%surface_flux_saved_state(m)%field_2d(i,j) = & + MARBL_instances%surface_flux_saved_state%state(m)%field_2d(1) + enddo + + ! * diagnostics + do m=1,size(MARBL_instances%surface_flux_diags%diags) + ! All diags are 2D coming from surface + if (CS%surface_flux_diags(m)%id > 0) & + CS%surface_flux_diags(m)%field_2d(i,j) = & + real(MARBL_instances%surface_flux_diags%diags(m)%field_2d(1)) + enddo + + ! * Surface tracer flux + CS%STF(i,j,:) = MARBL_instances%surface_fluxes(1,:) * (US%m_to_Z * US%T_to_s) + + ! * Surface flux output + do m=1,CS%sfo_cnt + CS%SFO(i,j,m) = MARBL_instances%surface_flux_output%outputs_for_GCM(m)%forcing_field_0d(1) + enddo + + enddo + enddo + + if (CS%debug) then + do m=1,CS%ntr + call hchksum(CS%STF(:,:,m), & + trim(MARBL_instances%tracer_metadata(m)%short_name)//" sfc_flux", G%HI, & + scale=US%Z_to_m*US%s_to_T) + enddo + endif + + ! (2) Post surface fluxes and their diagnostics (currently all 2D) + do m=1,CS%ntr + if (CS%id_surface_flux_out(m) > 0) & + call post_data(CS%id_surface_flux_out(m), CS%STF(:,:,m), CS%diag) + enddo + do m=1,size(CS%surface_flux_diags) + if (CS%surface_flux_diags(m)%id > 0) & + call post_data(CS%surface_flux_diags(m)%id, CS%surface_flux_diags(m)%field_2d(:,:), CS%diag) + enddo + + ! (3) Apply surface fluxes via vertical diffusion + ! Compute KPP nonlocal term if necessary + if (present(KPP_CSp)) then + if (associated(KPP_CSp) .and. present(nonLocalTrans)) then + do m=1,CS%ntr + call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, CS%STF(:,:,m), dt, & + CS%diag, CS%tracer_data(m)%tr_ptr, CS%tracer_data(m)%tr(:,:,:), & + flux_scale=GV%Z_to_H) + enddo + endif + if (CS%debug) then + do m=1,CS%ntr + call hchksum(CS%tracer_data(m)%tr(:,:,m), & + trim(MARBL_instances%tracer_metadata(m)%short_name)//' post KPP', G%HI) + enddo + endif + endif + + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + do m=1,CS%ntr + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + ! CS%RIV_FLUXES is conc m/s, in_flux_optional expects time-integrated flux (conc H) + do j=js,je ; do i=is,ie + riv_flux_loc(i,j) = (CS%RIV_FLUXES(i,j,m) * (dt*US%T_to_s)) * GV%m_to_H + enddo ; enddo + if (CS%debug) & + call hchksum(riv_flux_loc(:,:), & + trim(MARBL_instances%tracer_metadata(m)%short_name)//' riv flux', G%HI, scale=GV%H_to_m) + call applyTracerBoundaryFluxesInOut(G, GV, CS%tracer_data(m)%tr(:,:,:) , dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth, in_flux_optional=riv_flux_loc) + call tracer_vertdiff(h_work, ea, eb, dt, CS%tracer_data(m)%tr(:,:,:), G, GV, & + sfc_flux=GV%Rho0 * CS%STF(:,:,m)) + enddo + else + do m=1,CS%ntr + call tracer_vertdiff(h_old, ea, eb, dt, CS%tracer_data(m)%tr(:,:,:), G, GV, & + sfc_flux=GV%Rho0 * CS%STF(:,:,m)) + enddo + endif + + if (CS%debug) then + do m=1,CS%ntr + call hchksum(CS%tracer_data(m)%tr(:,:,m), & + trim(MARBL_instances%tracer_metadata(m)%short_name)//' post tracer_vertdiff', G%HI) + enddo + endif + + ! (4) Compute interior tendencies + + bot_flux_to_tend(:, :, :) = 0. + do j=js,je + do i=is,ie + ! i. only want ocean points in this loop + if (G%mask2dT(i,j) == 0) cycle + + ! ii. Set up vertical domain and bot_flux_to_tend + ! Calculate depth of interface by building up thicknesses from the bottom (top interface is always 0) + ! MARBL wants this to be positive-down + zi(GV%ke) = G%bathyT(i,j) + MARBL_instances%bot_flux_to_tend(:) = 0. + cum_bftt_dz = 0. + do k = GV%ke, 1, -1 + ! TODO: if we move this above vertical mixing, use h_old + dz(k) = h_new(i,j,k) ! cell thickness + zc(k) = zi(k) - 0.5 * (dz(k)*GV%H_to_Z) + zi(k-1) = zi(k) - (dz(k)*GV%H_to_Z) + if (G%bathyT(i,j) - zi(k-1) <= CS%bot_flux_mix_thickness) then + MARBL_instances%bot_flux_to_tend(k) = US%m_to_Z * CS%Ibfmt + cum_bftt_dz = cum_bftt_dz + MARBL_instances%bot_flux_to_tend(k) * (GV%H_to_m * dz(k)) + elseif (G%bathyT(i,j) - zi(k) < CS%bot_flux_mix_thickness) then + ! MARBL_instances%bot_flux_to_tend(k) = (1. - (G%bathyT(i,j) - zi(k)) * CS%Ibfmt) / dz(k) + MARBL_instances%bot_flux_to_tend(k) = (1. - cum_bftt_dz) / (GV%H_to_m * dz(k)) + endif + enddo + if (G%bathyT(i,j) - zi(0) < CS%bot_flux_mix_thickness) & + MARBL_instances%bot_flux_to_tend(:) = MARBL_instances%bot_flux_to_tend(:) * & + CS%bot_flux_mix_thickness / (G%bathyT(i,j) - zi(0)) + if (CS%bot_flux_to_tend_id > 0) & + bot_flux_to_tend(i, j, :) = MARBL_instances%bot_flux_to_tend(:) + + ! zw(1:nz) is bottom cell depth so no element of zw = 0, it is assumed to be top layer depth + MARBL_instances%domain%zw(:) = US%Z_to_m * zi(1:GV%ke) + MARBL_instances%domain%zt(:) = US%Z_to_m * zc(:) + MARBL_instances%domain%delta_z(:) = GV%H_to_m * dz(:) + + ! iii. Load proper column data + ! * Forcing Fields + ! These fields are getting the correct data + if (CS%potemp_ind > 0) & + MARBL_instances%interior_tendency_forcings(CS%potemp_ind)%field_1d(1,:) = tv%T(i,j,:) * US%C_to_degC + if (CS%salinity_ind > 0) & + MARBL_instances%interior_tendency_forcings(CS%salinity_ind)%field_1d(1,:) = tv%S(i,j,:) * US%S_to_ppt + + ! This are okay, but need option to read in from file + ! (Same as dust_dep_ind for surface_flux_forcings) + if (CS%dustflux_ind > 0) & + MARBL_instances%interior_tendency_forcings(CS%dustflux_ind)%field_0d(1) = & + fluxes%dust_flux(i,j) * US%RZ_T_to_kg_m2s + + ! TODO: Support PAR (currently just using single subcolumn) + ! (Look for Pen_sw_bnd?) + if (CS%PAR_col_frac_ind > 0) then + ! second index is num_subcols, not depth + !MARBL_instances%interior_tendency_forcings(CS%PAR_col_frac_ind)%field_1d(1,:) = fluxes%fracr_cat(i,j,:) + if (CS%use_ice_category_fields) then + MARBL_instances%interior_tendency_forcings(CS%PAR_col_frac_ind)%field_1d(1,:) = & + fluxes%fracr_cat(i,j,:) + else + MARBL_instances%interior_tendency_forcings(CS%PAR_col_frac_ind)%field_1d(1,1) = 1. + endif + endif + + if (CS%surf_shortwave_ind > 0) then + ! second index is num_subcols, not depth + if (CS%use_ice_category_fields) then + MARBL_instances%interior_tendency_forcings(CS%surf_shortwave_ind)%field_1d(1,:) = & + fluxes%qsw_cat(i,j,:) + else + MARBL_instances%interior_tendency_forcings(CS%surf_shortwave_ind)%field_1d(1,1) = & + fluxes%sw(i,j) * US%QRZ_T_to_W_m2 + endif + endif + ! Tracer restoring + do m=1,CS%restore_count + MARBL_instances%interior_tendency_forcings(CS%tracer_restoring_ind(m))%field_1d(1,:) = 0. + call remapping_core_h(CS%restoring_remapCS, CS%restoring_nz, CS%restoring_dz(:), & + CS%restoring_in(i,j,:,m), GV%ke, dz(:), & + MARBL_instances%interior_tendency_forcings(CS%tracer_restoring_ind(m))%field_1d(1,:)) + if (m==1) then + call remapping_core_h(CS%restoring_remapCS, CS%restoring_timescale_nz, & + CS%restoring_timescale_dz(:), CS%I_tau(i,j,:), GV%ke, dz(:), & + MARBL_instances%interior_tendency_forcings(CS%tracer_I_tau_ind(m))%field_1d(1,:)) + else + MARBL_instances%interior_tendency_forcings(CS%tracer_I_tau_ind(m))%field_1d(1,:) = & + MARBL_instances%interior_tendency_forcings(CS%tracer_I_tau_ind(1))%field_1d(1,:) + endif + enddo + + ! TODO: In POP, pressure comes from a function in state_mod.F90; I don't see a similar function here + ! This formulation is from Levitus 1994, and I think it belongs in MOM_EOS.F90? + ! Converts depth [m] -> pressure [bars] + ! NOTE: Andrew recommends using GV%H_to_Pa + if (CS%pressure_ind > 0) & + MARBL_instances%interior_tendency_forcings(CS%pressure_ind)%field_1d(1,:) = & + (0.0598088 * (exp(-0.025*US%Z_to_m * zc(:)) - 1.)) + & + (0.100766 * US%Z_to_m * zc(:)) + (2.28405e-7*((US%Z_to_m * zc(:))**2)) + + if (CS%fesedflux_ind > 0) then + MARBL_instances%interior_tendency_forcings(CS%fesedflux_ind)%field_1d(1,:) = 0. + call reintegrate_column(CS%fesedflux_nz, & + CS%fesedflux_dz(i,j,:) * (sum(dz(:) * GV%H_to_Z) / G%bathyT(i,j)), & + CS%fesedflux_in(i,j,:) + CS%feventflux_in(i,j,:), GV%ke, dz(:), & + MARBL_instances%interior_tendency_forcings(CS%fesedflux_ind)%field_1d(1,:)) + endif + + ! TODO: add ability to read these fields from file + ! also, add constant values to CS + if (CS%o2_scalef_ind > 0) & + MARBL_instances%interior_tendency_forcings(CS%o2_scalef_ind)%field_1d(1,:) = 1. + if (CS%remin_scalef_ind > 0) & + MARBL_instances%interior_tendency_forcings(CS%remin_scalef_ind)%field_1d(1,:) = 1. + + ! * Column Tracers + do m=1,CS%ntr + MARBL_instances%tracers(m, :) = CS%tracer_data(m)%tr(i,j,:) + enddo + + ! * interior tendency saved state + ! (currently only 3D fields are saved from interior_tendency_compute()) + do m=1,size(MARBL_instances%interior_tendency_saved_state%state) + MARBL_instances%interior_tendency_saved_state%state(m)%field_3d(:,1) = & + CS%interior_tendency_saved_state(m)%field_3d(i,j,:) + enddo + + ! iv. Compute interior tendencies in MARBL + call MARBL_instances%interior_tendency_compute() + if (MARBL_instances%StatusLog%labort_marbl) then + call MARBL_instances%StatusLog%log_error_trace(& + "MARBL_instances%interior_tendency_compute()", "MARBL_tracers_column_physics") + endif + call print_marbl_log(MARBL_instances%StatusLog, G, i, j) + call MARBL_instances%StatusLog%erase() + + ! v. Apply tendencies immediately + ! First pass - Euler step; if stability issues, we can do something different (subcycle?) + do m=1,CS%ntr + CS%tracer_data(m)%tr(i,j,:) = CS%tracer_data(m)%tr(i,j,:) + (dt * US%T_to_s) * & + MARBL_instances%interior_tendencies(m,:) + enddo + + ! vi. Copy output that MOM6 needs to hold on to + ! * saved state + do m=1,size(MARBL_instances%interior_tendency_saved_state%state) + CS%interior_tendency_saved_state(m)%field_3d(i,j,:) = & + MARBL_instances%interior_tendency_saved_state%state(m)%field_3d(:,1) + enddo + + ! * diagnostics + do m=1,size(MARBL_instances%interior_tendency_diags%diags) + if (CS%interior_tendency_diags(m)%id > 0) then + if (allocated(CS%interior_tendency_diags(m)%field_2d)) then + ! Only copy values if ref_depth < bathyT + if (G%bathyT(i,j) > real(MARBL_instances%interior_tendency_diags%diags(m)%ref_depth)) then + CS%interior_tendency_diags(m)%field_2d(i,j) = & + real(MARBL_instances%interior_tendency_diags%diags(m)%field_2d(1)) + endif + else ! not a 2D diagnostic + CS%interior_tendency_diags(m)%field_3d(i,j,:) = & + real(MARBL_instances%interior_tendency_diags%diags(m)%field_3d(:,1)) + endif + endif + enddo + + ! * tendency values themselves (and vertical integrals of them) + do m=1,CS%ntr + if (allocated(CS%interior_tendency_out(m)%field_3d)) & + CS%interior_tendency_out(m)%field_3d(i,j,:) = MARBL_instances%interior_tendencies(m,:) + + if (allocated(CS%interior_tendency_out_zint(m)%field_2d)) & + CS%interior_tendency_out_zint(m)%field_2d(i,j) = (sum(dz(:) * & + MARBL_instances%interior_tendencies(m,:))) + + if (allocated(CS%interior_tendency_out_zint_100m(m)%field_2d)) then + CS%interior_tendency_out_zint_100m(m)%field_2d(i,j) = 0. + do k=1,GV%ke + if (zi(k) < US%m_to_Z * 100.) then + CS%interior_tendency_out_zint_100m(m)%field_2d(i,j) = & + CS%interior_tendency_out_zint_100m(m)%field_2d(i,j) + GV%H_to_m * dz(k) * & + MARBL_instances%interior_tendencies(m,k) + elseif (zi(k-1) < US%m_to_Z * 100.) then + CS%interior_tendency_out_zint_100m(m)%field_2d(i,j) = & + CS%interior_tendency_out_zint_100m(m)%field_2d(i,j) + GV%H_to_m * dz(k) * & + ((US%m_to_Z * 100. - zi(k-1)) / (zi(k) - zi(k-1))) * & + MARBL_instances%interior_tendencies(m,k) + else + exit + endif + enddo + endif + enddo + + ! * Interior tendency output + do m=1,CS%ito_cnt + CS%ITO(i,j,:,m) = & + MARBL_instances%interior_tendency_output%outputs_for_GCM(m)%forcing_field_1d(1,:) + enddo + + enddo + enddo + + if (CS%debug) then + do m=1,CS%ntr + call hchksum(CS%tracer_data(m)%tr(:,:,m), & + trim(MARBL_instances%tracer_metadata(m)%short_name)//' post source-sink', G%HI) + enddo + endif + + ! (5) Post diagnostics from our buffer + ! i. Interior tendency diagnostics (mix of 2D and 3D) + ! ii. Interior tendencies themselves + ! iii. Forcing fields + if (CS%bot_flux_to_tend_id > 0) & + call post_data(CS%bot_flux_to_tend_id, bot_flux_to_tend(:, :, :), CS%diag) + + do m=1,size(CS%interior_tendency_diags) + if (CS%interior_tendency_diags(m)%id > 0) then + if (allocated(CS%interior_tendency_diags(m)%field_2d)) then + if (real(MARBL_instances%interior_tendency_diags%diags(m)%ref_depth) == 0.) then + call post_data(CS%interior_tendency_diags(m)%id, & + CS%interior_tendency_diags(m)%field_2d(:,:), CS%diag) + else ! non-zero ref-depth + ref_mask(:, :) = 0. + do j=js,je ; do i=is,ie + if (G%bathyT(i,j) > real(MARBL_instances%interior_tendency_diags%diags(m)%ref_depth)) & + ref_mask(i,j) = 1. + enddo ; enddo + call post_data(CS%interior_tendency_diags(m)%id, & + CS%interior_tendency_diags(m)%field_2d(:,:), CS%diag, mask=ref_mask(:,:)) + endif + elseif (allocated(CS%interior_tendency_diags(m)%field_3d)) then + call post_data(CS%interior_tendency_diags(m)%id, & + CS%interior_tendency_diags(m)%field_3d(:,:,:), CS%diag) + else + write(log_message, "(A, I0, A, I0, A)") "Diagnostic number ", m, " post id ", & + CS%interior_tendency_diags(m)%id," did not allocate 2D or 3D array" + call MOM_error(FATAL, log_message) + endif + endif + enddo + + do m=1,CS%ntr + if (allocated(CS%interior_tendency_out(m)%field_3d)) & + call post_data(CS%interior_tendency_out(m)%id, & + CS%interior_tendency_out(m)%field_3d(:,:,:), CS%diag) + if (allocated(CS%interior_tendency_out_zint(m)%field_2d)) & + call post_data(CS%interior_tendency_out_zint(m)%id, & + CS%interior_tendency_out_zint(m)%field_2d(:,:), CS%diag) + if (allocated(CS%interior_tendency_out_zint_100m(m)%field_2d)) & + call post_data(CS%interior_tendency_out_zint_100m(m)%id, & + CS%interior_tendency_out_zint_100m(m)%field_2d(:,:), CS%diag) + enddo + + if (CS%ice_ncat > 0) then + do m=1,CS%ice_ncat+1 + if (CS%fracr_cat_id(m) > 0) & + call post_data(CS%fracr_cat_id(m), fluxes%fracr_cat(:,:,m), CS%diag) + if (CS%qsw_cat_id(m) > 0) & + call post_data(CS%qsw_cat_id(m), fluxes%qsw_cat(:,:,m), CS%diag) + enddo + endif + + +end subroutine MARBL_tracers_column_physics + +!> This subroutine reads time-varying forcing from files +subroutine MARBL_tracers_set_forcing(day_start, G, CS) + + type(time_type), intent(in) :: day_start !< Start time of the fluxes. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(MARBL_tracers_CS), pointer :: CS !< The control structure returned by a + + ! Fraction of river nutrients in refractory pools + real, parameter :: DONriv_refract = 0.1 + real, parameter :: DOCriv_refract = 0.2 + real, parameter :: DOPriv_refract = 0.025 + + real, dimension(SZI_(G),SZJ_(G)) :: riv_flux_in !< The field read in from forcing file with time dimension + type(time_type) :: Time_forcing !< For reading river flux fields, we use a modified version of Time + integer :: i, j, k, is, ie, js, je, m + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + ! Abiotic DIC forcing + if (CS%abio_dic_on) then + ! Read d14c bands + do m=1,3 + Time_forcing = map_model_time_to_forcing_time(day_start, CS%d14c_dataset(m)) + call time_interp_external(CS%id_d14c(m),Time_forcing,CS%d14c_bands(m)) + enddo + + ! Set d14c according to the bands + do j=js,je ; do i=is,ie + if (G%geoLatT(i,j) > 30.) then + CS%d14c(i,j) = CS%d14c_bands(1) + elseif (G%geoLatT(i,j) > -30.) then + CS%d14c(i,j) = CS%d14c_bands(2) + else + CS%d14c(i,j) = CS%d14c_bands(3) + endif + enddo ; enddo + endif + + ! River fluxes + if (CS%read_riv_fluxes) then + CS%RIV_FLUXES(:,:,:) = 0. + Time_forcing = map_model_time_to_forcing_time(day_start, CS%riv_flux_dataset) + + ! DIN river flux affects NO3, ALK, and ALK_ALT_CO2 + call time_interp_external(CS%id_din_riv,Time_forcing,riv_flux_in) + + if (CS%tracer_inds%no3_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%no3_ind) = G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%alk_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_ind) = CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_ind) - & + G%mask2dT(i,j) *riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%alk_alt_co2_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_alt_co2_ind) = & + CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_alt_co2_ind) - G%mask2dT(i,j) *riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_dip_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%po4_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%po4_ind) = G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_don_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%don_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%don_ind) = G%mask2dT(i,j) * (1. - DONriv_refract) * & + riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%donr_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%donr_ind) = G%mask2dT(i,j) * DONriv_refract * & + riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_dop_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%dop_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%dop_ind) = G%mask2dT(i,j) * (1. - DOPriv_refract) * & + riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%dopr_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%dopr_ind) = G%mask2dT(i,j) * DOPriv_refract * & + riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_dsi_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%sio3_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%sio3_ind) = G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_dfe_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%fe_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%fe_ind) = G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_dic_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%dic_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%dic_ind) = G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%dic_alt_co2_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%dic_alt_co2_ind) = G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_alk_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%alk_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_ind) = CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_ind) + & + G%mask2dT(i,j) *riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%alk_alt_co2_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_alt_co2_ind) = & + CS%RIV_FLUXES(i,j,CS%tracer_inds%alk_alt_co2_ind) + G%mask2dT(i,j) * riv_flux_in(i,j) + enddo ; enddo + endif + + call time_interp_external(CS%id_doc_riv,Time_forcing,riv_flux_in) + if (CS%tracer_inds%doc_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%doc_ind) = G%mask2dT(i,j) * (1. - DOCriv_refract) * & + riv_flux_in(i,j) + enddo ; enddo + endif + if (CS%tracer_inds%docr_ind > 0) then + do j=js,je ; do i=is,ie + CS%RIV_FLUXES(i,j,CS%tracer_inds%docr_ind) = G%mask2dT(i,j) * DOCriv_refract * & + riv_flux_in(i,j) + enddo ; enddo + endif + endif + + ! Tracer restoring + do m=1,CS%restore_count + call time_interp_external(CS%id_tracer_restoring(m),day_start,CS%restoring_in(:,:,:,m)) + do k=1,CS%restoring_nz ; do j=js,je ; do i=is,ie + CS%restoring_in(i,j,k,m) = G%mask2dT(i,j) * CS%restoring_in(i,j,k,m) + enddo ; enddo ; enddo + enddo + + ! Post Forcing to Diagnostics + if (CS%read_riv_fluxes) then + if (CS%no3_riv_flux > 0 .and. CS%tracer_inds%no3_ind > 0) & + call post_data(CS%no3_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%no3_ind), CS%diag) + if (CS%po4_riv_flux > 0 .and. CS%tracer_inds%po4_ind > 0) & + call post_data(CS%po4_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%po4_ind), CS%diag) + if (CS%don_riv_flux > 0 .and. CS%tracer_inds%don_ind > 0) & + call post_data(CS%don_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%don_ind), CS%diag) + if (CS%donr_riv_flux > 0 .and. CS%tracer_inds%donr_ind > 0) & + call post_data(CS%donr_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%donr_ind), CS%diag) + if (CS%dop_riv_flux > 0 .and. CS%tracer_inds%dop_ind > 0) & + call post_data(CS%dop_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%dop_ind), CS%diag) + if (CS%dopr_riv_flux > 0 .and. CS%tracer_inds%dopr_ind > 0) & + call post_data(CS%dopr_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%dopr_ind), CS%diag) + if (CS%sio3_riv_flux > 0 .and. CS%tracer_inds%sio3_ind > 0) & + call post_data(CS%sio3_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%sio3_ind), CS%diag) + if (CS%fe_riv_flux > 0 .and. CS%tracer_inds%fe_ind > 0) & + call post_data(CS%fe_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%fe_ind), CS%diag) + if (CS%doc_riv_flux > 0 .and. CS%tracer_inds%doc_ind > 0) & + call post_data(CS%doc_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%doc_ind), CS%diag) + if (CS%docr_riv_flux > 0 .and. CS%tracer_inds%docr_ind > 0) & + call post_data(CS%docr_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%docr_ind), CS%diag) + if (CS%alk_riv_flux > 0 .and. CS%tracer_inds%alk_ind > 0) & + call post_data(CS%alk_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%alk_ind), CS%diag) + if (CS%alk_alt_co2_riv_flux > 0 .and. CS%tracer_inds%alk_alt_co2_ind > 0) & + call post_data(CS%alk_alt_co2_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%alk_alt_co2_ind), & + CS%diag) + if (CS%dic_riv_flux > 0 .and. CS%tracer_inds%dic_ind > 0) & + call post_data(CS%dic_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%dic_ind), CS%diag) + if (CS%dic_alt_co2_riv_flux > 0 .and. CS%tracer_inds%dic_alt_co2_ind > 0) & + call post_data(CS%dic_alt_co2_riv_flux, CS%RIV_FLUXES(:,:,CS%tracer_inds%dic_alt_co2_ind), & + CS%diag) + endif + if (CS%abio_dic_on) then + if (CS%d14c_id > 0) & + call post_data(CS%d14c_id, CS%d14c, CS%diag) + endif + +end subroutine MARBL_tracers_set_forcing + +!> This function calculates the mass-weighted integral of all tracer stocks, +!! returning the number of stocks it has calculated. If the stock_index +!! is present, only the stock corresponding to that coded index is returned. +function MARBL_tracers_stock(h, stocks, G, GV, CS, names, units, stock_index) + real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(EFP_type), dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of + !! each tracer, in kg times concentration units + !! [kg conc]. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(MARBL_tracers_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_MARBL_tracers. + character(len=*), dimension(:), intent(out) :: names !< the names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< the units of the stocks calculated. + integer, optional, intent(in) :: stock_index !< the coded index of a specific stock + !! being sought. + integer :: MARBL_tracers_stock !< Return value: the number of stocks + !! calculated here. + +! Local variables + integer :: i, j, k, is, ie, js, je, nz, m + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + MARBL_tracers_stock = 0 + if (.not.associated(CS)) return + if (CS%ntr < 1) return + + if (present(stock_index)) then ; if (stock_index > 0) then + ! Check whether this stock is available from this routine. + + ! No stocks from this routine are being checked yet. Return 0. + return + endif ; endif + + do m=1,CS%ntr + call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="MARBL_tracers_stock") + units(m) = trim(units(m))//" kg" + stocks(m) = global_mass_int_EFP(h, G, GV, CS%tracer_data(m)%tr(:,:,:), on_PE_only=.true.) + enddo + MARBL_tracers_stock = CS%ntr + +end function MARBL_tracers_stock + +!> This subroutine extracts the surface fields from this tracer package that +!! are to be shared with the atmosphere in coupled configurations. +subroutine MARBL_tracers_surface_state(sfc_state, G, US, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(surface), intent(inout) :: sfc_state !< A structure containing fields that + !! describe the surface state of the ocean. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(MARBL_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_MARBL_tracers. + + ! Local variables + integer :: i, j, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + if (.not.associated(CS)) return + + if (allocated(sfc_state%fco2)) then + do j=js,je ; do i=is,ie + ! 44e-6 converts mmol/m^2/s (positive down) to kg CO2/m^2/s (positive down) + sfc_state%fco2(i,j) = US%kg_m2s_to_RZ_T * (44.0e-6 * CS%SFO(i,j,CS%flux_co2_ind)) + enddo ; enddo + endif + +end subroutine MARBL_tracers_surface_state + +!> Copy the requested interior tendency output field into an array. +subroutine MARBL_tracers_get(name, G, GV, array, CS) + + character(len=*), intent(in) :: name !< Name of requested tracer. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: array !< Array filled by this routine. + type(MARBL_tracers_CS), pointer :: CS !< Pointer to the control structure for this module. + + character(len=128), parameter :: sub_name = 'MARBL_tracers_get' + character(len=128) :: log_message + + array(:,:,:) = 0.0 + select case(trim(name)) + case ('Chl') + array(:,:,:) = CS%ITO(:,:,:,CS%total_Chl_ind) + case DEFAULT + write(log_message, "(3A)") "'", trim(name), & + "' is not a valid interior tendency output field name" + call MOM_error(FATAL, log_message) + end select + +end subroutine MARBL_tracers_get + +!> Clean up any allocated memory after the run. +subroutine MARBL_tracers_end(CS) + type(MARBL_tracers_CS), pointer, intent(inout) :: CS !< The control structure returned by a previous + !! call to register_MARBL_tracers. + + integer :: m + + call print_marbl_log(MARBL_instances%StatusLog) + call MARBL_instances%StatusLog%erase() + call MARBL_instances%shutdown() + ! TODO: print MARBL timers to stdout as well + + if (associated(CS)) then + if (allocated(CS%tracer_data)) then + do m=1,CS%ntr + if (associated(CS%tracer_data(m)%tr)) deallocate(CS%tracer_data(m)%tr) + enddo + deallocate(CS%tracer_data) + endif + if (allocated(CS%ind_tr)) deallocate(CS%ind_tr) + if (allocated(CS%id_surface_flux_out)) deallocate(CS%id_surface_flux_out) + if (allocated(CS%interior_tendency_out)) deallocate(CS%interior_tendency_out) + if (allocated(CS%interior_tendency_out_zint)) deallocate(CS%interior_tendency_out_zint) + if (allocated(CS%interior_tendency_out_zint_100m)) & + deallocate(CS%interior_tendency_out_zint_100m) + if (allocated(CS%fracr_cat_id)) deallocate(CS%fracr_cat_id) + if (allocated(CS%qsw_cat_id)) deallocate(CS%qsw_cat_id) + if (allocated(CS%STF)) deallocate(CS%STF) + if (allocated(CS%RIV_FLUXES)) deallocate(CS%RIV_FLUXES) + if (allocated(CS%SFO)) deallocate(CS%SFO) + if (allocated(CS%tracer_restoring_ind)) deallocate(CS%tracer_restoring_ind) + if (allocated(CS%tracer_I_tau_ind)) deallocate(CS%tracer_I_tau_ind) + if (allocated(CS%fesedflux_in)) deallocate(CS%fesedflux_in) + if (allocated(CS%feventflux_in)) deallocate(CS%feventflux_in) + if (allocated(CS%I_tau)) deallocate(CS%I_tau) + deallocate(CS) + endif +end subroutine MARBL_tracers_end + +subroutine set_riv_flux_tracer_inds(CS) + + type(MARBL_tracers_CS), pointer, intent(inout) :: CS !< The MARBL tracers control structure + + character(len=256) :: log_message + character(len=48) :: name ! A variable's name in a NetCDF file. + integer :: m + + ! Initialize tracers from file (unless they were initialized by restart file) + ! Also save indices of tracers that have river fluxes + CS%tracer_inds%no3_ind = 0 + CS%tracer_inds%po4_ind = 0 + CS%tracer_inds%don_ind = 0 + CS%tracer_inds%donr_ind = 0 + CS%tracer_inds%dop_ind = 0 + CS%tracer_inds%dopr_ind = 0 + CS%tracer_inds%sio3_ind = 0 + CS%tracer_inds%fe_ind = 0 + CS%tracer_inds%doc_ind = 0 + CS%tracer_inds%docr_ind = 0 + CS%tracer_inds%alk_ind = 0 + CS%tracer_inds%alk_alt_co2_ind = 0 + CS%tracer_inds%dic_ind = 0 + CS%tracer_inds%dic_alt_co2_ind = 0 + do m=1,CS%ntr + name = MARBL_instances%tracer_metadata(m)%short_name + if (trim(name) == "NO3") then + CS%tracer_inds%no3_ind = m + elseif (trim(name) == "PO4") then + CS%tracer_inds%po4_ind = m + elseif (trim(name) == "DON") then + CS%tracer_inds%don_ind = m + elseif (trim(name) == "DONr") then + CS%tracer_inds%donr_ind = m + elseif (trim(name) == "DOP") then + CS%tracer_inds%dop_ind = m + elseif (trim(name) == "DOPr") then + CS%tracer_inds%dopr_ind = m + elseif (trim(name) == "SiO3") then + CS%tracer_inds%sio3_ind = m + elseif (trim(name) == "Fe") then + CS%tracer_inds%fe_ind = m + elseif (trim(name) == "DOC") then + CS%tracer_inds%doc_ind = m + elseif (trim(name) == "DOCr") then + CS%tracer_inds%docr_ind = m + elseif (trim(name) == "ALK") then + CS%tracer_inds%alk_ind = m + elseif (trim(name) == "ALK_ALT_CO2") then + CS%tracer_inds%alk_alt_co2_ind = m + elseif (trim(name) == "DIC") then + CS%tracer_inds%dic_ind = m + elseif (trim(name) == "DIC_ALT_CO2") then + CS%tracer_inds%dic_alt_co2_ind = m + endif + enddo + + ! Log indices for each tracer to ensure we set them all correctly + write(log_message, "(A,I0)") "NO3 index: ", CS%tracer_inds%no3_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "PO4 index: ", CS%tracer_inds%po4_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DON index: ", CS%tracer_inds%don_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DONr index: ", CS%tracer_inds%donr_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DOP index: ", CS%tracer_inds%dop_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DOPr index: ", CS%tracer_inds%dopr_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "SiO3 index: ", CS%tracer_inds%sio3_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "Fe index: ", CS%tracer_inds%fe_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DOC index: ", CS%tracer_inds%doc_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DOCr index: ", CS%tracer_inds%docr_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "ALK index: ", CS%tracer_inds%alk_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "ALK_ALT_CO2 index: ", CS%tracer_inds%alk_alt_co2_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DIC index: ", CS%tracer_inds%dic_ind + call MOM_error(NOTE, log_message) + write(log_message, "(A,I0)") "DIC_ALT_CO2 index: ", CS%tracer_inds%dic_alt_co2_ind + call MOM_error(NOTE, log_message) + +end subroutine set_riv_flux_tracer_inds + +! TODO: some log messages come from a specific grid point, and this routine +! needs to include the location in the preamble +!> This subroutine writes the contents of the MARBL log using MOM_error(NOTE, ...). +subroutine print_marbl_log(log_to_print, G, i, j) + + use marbl_logging, only : marbl_status_log_entry_type + use marbl_logging, only : marbl_log_type + use MOM_coms, only : PE_here + + class(marbl_log_type), intent(in) :: log_to_print !< MARBL log to include in MOM6 logfile + type(ocean_grid_type), optional, intent(in) :: G !< The ocean's grid structure + integer, optional, intent(in) :: i !< i of (i,j) index of column providing the log + integer, optional, intent(in) :: j !< j of (i,j) index of column providing the log + + character(len=*), parameter :: subname = 'MARBL_tracers:print_marbl_log' + character(len=256) :: message_prefix, message_location, log_message + type(marbl_status_log_entry_type), pointer :: tmp + integer :: msg_lev, elem_old + + ! elem_old is used to keep track of whether all messages are coming from the same point + elem_old = -1 + write(message_prefix, "(A,I0,A)") '(Task ', PE_here(), ')' + + tmp => log_to_print%FullLog + do while (associated(tmp)) + ! 1) Do I need to write this message? Yes, if all tasks should write this + ! or if I am master_task + if ((.not. tmp%lonly_master_writes) .or. is_root_PE()) then + ! 2) Print message location? (only if ElementInd changed and is positive; requires G) + if ((present(G)) .and. (tmp%ElementInd .ne. elem_old)) then + if (tmp%ElementInd .gt. 0) then + if (present(i) .and. present(j)) then + write(message_location, "(A,F8.3,A,F7.3,A,I0,A,I0,A,I0)") & + 'Message from (lon, lat) (', G%geoLonT(i,j), ', ', G%geoLatT(i,j), & + '), which is global (i,j) (', i + G%HI%idg_offset, ', ', j + G%HI%jdg_offset, & + '). Level: ', tmp%ElementInd + else + write(message_location, "(A)") "Grid cell responsible for message is unknown" + endif ! i,j present + ! master task does not need prefix + if (is_root_PE()) then + write(log_message, "(A)") trim(message_location) + msg_lev = NOTE + else + write(log_message, "(A,1X,A)") trim(message_prefix), trim(message_location) + msg_lev = WARNING + endif ! print message prefix? + call MOM_error(msg_lev, log_message, all_print=.true.) + endif ! ElementInd > 0 + elem_old = tmp%ElementInd + endif ! ElementInd /= elem_old + + ! 3) Write message from the log + ! master task does not need prefix + if (is_root_PE()) then + write(log_message, "(A)") trim(tmp%LogMessage) + msg_lev = NOTE + else + write(log_message, "(A,1X,A)") trim(message_prefix), trim(tmp%LogMessage) + msg_lev = WARNING + endif ! print message prefix? + call MOM_error(msg_lev, log_message, all_print=.true.) + endif ! write the message? + tmp => tmp%next + enddo + + if (log_to_print%labort_marbl) then + call MOM_error(WARNING, 'ERROR reported from MARBL library', all_print=.true.) + call MOM_error(FATAL, 'Stopping in ' // subname) + endif + +end subroutine print_marbl_log + +!> \namespace MARBL_tracers +!! +!! This module contains the code that is needed to provide +!! the MARBL BGC tracer library with necessary forcings and +!! apply the resulting surface fluxes and tendencies to the +!! requested tracers. + +end module MARBL_tracers diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index 5b9af238d6..a81a42b428 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -421,7 +421,7 @@ integer function find_minimum(x, s, e) if (x(i) < minimum) then ! if x(i) less than the min? minimum = x(i) ! Yes, a new minimum found location = i ! record its position - end if + endif enddo find_minimum = location ! return the position end function find_minimum diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index fab7da3917..caa2d10a04 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -16,7 +16,7 @@ module MOM_tracer_Z_init #include -public tracer_Z_init, tracer_Z_init_array, determine_temperature +public tracer_Z_init, read_Z_edges, tracer_Z_init_array, determine_temperature ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index c8ce2f5f75..ef80f9d23c 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -38,6 +38,10 @@ module MOM_tracer_flow_control use ideal_age_example, only : register_ideal_age_tracer, initialize_ideal_age_tracer use ideal_age_example, only : ideal_age_tracer_column_physics, ideal_age_tracer_surface_state use ideal_age_example, only : ideal_age_stock, ideal_age_example_end, ideal_age_tracer_CS +use MARBL_tracers, only : register_MARBL_tracers, initialize_MARBL_tracers +use MARBL_tracers, only : MARBL_tracers_column_physics, MARBL_tracers_set_forcing +use MARBL_tracers, only : MARBL_tracers_surface_state, MARBL_tracers_get +use MARBL_tracers, only : MARBL_tracers_stock, MARBL_tracers_end, MARBL_tracers_CS use regional_dyes, only : register_dye_tracer, initialize_dye_tracer use regional_dyes, only : dye_tracer_column_physics, dye_tracer_surface_state use regional_dyes, only : dye_stock, regional_dyes_end, dye_tracer_CS @@ -85,6 +89,7 @@ module MOM_tracer_flow_control logical :: use_ISOMIP_tracer = .false. !< If true, use the ISOMPE_tracer package logical :: use_RGC_tracer =.false. !< If true, use the RGC_tracer package logical :: use_ideal_age = .false. !< If true, use the ideal age tracer package + logical :: use_MARBL_tracers = .false. !< If true, use the MARBL tracer package logical :: use_regional_dyes = .false. !< If true, use the regional dyes tracer package logical :: use_oil = .false. !< If true, use the oil tracer package logical :: use_advection_test_tracer = .false. !< If true, use the advection_test_tracer package @@ -95,12 +100,14 @@ module MOM_tracer_flow_control logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package logical :: use_dyed_obc_tracer = .false. !< If true, use the dyed OBC tracer package logical :: use_nw2_tracers = .false. !< If true, use the NW2 tracer package + logical :: get_chl_from_MARBL = .false. !< If true, use the MARBL-provided Chl for shortwave penetration !>@{ Pointers to the control strucures for the tracer packages type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL() type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL() type(ISOMIP_tracer_CS), pointer :: ISOMIP_tracer_CSp => NULL() type(RGC_tracer_CS), pointer :: RGC_tracer_CSp => NULL() type(ideal_age_tracer_CS), pointer :: ideal_age_tracer_CSp => NULL() + type(MARBL_tracers_CS), pointer :: MARBL_tracers_CSp => NULL() type(dye_tracer_CS), pointer :: dye_tracer_CSp => NULL() type(oil_tracer_CS), pointer :: oil_tracer_CSp => NULL() type(advection_test_tracer_CS), pointer :: advection_test_tracer_CSp => NULL() @@ -193,6 +200,9 @@ subroutine call_tracer_register(G, GV, US, param_file, CS, tr_Reg, restart_CS) call get_param(param_file, mdl, "USE_IDEAL_AGE_TRACER", CS%use_ideal_age, & "If true, use the ideal_age_example tracer package.", & default=.false.) + call get_param(param_file, mdl, "USE_MARBL_TRACERS", CS%use_marbl_tracers, & + "If true, use the MARBL tracer package.", & + default=.false.) call get_param(param_file, mdl, "USE_REGIONAL_DYES", CS%use_regional_dyes, & "If true, use the regional_dyes tracer package.", & default=.false.) @@ -243,6 +253,9 @@ subroutine call_tracer_register(G, GV, US, param_file, CS, tr_Reg, restart_CS) if (CS%use_ideal_age) CS%use_ideal_age = & register_ideal_age_tracer(G%HI, GV, param_file, CS%ideal_age_tracer_CSp, & tr_Reg, restart_CS) + if (CS%use_MARBL_tracers) CS%use_MARBL_tracers = & + register_MARBL_tracers(G%HI, GV, US, param_file, CS%MARBL_tracers_CSp, & + tr_Reg, restart_CS, CS%get_chl_from_MARBL) if (CS%use_regional_dyes) CS%use_regional_dyes = & register_dye_tracer(G%HI, GV, US, param_file, CS%dye_tracer_CSp, & tr_Reg, restart_CS) @@ -327,6 +340,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag if (CS%use_ideal_age) & call initialize_ideal_age_tracer(restart, day, G, GV, US, h, diag, OBC, CS%ideal_age_tracer_CSp, & sponge_CSp) + if (CS%use_MARBL_tracers) & + call initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag, OBC, CS%MARBL_tracers_CSp, & + sponge_CSp) if (CS%use_regional_dyes) & call initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS%dye_tracer_CSp, sponge_CSp, tv) if (CS%use_oil) & @@ -386,7 +402,9 @@ subroutine get_chl_from_model(Chl_array, G, GV, CS) type(tracer_flow_control_CS), pointer :: CS !< The control structure returned by a !! previous call to call_tracer_register. - if (CS%use_MOM_generic_tracer) then + if (CS%get_chl_from_MARBL) then + call MARBL_tracers_get('Chl', G, GV, Chl_array, CS%MARBL_tracers_CSp) + elseif (CS%use_MOM_generic_tracer) then call MOM_generic_tracer_get('chl', 'field', Chl_array, CS%MOM_generic_tracer_CSp) else call MOM_error(FATAL, "get_chl_from_model was called in a configuration "// & @@ -424,6 +442,9 @@ subroutine call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G call CFC_cap_set_forcing(sfc_state, fluxes, day_start, day_interval, G, US, Rho0, & CS%CFC_cap_CSp) + if (CS%use_MARBL_tracers) & + call MARBL_tracers_set_forcing(day_start, G, CS%MARBL_tracers_CSp) + end subroutine call_tracer_set_forcing !> This subroutine calls all registered tracer column physics subroutines. @@ -494,6 +515,13 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, evap_CFL_limit=evap_CFL_limit, & minimum_forcing_depth=minimum_forcing_depth, & Hbl=Hml) + if (CS%use_MARBL_tracers) & + call MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%MARBL_tracers_CSp, tv, & + KPP_CSp=KPP_CSp, & + nonLocalTrans=nonLocalTrans, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) if (CS%use_regional_dyes) & call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, tv, CS%dye_tracer_CSp, & @@ -570,6 +598,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, if (CS%use_ideal_age) & call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, tv, CS%ideal_age_tracer_CSp, Hbl=Hml) + if (CS%use_MARBL_tracers) & + call MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%MARBL_tracers_CSp, tv, & + KPP_CSp=KPP_CSp, & + nonLocalTrans=nonLocalTrans) if (CS%use_regional_dyes) & call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, tv, CS%dye_tracer_CSp) @@ -691,6 +724,12 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, US, CS, stock_names, stock call store_stocks("ideal_age_example", ns, names, units, values_EFP, index, stock_val_EFP, & set_pkg_name, max_ns, ns_tot, stock_names, stock_units) endif + if (CS%use_MARBL_tracers) then + ns = MARBL_tracers_stock(h, values_EFP, G, GV, CS%MARBL_tracers_CSp, & + names, units, stock_index) + call store_stocks("MARBL_tracers", ns, names, units, values_EFP, index, stock_val_EFP, & + set_pkg_name, max_ns, ns_tot, stock_names, stock_units) + endif if (CS%use_regional_dyes) then ns = dye_stock(h, values_EFP, G, GV, CS%dye_tracer_CSp, names, units, stock_index) call store_stocks("regional_dyes", ns, names, units, values_EFP, index, stock_val_EFP, & @@ -844,6 +883,8 @@ subroutine call_tracer_surface_state(sfc_state, h, G, GV, US, CS) call ISOMIP_tracer_surface_state(sfc_state, h, G, GV, CS%ISOMIP_tracer_CSp) if (CS%use_ideal_age) & call ideal_age_tracer_surface_state(sfc_state, h, G, GV, CS%ideal_age_tracer_CSp) + if (CS%use_MARBL_tracers) & + call MARBL_tracers_surface_state(sfc_state, G, US, CS%MARBL_tracers_CSp) if (CS%use_regional_dyes) & call dye_tracer_surface_state(sfc_state, h, G, GV, CS%dye_tracer_CSp) if (CS%use_oil) & @@ -867,6 +908,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_ISOMIP_tracer) call ISOMIP_tracer_end(CS%ISOMIP_tracer_CSp) if (CS%use_RGC_tracer) call RGC_tracer_end(CS%RGC_tracer_CSp) if (CS%use_ideal_age) call ideal_age_example_end(CS%ideal_age_tracer_CSp) + if (CS%use_MARBL_tracers) call MARBL_tracers_end(CS%MARBL_tracers_CSp) if (CS%use_regional_dyes) call regional_dyes_end(CS%dye_tracer_CSp) if (CS%use_oil) call oil_tracer_end(CS%oil_tracer_CSp) if (CS%use_advection_test_tracer) call advection_test_tracer_end(CS%advection_test_tracer_CSp) diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index c01419f3f8..c7d11b6030 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -390,6 +390,16 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') endif + Tr%id_zint = register_diag_field("ocean_model", trim(shortnm)//"_zint", & + diag%axesT1, Time, & + "Thickness-weighted integral of " // trim(longname), & + trim(units) // " m") + Tr%id_zint_100m = register_diag_field("ocean_model", trim(shortnm)//"_zint_100m", & + diag%axesT1, Time, & + "Thickness-weighted integral of "// trim(longname) // " over top 100m", & + trim(units) // " m") + Tr%id_surf = register_diag_field("ocean_model", trim(shortnm)//"_SURF", & + diag%axesT1, Time, "Surface values of "// trim(longname), trim(units)) if (Tr%id_adx > 0) call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_ady > 0) call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz) if (Tr%id_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz) @@ -592,7 +602,7 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u conversion = GV%H_to_kg_m2 else conversion = Tr%conv_scale - end if + endif ! We actually want conversion=Tr%conv_scale for all tracers, but introducing the local variable ! 'conversion' and setting it to GV%H_to_kg_m2 instead of 0.001*GV%H_to_kg_m2 for salt tracers ! keeps changes introduced by this refactoring limited to round-off level; as it turns out, @@ -716,12 +726,42 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) intent(in) :: h_diag !< Layer thicknesses on which to post fields [H ~> m or kg m-2] type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output - integer :: i, j, k, is, ie, js, je, nz, m - real :: work2d(SZI_(G),SZJ_(G)) + integer :: i, j, k, is, ie, js, je, nz, m, khi + real :: frac_under_100m(SZI_(G),SZJ_(G),SZK_(GV)) + real :: work2d(SZI_(G),SZJ_(G)), ztop(SZI_(G),SZJ_(G)), zbot(SZI_(G),SZJ_(G)) type(tracer_type), pointer :: Tr=>NULL() is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + ! If any tracers are posting 100m vertical integrals, compute weights + frac_under_100m(:,:,:) = 0.0 + ! khi will be the largest layer index corresponding where ztop < 100m and ztop >= 100m + ! in any column (we can reduce computation of 100m integrals by only looping through khi + ! rather than GV%ke) + khi = 0 + do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then + Tr => Reg%Tr(m) + if (Tr%id_zint_100m > 0) then + zbot(:,:) = 0.0 + do k=1, nz + do j=js,je ; do i=is,ie + ztop(i,j) = zbot(i,j) + zbot(i,j) = ztop(i,j) + h_diag(i,j,k)*GV%H_to_m + if (zbot(i,j) <= 100.0) then + frac_under_100m(i,j,k) = 1.0 + elseif (ztop(i,j) < 100.0) then + frac_under_100m(i,j,k) = (100.0 - ztop(i,j)) / (zbot(i,j) - ztop(i,j)) + else + frac_under_100m(i,j,k) = 0.0 + endif + ! frac_under_100m(i,j,k) = max(0, min(1.0, (100.0 - ztop(i,j)) / (zbot(i,j) - ztop(i,j)))) + enddo ; enddo + if (any(frac_under_100m(:,:,k) > 0)) khi = k + enddo + exit + endif + endif; enddo + do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then Tr => Reg%Tr(m) if (Tr%id_tr_post_horzn> 0) call post_data(Tr%id_tr_post_horzn, Tr%t, diag) @@ -741,6 +781,28 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) enddo ; enddo ; enddo call post_data(Tr%id_adv_xy_2d, work2d, diag) endif + + ! A few diagnostics introduce with MARBL driver + ! Compute full-depth vertical integral + if (Tr%id_zint > 0) then + work2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do i=is,ie + work2d(i,j) = work2d(i,j) + (h_diag(i,j,k)*GV%H_to_m)*tr%t(i,j,k) + enddo ; enddo ; enddo + call post_data(Tr%id_zint, work2d, diag) + endif + + ! Compute 100m vertical integral + if (Tr%id_zint_100m > 0) then + work2d(:,:) = 0.0 + do k=1,khi ; do j=js,je ; do i=is,ie + work2d(i,j) = work2d(i,j) + frac_under_100m(i,j,k)*((h_diag(i,j,k)*GV%H_to_m)*tr%t(i,j,k)) + enddo ; enddo ; enddo + call post_data(Tr%id_zint_100m, work2d, diag) + endif + + ! Surface values of tracers + if (Tr%id_SURF > 0) call post_data(Tr%id_SURF, Tr%t(:,:,1), diag) endif ; enddo end subroutine post_tracer_transport_diagnostics diff --git a/src/tracer/MOM_tracer_types.F90 b/src/tracer/MOM_tracer_types.F90 index bdae8bcee9..55326a0b1b 100644 --- a/src/tracer/MOM_tracer_types.F90 +++ b/src/tracer/MOM_tracer_types.F90 @@ -111,6 +111,7 @@ module MOM_tracer_types integer :: id_remap_conc = -1, id_remap_cont = -1, id_remap_cont_2d = -1 integer :: id_tendency = -1, id_trxh_tendency = -1, id_trxh_tendency_2d = -1 integer :: id_tr_vardec = -1 + integer :: id_zint = -1, id_zint_100m = -1, id_surf = -1 integer :: id_net_surfflux = -1, id_NLT_tendency = -1, id_NLT_budget = -1 !>@} end type tracer_type