diff --git a/dglc/cime_config/config_component.xml b/dglc/cime_config/config_component.xml index 15271a5b8..a47b9d5e9 100644 --- a/dglc/cime_config/config_component.xml +++ b/dglc/cime_config/config_component.xml @@ -43,7 +43,7 @@ - + logical FALSE run_component_dglc @@ -58,7 +58,7 @@ - + logical FALSE run_component_dglc diff --git a/dglc/cime_config/testdefs/testlist_dglc.xml b/dglc/cime_config/testdefs/testlist_dglc.xml index 194773fa7..8383469c7 100644 --- a/dglc/cime_config/testdefs/testlist_dglc.xml +++ b/dglc/cime_config/testdefs/testlist_dglc.xml @@ -1,7 +1,7 @@ - + @@ -10,7 +10,7 @@ - + @@ -19,7 +19,7 @@ - + diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 9a5431e62..e4392105e 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -30,7 +30,8 @@ module dglc_datamode_noevolve_mod ! Data structure to enable multiple ice sheets type icesheet_ptr_t - real(r8), pointer :: ptr(:) => null() ! pointer to array + real(r8), pointer :: ptr(:) => null() ! pointer to array + real(r8), pointer :: ptr2d(:,:) => null() ! pointer to 2d array endtype icesheet_ptr_t ! Export fields @@ -39,10 +40,14 @@ module dglc_datamode_noevolve_mod type(icesheet_ptr_t), allocatable :: Sg_ice_covered(:) type(icesheet_ptr_t), allocatable :: Sg_icemask(:) type(icesheet_ptr_t), allocatable :: Sg_icemask_coupled_fluxes(:) + type(icesheet_ptr_t), allocatable :: Fogg_rofi(:) ! Import fields + integer, parameter :: nlev_import = 30 type(icesheet_ptr_t), allocatable :: Sl_tsrf(:) type(icesheet_ptr_t), allocatable :: Flgl_qice(:) +! type(icesheet_ptr_t), allocatable :: So_t(:) +! type(icesheet_ptr_t), allocatable :: So_q(:) ! Export Field names character(len=*), parameter :: field_out_area = 'Sg_area' @@ -50,10 +55,13 @@ module dglc_datamode_noevolve_mod character(len=*), parameter :: field_out_ice_covered = 'Sg_ice_covered' character(len=*), parameter :: field_out_icemask = 'Sg_icemask' character(len=*), parameter :: field_out_icemask_coupled_fluxes = 'Sg_icemask_coupled_fluxes' + character(len=*), parameter :: field_out_rofi = 'Fogg_rofi' ! Import Field names character(len=*), parameter :: field_in_tsrf = 'Sl_tsrf' character(len=*), parameter :: field_in_qice = 'Flgl_qice' + character(len=*), parameter :: field_in_so_t_depth = 'So_t_depth' + character(len=*), parameter :: field_in_so_s_depth = 'So_s_depth' character(*) , parameter :: rpfile = 'rpointer.glc' character(*) , parameter :: u_FILE_u = & @@ -96,6 +104,7 @@ subroutine dglc_datamode_noevolve_advertise(NStateExp, fldsexport, NStateImp, fl call dshr_fldList_add(fldsExport, field_out_topo) call dshr_fldList_add(fldsExport, field_out_icemask) call dshr_fldList_add(fldsExport, field_out_icemask_coupled_fluxes) + call dshr_fldList_add(fldsExport, field_out_rofi) do ns = 1,num_icesheets write(cnum,'(i0)') ns @@ -112,6 +121,9 @@ subroutine dglc_datamode_noevolve_advertise(NStateExp, fldsexport, NStateImp, fl call dshr_fldList_add(fldsImport, trim(flds_scalar_name)) call dshr_fldList_add(fldsImport, field_in_tsrf) call dshr_fldList_add(fldsImport, field_in_qice) + ! TODO: Add namelist for this + call dshr_fldList_add(fldsImport, field_in_so_t_depth, ungridded_lbound=1, ungridded_ubound=nlev_import) + call dshr_fldList_add(fldsImport, field_in_so_s_depth, ungridded_lbound=1, ungridded_ubound=nlev_import) do ns = 1,num_icesheets write(cnum,'(i0)') ns @@ -150,6 +162,7 @@ subroutine dglc_datamode_noevolve_init_pointers(NStateExp, NstateImp, rc) allocate(Sg_ice_covered(num_icesheets)) allocate(Sg_icemask(num_icesheets)) allocate(Sg_icemask_coupled_fluxes(num_icesheets)) + allocate(Fogg_rofi(num_icesheets)) do ns = 1,num_icesheets call dshr_state_getfldptr(NStateExp(ns), field_out_area, fldptr1=Sg_area(ns)%ptr, rc=rc) @@ -162,6 +175,8 @@ subroutine dglc_datamode_noevolve_init_pointers(NStateExp, NstateImp, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call dshr_state_getfldptr(NStateExp(ns), field_out_icemask_coupled_fluxes, fldptr1=Sg_icemask_coupled_fluxes(ns)%ptr, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateExp(ns), field_out_rofi, fldptr1=Fogg_rofi(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return end do ! initialize pointers to import fields if appropriate @@ -169,16 +184,15 @@ subroutine dglc_datamode_noevolve_init_pointers(NStateExp, NstateImp, rc) allocate(Flgl_qice(num_icesheets)) do ns = 1,num_icesheets - ! NOTE: the field is connected ONLY if the MED->GLC entry is in the nuopc.runconfig file - ! This restriction occurs even if the field was advertised - if (NUOPC_IsConnected(NStateImp(ns), fieldName=field_in_tsrf)) then - call dshr_state_getfldptr(NStateImp(ns), field_in_tsrf, fldptr1=Sl_tsrf(ns)%ptr, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - end if - if (NUOPC_IsConnected(NStateImp(ns), fieldName=field_in_qice)) then - call dshr_state_getfldptr(NStateImp(ns), field_in_qice, fldptr1=Flgl_qice(ns)%ptr, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (.not. NUOPC_IsConnected(NStateImp(ns), fieldName=field_in_tsrf)) then + ! NOTE: the field is connected ONLY if the MED->GLC entry is in the nuopc.runconfig file + ! This restriction occurs even if the field was advertised + call shr_sys_abort(trim(subname)//": MED->GLC must appear in run sequence") end if + call dshr_state_getfldptr(NStateImp(ns), field_in_tsrf, fldptr1=Sl_tsrf(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(NStateImp(ns), field_in_qice, fldptr1=Flgl_qice(ns)%ptr, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return end do end subroutine dglc_datamode_noevolve_init_pointers @@ -226,124 +240,140 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & rc = ESMF_SUCCESS - if (initialized_noevolve) then - RETURN - end if + if (.not. initialized_noevolve) then - do ns = 1,num_icesheets + do ns = 1,num_icesheets - ! Get mesh info - call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - allocate(gindex(lsize)) - call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=gindex, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - ! Determine "glc_area" ; - ! Sg_areas is in radians - ! SHR_CONST_REARTH is the radius of earth in m - ! model_internal_gridsize is the internal model gridsize in m - do ng = 1,lsize - Sg_area(ns)%ptr(ng) = (model_internal_gridsize(ns)/SHR_CONST_REARTH)**2 - end do + ! Get mesh info + call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + allocate(gindex(lsize)) + call ESMF_DistGridGet(distGrid, localDe=0, seqIndexList=gindex, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Create module level field bundle - fldbun_noevolve = ESMF_FieldBundleCreate(rc=rc) ! input field bundle - - ! "ice thickness" ; - field_noevolve = ESMF_FieldCreate(model_meshes(ns), ESMF_TYPEKIND_R8, name='thk', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_FieldBundleAdd(fldbun_noevolve, (/field_noevolve/), rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - - ! "bed topography" ; - field_noevolve = ESMF_FieldCreate(model_meshes(ns), ESMF_TYPEKIND_R8, name='topg', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_FieldBundleAdd(fldbun_noevolve, (/field_noevolve/), rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - - ! Create pioid and pio_iodesc at the module level - inquire(file=trim(model_datafiles(ns)), exist=exists) - if (.not.exists) then - write(6,'(a)')' ERROR: model input file '//trim(model_datafiles(ns))//' does not exist' - call shr_sys_abort() - end if - rcode = pio_openfile(pio_subsystem, pioid, io_type, trim(model_datafiles(ns)), pio_nowrite) - call pio_seterrorhandling(pioid, PIO_BCAST_ERROR) - rcode = pio_inq_varid(pioid, 'thk', varid) - rcode = pio_inq_varndims(pioid, varid, ndims) - allocate(dimid(ndims)) - rcode = pio_inq_vardimid(pioid, varid, dimid(1:ndims)) - rcode = pio_inq_dimlen(pioid, dimid(1), nxg) - rcode = pio_inq_dimlen(pioid, dimid(2), nyg) - call pio_initdecomp(pio_subsystem, pio_double, (/nxg,nyg/), gindex, pio_iodesc) - deallocate(dimid) - deallocate(gindex) - - ! Read in the data into the appropriate field bundle pointers - ! Note that Sg_ice_covered(ns)%ptr points into the data for - ! the Sg_ice_covered field in NStateExp(ns) - ! Note that Sg_topo(ns)%ptr points into the data for - ! the Sg_topon NStateExp(ns) - ! Note that topog is bedrock topography - - call dshr_fldbun_getFldPtr(fldbun_noevolve, 'topg', topog, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - rcode = pio_inq_varid(pioid, 'topg', varid) - call pio_read_darray(pioid, varid, pio_iodesc, topog, rcode) - - call dshr_fldbun_getFldPtr(fldbun_noevolve, 'thk', thck, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - rcode = pio_inq_varid(pioid, 'thk', varid) - call pio_read_darray(pioid, varid, pio_iodesc, thck, rcode) - - allocate(lsrf(lsize)) - allocate(usrf(lsize)) - - rhoi = SHR_CONST_RHOICE ! 0.917e3 - rhoo = SHR_CONST_RHOSW ! 1.026e3 - eus = 0 - do ng = 1,lsize - if (topog(ng) - eus < (-rhoi/rhoo) * thck(ng)) then - lsrf(ng) = (-rhoi/rhoo) * thck(ng) - else - lsrf(ng) = topog(ng) - end if - usrf(ng) = max(0.d0, thck(ng) + lsrf(ng)) - - if (is_in_active_grid(usrf(ng))) then - Sg_icemask(ns)%ptr(ng) = 1.d0 - Sg_icemask_coupled_fluxes(ns)%ptr(ng) = 0.d0 - if (is_ice_covered(thck(ng))) then - Sg_ice_covered(ns)%ptr(ng) = 1.d0 - else - Sg_ice_covered(ns)%ptr(ng) = 0.d0 - end if - ! Note that we use the same method for computing topo whether this point is - ! ice-covered or ice-free. This is in contrast to the method for computing - ! ice-free topo in glint_upscaling_gcm. - Sg_topo(ns)%ptr(ng) = thk0 * usrf(ng) - else - ! Note that this logic implies that if (in theory) we had an ice-covered - ! point outside the "active grid", it will get classified as ice-free for - ! these purposes. This mimics the logic currently in glint_upscaling_gcm. - Sg_icemask(ns)%ptr(ng) = 0.d0 - Sg_icemask_coupled_fluxes(ns)%ptr(ng) = 0.d0 - Sg_ice_covered(ns)%ptr(ng) = 0.d0 - Sg_topo(ns)%ptr(ng) = 0.d0 - end if - end do + ! Determine "glc_area" ; + ! Sg_areas is in radians + ! SHR_CONST_REARTH is the radius of earth in m + ! model_internal_gridsize is the internal model gridsize in m + do ng = 1,lsize + Sg_area(ns)%ptr(ng) = (model_internal_gridsize(ns)/SHR_CONST_REARTH)**2 + end do - deallocate(lsrf) - deallocate(usrf) + ! Create module level field bundle + fldbun_noevolve = ESMF_FieldBundleCreate(rc=rc) ! input field bundle - call pio_closefile(pioid) - call pio_freedecomp(pio_subsystem, pio_iodesc) + ! "ice thickness" ; + field_noevolve = ESMF_FieldCreate(model_meshes(ns), ESMF_TYPEKIND_R8, name='thk', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldBundleAdd(fldbun_noevolve, (/field_noevolve/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! "bed topography" ; + field_noevolve = ESMF_FieldCreate(model_meshes(ns), ESMF_TYPEKIND_R8, name='topg', meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_FieldBundleAdd(fldbun_noevolve, (/field_noevolve/), rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Create pioid and pio_iodesc at the module level + inquire(file=trim(model_datafiles(ns)), exist=exists) + if (.not.exists) then + write(6,'(a)')' ERROR: model input file '//trim(model_datafiles(ns))//' does not exist' + call shr_sys_abort() + end if + rcode = pio_openfile(pio_subsystem, pioid, io_type, trim(model_datafiles(ns)), pio_nowrite) + call pio_seterrorhandling(pioid, PIO_BCAST_ERROR) + rcode = pio_inq_varid(pioid, 'thk', varid) + rcode = pio_inq_varndims(pioid, varid, ndims) + allocate(dimid(ndims)) + rcode = pio_inq_vardimid(pioid, varid, dimid(1:ndims)) + rcode = pio_inq_dimlen(pioid, dimid(1), nxg) + rcode = pio_inq_dimlen(pioid, dimid(2), nyg) + call pio_initdecomp(pio_subsystem, pio_double, (/nxg,nyg/), gindex, pio_iodesc) + deallocate(dimid) + deallocate(gindex) + + ! Read in the data into the appropriate field bundle pointers + ! Note that Sg_ice_covered(ns)%ptr points into the data for + ! the Sg_ice_covered field in NStateExp(ns) + ! Note that Sg_topo(ns)%ptr points into the data for + ! the Sg_topon NStateExp(ns) + ! Note that topog is bedrock topography + + call dshr_fldbun_getFldPtr(fldbun_noevolve, 'topg', topog, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + rcode = pio_inq_varid(pioid, 'topg', varid) + call pio_read_darray(pioid, varid, pio_iodesc, topog, rcode) - end do ! end loop over ice sheets + call dshr_fldbun_getFldPtr(fldbun_noevolve, 'thk', thck, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + rcode = pio_inq_varid(pioid, 'thk', varid) + call pio_read_darray(pioid, varid, pio_iodesc, thck, rcode) + + allocate(lsrf(lsize)) + allocate(usrf(lsize)) + + rhoi = SHR_CONST_RHOICE ! 0.917e3 + rhoo = SHR_CONST_RHOSW ! 1.026e3 + eus = 0 + do ng = 1,lsize + if (topog(ng) - eus < (-rhoi/rhoo) * thck(ng)) then + lsrf(ng) = (-rhoi/rhoo) * thck(ng) + else + lsrf(ng) = topog(ng) + end if + usrf(ng) = max(0.d0, thck(ng) + lsrf(ng)) + + ! The export field 'ice_mask_coupled_fluxes' determines who is handling the + ! runoff associated with the surface mass balance + ! If its 0 -then ctsm needs to handle it. + ! Since we want dglc to handle it no evolve mode - then + ! ice_mask_coupled_fluxes to be identical to the mask + + if (is_in_active_grid(usrf(ng))) then + Sg_icemask(ns)%ptr(ng) = 1.d0 + Sg_icemask_coupled_fluxes(ns)%ptr(ng) = 1.d0 + if (is_ice_covered(thck(ng))) then + Sg_ice_covered(ns)%ptr(ng) = 1.d0 + else + Sg_ice_covered(ns)%ptr(ng) = 0.d0 + end if + ! Note that we use the same method for computing topo whether this point is + ! ice-covered or ice-free. This is in contrast to the method for computing + ! ice-free topo in glint_upscaling_gcm. + Sg_topo(ns)%ptr(ng) = thk0 * usrf(ng) + else + ! Note that this logic implies that if (in theory) we had an ice-covered + ! point outside the "active grid", it will get classified as ice-free for + ! these purposes. This mimics the logic currently in glint_upscaling_gcm. + Sg_icemask(ns)%ptr(ng) = 0.d0 + Sg_icemask_coupled_fluxes(ns)%ptr(ng) = 0.d0 + Sg_ice_covered(ns)%ptr(ng) = 0.d0 + Sg_topo(ns)%ptr(ng) = 0.d0 + end if + end do + + deallocate(lsrf) + deallocate(usrf) + + call pio_closefile(pioid) + call pio_freedecomp(pio_subsystem, pio_iodesc) + + end do ! end loop over ice sheets + + end if + + if (initialized_noevolve) then + ! Compute Fogg_rofi + do ns = 1,num_icesheets + do ng = 1,size(Fogg_rofi(ns)%ptr) + Fogg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + end do + end do + end if + ! Set initialized flag initialized_noevolve = .true. end subroutine dglc_datamode_noevolve_advance diff --git a/dglc/glc_comp_nuopc.F90 b/dglc/glc_comp_nuopc.F90 index 233794fd3..09c858ff1 100644 --- a/dglc/glc_comp_nuopc.F90 +++ b/dglc/glc_comp_nuopc.F90 @@ -9,11 +9,13 @@ module cdeps_dglc_comp !---------------------------------------------------------------------------- use ESMF , only : ESMF_VM, ESMF_VmGet, ESMF_VMBroadcast, ESMF_GridCompGet use ESMF , only : ESMF_Mesh, ESMF_GridComp, ESMF_State, ESMF_Clock, ESMF_Time - use ESMF , only : ESMF_SUCCESS, ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_METHOD_INITIALIZE - use ESMF , only : ESMF_TraceRegionEnter, ESMF_TraceRegionExit, ESMF_ClockGet - use ESMF , only : ESMF_TimeGet, ESMF_TimeInterval, ESMF_Field, ESMF_MAXSTR - use ESMF , only : ESMF_Alarm, ESMF_MethodRemove, ESMF_MethodAdd - use ESMF , only : ESMF_GridCompSetEntryPoint, ESMF_ClockGetAlarm, ESMF_AlarmIsRinging + use ESMF , only : ESMF_SUCCESS, ESMF_FAILURE, ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_METHOD_INITIALIZE + use ESMF , only : ESMF_TraceRegionEnter, ESMF_TraceRegionExit + use ESMF , only : ESMF_TimeGet, ESMF_TimeInterval, ESMF_TimeIntervalGet, ESMF_Field, ESMF_MAXSTR + use ESMF , only : ESMF_MethodRemove, ESMF_MethodAdd + use ESMF , only : ESMF_GridCompSetEntryPoint + use ESMF , only : ESMF_Alarm, ESMF_AlarmSet, ESMF_AlarmIsRinging, ESMF_ALARMLIST_ALL + use ESMF , only : ESMF_ClockGet, ESMF_ClockSet, ESMF_ClockAdvance, ESMF_ClockGetAlarm, ESMF_ClockGetAlarmList use ESMF , only : ESMF_StateGet, operator(+), ESMF_AlarmRingerOff, ESMF_LogWrite use ESMF , only : ESMF_Field, ESMF_FieldGet, ESMF_VmLogMemInfo use ESMF , only : ESMF_MeshCreate, ESMF_FILEFORMAT_ESMFMESH @@ -35,7 +37,7 @@ module cdeps_dglc_comp #endif use dshr_methods_mod , only : dshr_state_diagnose, chkerr, memcheck use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_advance, shr_strdata_init_from_config - use dshr_mod , only : dshr_model_initphase, dshr_init, dshr_mesh_init + use dshr_mod , only : dshr_model_initphase, dshr_init, dshr_mesh_init, dshr_alarm_init use dshr_mod , only : dshr_state_setscalar, dshr_set_runclock, dshr_check_restart_alarm use dshr_dfield_mod , only : dfield_type, dshr_dfield_add, dshr_dfield_copy use dshr_fldlist_mod , only : fldlist_type, dshr_fldlist_realize @@ -52,6 +54,7 @@ module cdeps_dglc_comp public :: SetVM private :: InitializeAdvertise private :: InitializeRealize + private :: ModelSetRunClock private :: ModelAdvance private :: dglc_comp_run private :: ModelFinalize @@ -157,7 +160,9 @@ subroutine SetServices(gcomp, rc) call ESMF_MethodRemove(gcomp, label=model_label_SetRunClock, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call NUOPC_CompSpecialize(gcomp, specLabel=model_label_SetRunClock, specRoutine=dshr_set_runclock, rc=rc) + ! The following specialization does not use dshr_set_runclock since we also need to add an alarm + ! to determine if the model inputs are valid + call NUOPC_CompSpecialize(gcomp, specLabel=model_label_SetRunClock, specRoutine=ModelSetRunClock, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call NUOPC_CompSpecialize(gcomp, specLabel=model_label_Finalize, specRoutine=ModelFinalize, rc=rc) @@ -425,7 +430,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) end do ! end loop over ice sheets ! Run dglc - call dglc_comp_run(clock, current_ymd, current_tod, restart_write=.false., rc=rc) + call dglc_comp_run(clock, current_ymd, current_tod, restart_write=.false., valid_inputs=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_TraceRegionExit('dglc_strdata_init') @@ -450,10 +455,12 @@ subroutine ModelAdvance(gcomp, rc) integer :: mon ! month integer :: day ! day in month logical :: restart_write + type(ESMF_Alarm) :: valid_alarm + logical :: valid_inputs ! if true, inputs from mediator are valid + character(len=CS) :: timestring character(len=*),parameter :: subname=trim(module_name)//':(ModelAdvance) ' !------------------------------------------------------------------------------- - rc = ESMF_SUCCESS call shr_log_setLogUnit(logunit) @@ -463,7 +470,17 @@ subroutine ModelAdvance(gcomp, rc) call NUOPC_ModelGet(gcomp, modelClock=clock, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! For nuopc - the component clock is advanced at the end of the time interval + ! Determine if inputs from mediator are valid + call ESMF_ClockGetAlarm(clock, alarmname='alarm_valid_inputs', alarm=valid_alarm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (ESMF_AlarmIsRinging(valid_alarm, rc=rc)) then + valid_inputs = .true. + call ESMF_AlarmRingerOff( valid_alarm, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + valid_inputs = .false. + endif + ! Need to advance nuopc one timestep ahead for shr_strdata time interpolation call ESMF_ClockGet( clock, currTime=currTime, timeStep=timeStep, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -471,18 +488,23 @@ subroutine ModelAdvance(gcomp, rc) call ESMF_TimeGet( nextTime, yy=yr, mm=mon, dd=day, s=next_tod, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return call shr_cal_ymd2date(yr, mon, day, next_ymd) + if (my_task == main_task) then + call ESMF_TimeGet(currTime, timestring=timestring, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + write(logunit,'(a,l)') trim(timestring)//': valid_input for dglc is ',valid_inputs + end if ! determine if will write restart restart_write = dshr_check_restart_alarm(clock, rc=rc) ! run dglc - call dglc_comp_run(clock, next_ymd, next_tod, restart_write, rc=rc) + call dglc_comp_run(clock, next_ymd, next_tod, restart_write, valid_inputs, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end subroutine ModelAdvance !=============================================================================== - subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, rc) + subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, valid_inputs, rc) ! -------------------------- ! advance dglc @@ -493,6 +515,7 @@ subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, rc) integer , intent(in) :: target_ymd ! model date integer , intent(in) :: target_tod ! model sec into model date logical , intent(in) :: restart_write + logical , intent(in) :: valid_inputs integer , intent(out) :: rc ! local variables @@ -565,9 +588,11 @@ subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, rc) ! Perform data mode specific calculations select case (trim(datamode)) case('noevolve') - call dglc_datamode_noevolve_advance(sdat(1)%pio_subsystem, sdat(1)%io_type, sdat(1)%io_format, & - model_meshes, model_internal_gridsize, model_datafiles, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (valid_inputs) then + call dglc_datamode_noevolve_advance(sdat(1)%pio_subsystem, sdat(1)%io_type, sdat(1)%io_format, & + model_meshes, model_internal_gridsize, model_datafiles, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if end select ! Write restarts if needed @@ -632,6 +657,150 @@ end subroutine dglc_init_dfields end subroutine dglc_comp_run + !=============================================================================== + + subroutine ModelSetRunClock(gcomp, rc) + + ! input/output variables + type(ESMF_GridComp) :: gcomp + integer, intent(out) :: rc + + ! local variables + type(ESMF_Clock) :: mclock, dclock + type(ESMF_Time) :: mcurrtime, dcurrtime + type(ESMF_Time) :: mstoptime + type(ESMF_TimeInterval) :: mtimestep, dtimestep + character(len=256) :: cvalue + character(len=256) :: restart_option ! Restart option units + integer :: restart_n ! Number until restart interval + integer :: restart_ymd ! Restart date (YYYYMMDD) + type(ESMF_ALARM) :: restart_alarm + character(len=256) :: stop_option ! Stop option units + integer :: stop_n ! Number until stop interval + integer :: stop_ymd ! Stop date (YYYYMMDD) + type(ESMF_ALARM) :: stop_alarm + integer :: alarmcount + character(len=CS) :: glc_avg_period ! averaging period in mediator + type(ESMF_ALARM) :: valid_alarm ! model alarm + integer :: dtime + character(len=*),parameter :: subname='glc_comp_nuopc:(dglc_set_runclock) ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! query the Component for its clocks + call NUOPC_ModelGet(gcomp, driverClock=dclock, modelClock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_ClockGet(dclock, currTime=dcurrtime, timeStep=dtimestep, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_ClockGet(mclock, currTime=mcurrtime, timeStep=mtimestep, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! force model clock currtime and timestep to match driver and set stoptime + mstoptime = mcurrtime + dtimestep + call ESMF_ClockSet(mclock, currTime=dcurrtime, timeStep=dtimestep, stopTime=mstoptime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! determine number of alarms + call ESMF_ClockGetAlarmList(mclock, alarmlistflag=ESMF_ALARMLIST_ALL, alarmCount=alarmCount, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (alarmCount == 0) then + + !---------------- + ! glc valid input alarm + !---------------- + call NUOPC_CompAttributeGet(gcomp, name="glc_avg_period", value=glc_avg_period, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (trim(glc_avg_period) == 'hour') then + call dshr_alarm_init(mclock, valid_alarm, 'nhours', opt_n=1, alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(glc_avg_period) == 'day') then + call dshr_alarm_init(mclock, valid_alarm, 'ndays' , opt_n=1, alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(glc_avg_period) == 'yearly') then + call dshr_alarm_init(mclock, valid_alarm, 'yearly', alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(glc_avg_period) == 'glc_coupling_period') then + call ESMF_TimeIntervalGet(mtimestep, s=dtime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call dshr_alarm_init(mclock, valid_alarm, 'nseconds', opt_n=dtime, alarmname='alarm_valid_inputs', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + call ESMF_LogWrite(trim(subname)// ": ERROR glc_avg_period = "//trim(glc_avg_period)//" not supported", & + ESMF_LOGMSG_INFO, rc=rc) + rc = ESMF_FAILURE + RETURN + end if + + call ESMF_AlarmSet(valid_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + !---------------- + ! Restart alarm + !---------------- + call ESMF_LogWrite(subname//'setting restart alarm for dglc' , ESMF_LOGMSG_INFO) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompAttributeGet(gcomp, name="restart_option", value=restart_option, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompAttributeGet(gcomp, name="restart_n", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) restart_n + + call NUOPC_CompAttributeGet(gcomp, name="restart_ymd", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) restart_ymd + + call dshr_alarm_init(mclock, restart_alarm, restart_option, & + opt_n = restart_n, & + opt_ymd = restart_ymd, & + RefTime = mcurrTime, & + alarmname = 'alarm_restart', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_AlarmSet(restart_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + !---------------- + ! Stop alarm + !---------------- + call ESMF_LogWrite(subname//'setting stop alarm for dglc' , ESMF_LOGMSG_INFO) + call NUOPC_CompAttributeGet(gcomp, name="stop_option", value=stop_option, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call NUOPC_CompAttributeGet(gcomp, name="stop_n", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) stop_n + + call NUOPC_CompAttributeGet(gcomp, name="stop_ymd", value=cvalue, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) stop_ymd + + call dshr_alarm_init(mclock, stop_alarm, stop_option, & + opt_n = stop_n, & + opt_ymd = stop_ymd, & + RefTime = mcurrTime, & + alarmname = 'alarm_stop', rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_AlarmSet(stop_alarm, clock=mclock, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + end if + + ! Advance model clock to trigger alarms then reset model clock back to currtime + call ESMF_ClockAdvance(mclock,rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_ClockSet(mclock, currTime=dcurrtime, timeStep=dtimestep, stopTime=mstoptime, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + end subroutine ModelSetRunClock + !=============================================================================== subroutine ModelFinalize(gcomp, rc) type(ESMF_GridComp) :: gcomp diff --git a/docn/CMakeLists.txt b/docn/CMakeLists.txt index 007d595ca..256dc304d 100644 --- a/docn/CMakeLists.txt +++ b/docn/CMakeLists.txt @@ -5,6 +5,7 @@ set(SRCFILES ocn_comp_nuopc.F90 docn_datamode_aquaplanet_mod.F90 docn_datamode_iaf_mod.F90 docn_datamode_cplhist_mod.F90 + docn_datamode_multilev_mod.F90 docn_import_data_mod.F90) foreach(FILE ${SRCFILES}) diff --git a/docn/cime_config/config_component.xml b/docn/cime_config/config_component.xml index 3a98b7878..ab88c6fed 100644 --- a/docn/cime_config/config_component.xml +++ b/docn/cime_config/config_component.xml @@ -13,7 +13,7 @@ This file may have ocn desc entries. --> - DOCN + DOCN prescribed ocean mode slab ocean mode aquaplanet slab ocean mode @@ -32,6 +32,7 @@ file input aquaplanet sst globally constant SST for idealized experiments, such as RCE mediator history output for ocean fields imported to mediator + input stream files have multi level data @@ -45,7 +46,7 @@ char - prescribed,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,som,som_aquap,sst_aquap_constant,interannual,cplhist + prescribed,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,som,som_aquap,sst_aquap_constant,interannual,cplhist,multilev prescribed prescribed @@ -65,6 +66,7 @@ sst_aquapfile sst_aquap_constant cplhist + multilev run_component_docn env_run.xml diff --git a/docn/cime_config/namelist_definition_docn.xml b/docn/cime_config/namelist_definition_docn.xml index a44ae344d..b2f788113 100644 --- a/docn/cime_config/namelist_definition_docn.xml +++ b/docn/cime_config/namelist_definition_docn.xml @@ -32,6 +32,7 @@ '' '' '' + sst_salinity_depth_blom @@ -39,7 +40,7 @@ char docn docn_nml - sstdata,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,sst_aquap_constant,som,som_aquap,iaf,cplhist + sstdata,sst_aquap1,sst_aquap2,sst_aquap3,sst_aquap4,sst_aquap5,sst_aquap6,sst_aquap7,sst_aquap8,sst_aquap9,sst_aquap10,sst_aquapfile,sst_aquap_constant,som,som_aquap,iaf,cplhist,multilev General method that operates on the data for a given docn_mode. ==> dataMode = "sstdata" @@ -107,6 +108,7 @@ sst_aquap_file sst_aquap_constant cplhist + multilev diff --git a/docn/cime_config/stream_definition_docn.xml b/docn/cime_config/stream_definition_docn.xml index 0dd4baf54..4e4dd6655 100644 --- a/docn/cime_config/stream_definition_docn.xml +++ b/docn/cime_config/stream_definition_docn.xml @@ -206,6 +206,104 @@ 1.e30 - single + single + + + + + $DIN_LOC_ROOT/share/meshes/gx1v7_151008_ESMFmesh.nc + + + /glade/collections/cmip/CMIP6/OMIP/NCAR/CESM2/omip2/r1i1p1f1/Omon/thetao/gn/v20190802/thetao_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc + + + thetao So_t_depth + + lev + + bilinear + + null + 1 + 245 + 245 + 0 + + linear + + + cycle + + + 1.5 + + single + + + + + $DIN_LOC_ROOT/share/meshes/gx1v7_151008_ESMFmesh.nc + + + /glade/collections/cmip/CMIP6/OMIP/NCAR/CESM2/omip2/r1i1p1f1/Omon/so/gn/v20190802/so_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc + + + so So_s_depth + + lev + + bilinear + + null + 1 + 245 + 245 + 0 + + linear + + + cycle + + + 1.5 + + single + + + + + $DIN_LOC_ROOT/share/meshes/tnx1v4_20170601_cdf5_ESMFmesh.nc + + + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2300.nc + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2301.nc + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2302.nc + $DIN_LOC_ROOT/ocn/docn7/MULTILEV/N1850frc2G_f09_tn14_gl4_SMB1_celo.blom.hm.2303.nc + + + templvl So_t_depth + salnlvl So_s_depth + + depth + + bilinear + + null + 1 + 2300 + 2303 + 0 + + linear + + + cycle + + + 1.5 + + single + diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 new file mode 100644 index 000000000..897797f3f --- /dev/null +++ b/docn/docn_datamode_multilev_mod.F90 @@ -0,0 +1,229 @@ +module docn_datamode_multilev_mod + use ESMF , only : ESMF_State, ESMF_LOGMSG_INFO, ESMF_LogWrite, ESMF_SUCCESS + use NUOPC , only : NUOPC_Advertise + use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs + use shr_const_mod , only : shr_const_TkFrz, shr_const_pi, shr_const_ocn_ref_sal, shr_const_spval + use shr_sys_mod , only : shr_sys_abort + use dshr_methods_mod , only : dshr_state_getfldptr, dshr_fldbun_getfldptr, chkerr + use dshr_fldlist_mod , only : fldlist_type, dshr_fldlist_add + use dshr_mod , only : dshr_restart_read, dshr_restart_write + use dshr_strdata_mod , only : shr_strdata_get_stream_pointer, shr_strdata_type + + implicit none + private ! except + + public :: docn_datamode_multilev_advertise + public :: docn_datamode_multilev_init_pointers + public :: docn_datamode_multilev_advance + public :: docn_datamode_multilev_restart_read + public :: docn_datamode_multilev_restart_write + + ! pointers to export fields + real(r8), pointer :: So_omask(:) => null() ! real ocean fraction sent to mediator + real(r8), pointer :: So_t_depth(:,:) => null() + real(r8), pointer :: So_s_depth(:,:) => null() + + ! pointers to stream fields + real(r8), pointer :: stream_So_t_depth(:,:) => null() + real(r8), pointer :: stream_So_s_depth(:,:) => null() + + integer, parameter :: nlev_export = 30 + real(r8) :: vertical_levels(nlev_export) = (/ & + 30., 90., 150., 210., 270., 330., 390., 450., 510., 570., & + 630., 690., 750., 810., 870., 930., 990., 1050., 1110., 1170., & + 1230., 1290., 1350., 1410., 1470., 1530., 1590., 1650., 1710., 1770. /) + + ! constants + character(*) , parameter :: nullstr = 'null' + character(*) , parameter :: rpfile = 'rpointer.ocn' + character(*) , parameter :: u_FILE_u = & + __FILE__ + +!=============================================================================== +contains +!=============================================================================== + + subroutine docn_datamode_multilev_advertise(exportState, fldsexport, flds_scalar_name, rc) + + ! input/output variables + type(esmf_State) , intent(inout) :: exportState + type(fldlist_type) , pointer :: fldsexport + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(out) :: rc + + ! local variables + type(fldlist_type), pointer :: fldList + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! Advertise export fields + call dshr_fldList_add(fldsExport, trim(flds_scalar_name)) + call dshr_fldList_add(fldsExport, 'So_omask') + call dshr_fldList_add(fldsExport, 'So_t_depth', ungridded_lbound=1, ungridded_ubound=nlev_export) + call dshr_fldList_add(fldsExport, 'So_s_depth', ungridded_lbound=1, ungridded_ubound=nlev_export) + + fldlist => fldsExport ! the head of the linked list + do while (associated(fldlist)) + call NUOPC_Advertise(exportState, standardName=fldlist%stdname, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_LogWrite('(docn_comp_advertise): Fr_ocn'//trim(fldList%stdname), ESMF_LOGMSG_INFO) + fldList => fldList%next + enddo + + end subroutine docn_datamode_multilev_advertise + + !=============================================================================== + subroutine docn_datamode_multilev_init_pointers(exportState, sdat, ocn_fraction, rc) + + ! input/output variables + type(ESMF_State) , intent(inout) :: exportState + type(shr_strdata_type) , intent(in) :: sdat + real(r8) , intent(in) :: ocn_fraction(:) + integer , intent(out) :: rc + + ! local variables + character(len=*), parameter :: subname='(docn_init_pointers): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! initialize pointers to stream fields + ! this has the full set of leveles in the stream data + call shr_strdata_get_stream_pointer( sdat, 'So_t_depth', stream_So_t_depth, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call shr_strdata_get_stream_pointer( sdat, 'So_s_depth', stream_So_s_depth, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! initialize pointers to export fields + ! the export state has only nlev_export levels + call dshr_state_getfldptr(exportState, 'So_omask' , fldptr1=So_omask , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_t_depth' , fldptr2=So_t_depth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call dshr_state_getfldptr(exportState, 'So_s_depth' , fldptr2=So_s_depth , rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! Initialize export state pointers to non-zero + So_t_depth(:,:) = shr_const_TkFrz + So_s_depth(:,:) = shr_const_ocn_ref_sal + + ! Set export state ocean fraction (So_omask) + So_omask(:) = ocn_fraction(:) + + end subroutine docn_datamode_multilev_init_pointers + + !=============================================================================== + subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) + + ! input/output variables + type(shr_strdata_type) , intent(in) :: sdat + integer , intent(in) :: logunit + logical , intent(in) :: mainproc + integer , intent(out) :: rc + + ! local variables + integer :: i,ki,ko + integer :: nlev_stream + integer :: stream_index + logical :: level_found + real(r8) :: factor + real(r8), allocatable :: stream_vlevs(:) + logical :: first_time = .true. + character(len=*), parameter :: subname='(docn_datamode_multilev): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + ! For now assume that all the streams have the same vertical levels + stream_index = 1 + + nlev_stream = sdat%pstrm(stream_index)%stream_nlev + allocate(stream_vlevs(nlev_stream)) + stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) + + do ko = 1,nlev_export + level_found = .false. + do ki = 1,nlev_stream-1 + if (vertical_levels(ko) > stream_vlevs(ki) .and. vertical_levels(ko) <= stream_vlevs(ki+1)) then + if (mainproc .and. first_time) then + write(logunit,'(a,3(i5,2x),3(f13.5,2x))') & + 'vertical interpolation: ki,ko,ki+1,lev(ki),lev(ko),lev(ki+1) = ',& + ki,ko,ki+1,stream_vlevs(ki), vertical_levels(ko), stream_vlevs(ki+1) + end if + level_found = .true. + do i = 1,size(So_omask) + if (So_omask(i) == 0.) then + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval + else + ! Assume input T forcing is in degrees C + if (stream_So_t_depth(ki+1,i) > 1.e10) then + if (stream_So_t_depth(ki,i) > 1.e10) then + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval + else + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + shr_const_tkfrz + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + end if + else + factor = (stream_So_t_depth(ki+1,i)-stream_So_t_depth(ki,i))/(stream_vlevs(ki+1)-stream_vlevs(ki)) + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + (vertical_levels(ko)-stream_vlevs(ki))*factor + So_t_depth(ko,i) = So_t_depth(ko,i) + shr_const_tkfrz + + factor = (stream_So_s_depth(ki+1,i)-stream_So_s_depth(ki,i))/(stream_vlevs(ki+1)-stream_vlevs(ki)) + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + (vertical_levels(ko)-stream_vlevs(ki))*factor + end if + end if + end do + end if + end do + if (.not. level_found) then + call shr_sys_abort("ERROR: could not find level bounds for vertical interpolation") + end if + end do + + first_time = .false. + + end subroutine docn_datamode_multilev_advance + + !=============================================================================== + subroutine docn_datamode_multilev_restart_write(case_name, inst_suffix, ymd, tod, & + logunit, my_task, sdat) + + ! write restart file + + ! input/output variables + character(len=*) , intent(in) :: case_name + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: ymd ! model date + integer , intent(in) :: tod ! model sec into model date + integer , intent(in) :: logunit + integer , intent(in) :: my_task + type(shr_strdata_type) , intent(inout) :: sdat + !------------------------------------------------------------------------------- + + call dshr_restart_write(rpfile, case_name, 'docn', inst_suffix, ymd, tod, & + logunit, my_task, sdat) + + end subroutine docn_datamode_multilev_restart_write + + !=============================================================================== + subroutine docn_datamode_multilev_restart_read(rest_filem, inst_suffix, logunit, my_task, mpicom, sdat) + + ! read restart file + + ! input/output arguments + character(len=*) , intent(inout) :: rest_filem + character(len=*) , intent(in) :: inst_suffix + integer , intent(in) :: logunit + integer , intent(in) :: my_task + integer , intent(in) :: mpicom + type(shr_strdata_type) , intent(inout) :: sdat + !------------------------------------------------------------------------------- + + call dshr_restart_read(rest_filem, rpfile, inst_suffix, nullstr, logunit, my_task, mpicom, sdat) + + end subroutine docn_datamode_multilev_restart_read + +end module docn_datamode_multilev_mod diff --git a/docn/ocn_comp_nuopc.F90 b/docn/ocn_comp_nuopc.F90 index 4bdeb81f6..ced8cb86d 100644 --- a/docn/ocn_comp_nuopc.F90 +++ b/docn/ocn_comp_nuopc.F90 @@ -26,7 +26,7 @@ module cdeps_docn_comp use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs use shr_sys_mod , only : shr_sys_abort use shr_cal_mod , only : shr_cal_ymd2date - use shr_log_mod , only : shr_log_setLogUnit + use shr_log_mod , only : shr_log_setLogUnit use dshr_methods_mod , only : dshr_state_diagnose, chkerr, memcheck use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_advance, shr_strdata_init_from_config use dshr_mod , only : dshr_model_initphase, dshr_init, dshr_mesh_init @@ -58,6 +58,11 @@ module cdeps_docn_comp use docn_datamode_cplhist_mod , only : docn_datamode_cplhist_advance use docn_datamode_cplhist_mod , only : docn_datamode_cplhist_restart_read use docn_datamode_cplhist_mod , only : docn_datamode_cplhist_restart_write + use docn_datamode_multilev_mod , only : docn_datamode_multilev_advertise + use docn_datamode_multilev_mod , only : docn_datamode_multilev_init_pointers + use docn_datamode_multilev_mod , only : docn_datamode_multilev_advance + use docn_datamode_multilev_mod , only : docn_datamode_multilev_restart_read + use docn_datamode_multilev_mod , only : docn_datamode_multilev_restart_write use docn_import_data_mod , only : docn_import_data_advertise implicit none @@ -179,6 +184,7 @@ end subroutine SetServices !=============================================================================== subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) use shr_nl_mod, only: shr_nl_find_group_name + ! input/output variables type(ESMF_GridComp) :: gcomp type(ESMF_State) :: importState, exportState @@ -300,7 +306,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) trim(datamode) == 'som_aquap' .or. & ! read stream, needs import data trim(datamode) == 'cplhist' .or. & ! read stream, needs import data trim(datamode) == 'sst_aquap_analytic' .or. & ! analytic, no streams, import or export data - trim(datamode) == 'sst_aquap_constant' ) then ! analytic, no streams, import or export data + trim(datamode) == 'sst_aquap_constant' .or. & ! analytic, no streams, import or export data + trim(datamode) == 'multilev') then ! multilevel ocean input ! success do nothing else call shr_sys_abort(' ERROR illegal docn datamode = '//trim(datamode)) @@ -323,6 +330,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) else if (trim(datamode) == 'cplhist') then call docn_datamode_cplhist_advertise(exportState, fldsExport, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + else if (trim(datamode) == 'multilev') then + call docn_datamode_multilev_advertise(exportState, fldsExport, flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end if if (trim(import_data_fields) /= 'none') then @@ -550,6 +560,9 @@ subroutine docn_comp_run(importState, exportState, clock, target_ymd, target_tod case('cplhist') call docn_datamode_cplhist_init_pointers(exportState, model_frac, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('multilev') + call docn_datamode_multilev_init_pointers(exportState, sdat, model_frac, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select ! Read restart if needed @@ -607,6 +620,9 @@ subroutine docn_comp_run(importState, exportState, clock, target_ymd, target_tod case('cplhist') call docn_datamode_cplhist_advance(rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + case('multilev') + call docn_datamode_multilev_advance(sdat, logunit, mainproc, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return end select ! Write restarts if needed (no restarts for aquaplanet analytic or aquaplanet input file) @@ -650,8 +666,10 @@ subroutine docn_init_dfields(importState, exportState, rc) ! local variables integer :: n integer :: fieldcount + integer :: dimcount type(ESMF_Field) :: lfield character(ESMF_MAXSTR) ,pointer :: lfieldnamelist(:) + character(ESMF_MAXSTR) :: fieldname(1) character(*), parameter :: subName = "(docn_init_dfields) " !------------------------------------------------------------------------------- @@ -668,9 +686,18 @@ subroutine docn_init_dfields(importState, exportState, rc) call ESMF_StateGet(exportState, itemName=trim(lfieldNameList(n)), field=lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (trim(lfieldnamelist(n)) /= flds_scalar_name) then - call dshr_dfield_add( dfields, sdat, trim(lfieldnamelist(n)), trim(lfieldnamelist(n)), exportState, & - logunit, mainproc, rc) + call ESMF_FieldGet(lfield, dimcount=dimCount, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (dimcount == 2) then + fieldname(1) = trim(lfieldnamelist(n)) + call dshr_dfield_add( dfields, sdat, trim(lfieldnamelist(n)), fieldname, exportState, & + logunit, mainproc, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + else + call dshr_dfield_add( dfields, sdat, trim(lfieldnamelist(n)), trim(lfieldnamelist(n)), exportState, & + logunit, mainproc, rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + endif end if end do end subroutine docn_init_dfields diff --git a/dshr/dshr_dfield_mod.F90 b/dshr/dshr_dfield_mod.F90 index 2c2862042..2f3b33d69 100644 --- a/dshr/dshr_dfield_mod.F90 +++ b/dshr/dshr_dfield_mod.F90 @@ -1,7 +1,7 @@ module dshr_dfield_mod use ESMF , only : ESMF_State, ESMF_FieldBundle, ESMF_MAXSTR, ESMF_SUCCESS - use ESMF , only : ESMF_FieldBundleGet, ESMF_ITEMORDER_ADDORDER, ESMF_Field + use ESMF , only : ESMF_FieldBundleGet, ESMF_ITEMORDER_ADDORDER, ESMF_Field, ESMF_FieldGet use shr_kind_mod , only : r8=>shr_kind_r8, cs=>shr_kind_cs, cl=>shr_kind_cl, cxx=>shr_kind_cxx use shr_sys_mod , only : shr_sys_abort use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_get_stream_count, shr_strdata_get_stream_fieldbundle @@ -438,9 +438,12 @@ subroutine dshr_dfield_copy(dfields, sdat, rc) type(ESMF_field) :: lfield type(dfield_type), pointer :: dfield real(r8), pointer :: data1d(:) + real(r8), pointer :: data2d(:,:) integer :: nf integer :: fldbun_index integer :: stream_index + integer :: ungriddedUBound_output(1) + integer :: ungriddedCount !------------------------------------------------------------------------------- rc = ESMF_SUCCESS @@ -464,13 +467,22 @@ subroutine dshr_dfield_copy(dfields, sdat, rc) do nf = 1,size(dfield%stream_indices) stream_index = dfield%stream_indices(nf) fldbun_index = dfield%fldbun_indices(nf) - if(stream_index > 0) then + if (stream_index > 0) then fldbun_model = shr_strdata_get_stream_fieldbundle(sdat, stream_index, 'model') call dshr_fldbun_getfieldn(fldbun_model, fldbun_index, lfield, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call dshr_field_getfldptr(lfield, fldptr1=data1d, rc=rc) + call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound_output, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - dfield%state_data2d(nf,:) = data1d(:) + ungriddedCount = ungriddedUBound_output(1) + if (ungriddedCount > 0) then + call dshr_field_getfldptr(lfield, fldptr2=data2d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + dfield%state_data2d(:,:) = data2d(:,:) + else + call dshr_field_getfldptr(lfield, fldptr1=data1d, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + dfield%state_data2d(nf,:) = data1d(:) + end if endif end do end if diff --git a/dshr/dshr_mod.F90 b/dshr/dshr_mod.F90 index 977b9d439..faa346b68 100644 --- a/dshr/dshr_mod.F90 +++ b/dshr/dshr_mod.F90 @@ -60,9 +60,9 @@ module dshr_mod public :: dshr_state_setscalar public :: dshr_orbital_update public :: dshr_orbital_init + public :: dshr_alarm_init ! initialize alarms private :: dshr_mesh_create_scol ! create mesh for single column mode - private :: dshr_alarm_init ! initialize alarms private :: dshr_time_init ! initialize time ! Note that gridTofieldMap = 2, therefore the ungridded dimension is innermost diff --git a/streams/dshr_methods_mod.F90 b/streams/dshr_methods_mod.F90 index 8e721bed1..0ce787596 100644 --- a/streams/dshr_methods_mod.F90 +++ b/streams/dshr_methods_mod.F90 @@ -568,7 +568,7 @@ subroutine dshr_field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return if (ungriddedUBound(1) > 0) then if (.not.present(fldptr2)) then - call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array ", & + call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array for "//trim(name), & ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) rc = ESMF_FAILURE return @@ -578,7 +578,7 @@ subroutine dshr_field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) lrank = 2 else if (.not.present(fldptr1)) then - call ESMF_LogWrite(trim(subname)//": ERROR missing rank=1 array ", & + call ESMF_LogWrite(trim(subname)//": ERROR missing rank=1 array for "//trim(name), & ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) rc = ESMF_FAILURE return diff --git a/streams/dshr_strdata_mod.F90 b/streams/dshr_strdata_mod.F90 index e9159f69b..60e3d0ddc 100644 --- a/streams/dshr_strdata_mod.F90 +++ b/streams/dshr_strdata_mod.F90 @@ -51,7 +51,7 @@ module dshr_strdata_mod use pio , only : pio_inquire, pio_inq_varid, pio_inq_varndims, pio_inq_vardimid use pio , only : pio_inq_dimlen, pio_inq_vartype, pio_inq_dimname, pio_inq_dimid use pio , only : pio_double, pio_real, pio_int, pio_offset_kind, pio_get_var - use pio , only : pio_read_darray, pio_setframe, pio_fill_double, pio_get_att + use pio , only : pio_read_darray, pio_setframe, pio_fill_double, pio_get_att, pio_inq_att use pio , only : PIO_BCAST_ERROR, PIO_RETURN_ERROR, PIO_NOERR, PIO_INTERNAL_ERROR, PIO_SHORT implicit none @@ -94,6 +94,7 @@ module dshr_strdata_mod character(CL), allocatable :: fldlist_stream(:) ! names of stream file fields character(CL), allocatable :: fldlist_model(:) ! names of stream model fields integer :: stream_nlev ! number of vertical levels in stream + real(r8), allocatable :: stream_vlevs(:) ! values of vertical levels in stream integer :: stream_lb ! index of the Lowerbound (LB) in fldlist_stream integer :: stream_ub ! index of the Upperbound (UB) in fldlist_stream type(ESMF_Field) :: field_stream ! a field on the stream data domain @@ -679,7 +680,10 @@ subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc) integer :: rcode character(CX) :: filename integer :: dimid + type(var_desc_t) :: varid integer :: stream_nlev + integer :: old_handle ! previous setting of pio error handling + character(CS) :: units character(*), parameter :: subname = '(shr_strdata_set_stream_domain) ' ! ---------------------------------------------- @@ -698,10 +702,25 @@ subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc) rcode = pio_openfile(sdat%pio_subsystem, pioid, sdat%io_type, trim(filename), pio_nowrite) rcode = pio_inq_dimid(pioid, trim(sdat%stream(stream_index)%lev_dimname), dimid) rcode = pio_inq_dimlen(pioid, dimid, stream_nlev) + allocate(sdat%pstrm(stream_index)%stream_vlevs(stream_nlev)) + rcode = pio_inq_varid(pioid, trim(sdat%stream(stream_index)%lev_dimname), varid) + rcode = pio_get_var(pioid, varid, sdat%pstrm(stream_index)%stream_vlevs) + + ! Determine vertical coordinates units - assume that default is m + call pio_seterrorhandling(pioid, PIO_BCAST_ERROR, old_handle) + rcode = pio_inq_att(pioid, varid, 'units') + call pio_seterrorhandling(pioid, old_handle) + if (rcode == PIO_NOERR) then + rcode = pio_get_att(pioid, varid, 'units', units) + if (trim(units) == 'centimeters' .or. trim(units) == 'cm') then + sdat%pstrm(stream_index)%stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) / 100. + end if + end if call pio_closefile(pioid) end if if (sdat%mainproc) then write(sdat%stream(1)%logunit,*) trim(subname)//' stream_nlev = ',stream_nlev + write(sdat%stream(1)%logunit,*)' stream vertical levels = ',sdat%pstrm(stream_index)%stream_vlevs end if ! Set stream_nlev in the per-stream sdat info