From c12a8d0c347e66e056c5f01f068fa563a05852ee Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 17 Apr 2024 11:04:49 -0600 Subject: [PATCH 01/15] add Fogg_rofi to dglc --- dglc/dglc_datamode_noevolve_mod.F90 | 261 +++++++++++++++------------- dglc/glc_comp_nuopc.F90 | 201 +++++++++++++++++++-- dshr/dshr_mod.F90 | 2 +- 3 files changed, 328 insertions(+), 136 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 9a5431e6..95a6d846 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -39,6 +39,7 @@ 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 type(icesheet_ptr_t), allocatable :: Sl_tsrf(:) @@ -50,6 +51,7 @@ 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' @@ -96,6 +98,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 @@ -150,6 +153,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 +166,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 +175,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 +231,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 233794fd..4597c406 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') @@ -449,11 +454,14 @@ subroutine ModelAdvance(gcomp, rc) integer :: yr ! year integer :: mon ! month integer :: day ! day in month + integer :: tod ! seconds in day 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 +471,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 +489,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 +516,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 +589,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 +658,151 @@ 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 + character(len=128) :: name + 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/dshr/dshr_mod.F90 b/dshr/dshr_mod.F90 index 977b9d43..faa346b6 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 From 6673b85419b09b7926ec35160189ba004ed12f44 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 17 Apr 2024 14:12:54 -0600 Subject: [PATCH 02/15] update to have antarctica work with clm --- dglc/cime_config/config_component.xml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dglc/cime_config/config_component.xml b/dglc/cime_config/config_component.xml index 15271a5b..a47b9d5e 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 From 1e9220b746984380e791490f3609c55b2a0628c6 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Fri, 19 Apr 2024 08:07:47 -0600 Subject: [PATCH 03/15] added ability to read in multi-level docn data --- docn/CMakeLists.txt | 1 + docn/cime_config/config_component.xml | 6 +- docn/cime_config/namelist_definition_docn.xml | 4 +- docn/cime_config/stream_definition_docn.xml | 65 ++++++++++++++++++- docn/ocn_comp_nuopc.F90 | 35 ++++++++-- dshr/dshr_dfield_mod.F90 | 7 +- streams/dshr_methods_mod.F90 | 4 +- 7 files changed, 109 insertions(+), 13 deletions(-) diff --git a/docn/CMakeLists.txt b/docn/CMakeLists.txt index 007d595c..256dc304 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 3a98b787..ab88c6fe 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 a44ae344..71109b1c 100644 --- a/docn/cime_config/namelist_definition_docn.xml +++ b/docn/cime_config/namelist_definition_docn.xml @@ -32,6 +32,7 @@ '' '' '' + sst_depth,salinity_depth @@ -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 0dd4baf5..4144e0e5 100644 --- a/docn/cime_config/stream_definition_docn.xml +++ b/docn/cime_config/stream_definition_docn.xml @@ -206,6 +206,69 @@ 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 + diff --git a/docn/ocn_comp_nuopc.F90 b/docn/ocn_comp_nuopc.F90 index 4bdeb81f..3a9c12a8 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 @@ -193,6 +198,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: bcasttmp(4) real(r8) :: rtmp(1) type(ESMF_VM) :: vm + integer :: nlev = 60 !DEBUG - remove this and put into namelist character(len=*),parameter :: subname=trim(module_name)//':(InitializeAdvertise) ' character(*) ,parameter :: F00 = "('(" // trim(module_name) // ") ',8a)" character(*) ,parameter :: F01 = "('(" // trim(module_name) // ") ',a,2x,i8)" @@ -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, nlev, 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, 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(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 2c286204..bfd5cdc1 100644 --- a/dshr/dshr_dfield_mod.F90 +++ b/dshr/dshr_dfield_mod.F90 @@ -438,6 +438,7 @@ 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 @@ -464,13 +465,13 @@ 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 dshr_field_getfldptr(lfield, fldptr2=data2d, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - dfield%state_data2d(nf,:) = data1d(:) + dfield%state_data2d(:,:) = data2d(:,:) endif end do end if diff --git a/streams/dshr_methods_mod.F90 b/streams/dshr_methods_mod.F90 index 8e721bed..0ce78759 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 From 40d63eb1d271f52553f525a43e0d420c01527fa9 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Fri, 19 Apr 2024 08:45:12 -0600 Subject: [PATCH 04/15] added new multilevel mode to docn --- docn/docn_datamode_multilev_mod.F90 | 159 ++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 docn/docn_datamode_multilev_mod.F90 diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 new file mode 100644 index 00000000..fa02eb26 --- /dev/null +++ b/docn/docn_datamode_multilev_mod.F90 @@ -0,0 +1,159 @@ +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 + 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_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 + + ! 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() + + ! 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, nlev, rc) + + ! input/output variables + type(esmf_State) , intent(inout) :: exportState + type(fldlist_type) , pointer :: fldsexport + character(len=*) , intent(in) :: flds_scalar_name + integer , intent(in) :: nlev + 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) + call dshr_fldList_add(fldsExport, 'So_s_depth', ungridded_lbound=1, ungridded_ubound=nlev) + + 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, ocn_fraction, rc) + + ! input/output variables + type(ESMF_State) , intent(inout) :: exportState + 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 export fields + 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(rc) + + ! input/output variables + integer, intent(out) :: rc + + ! local variables + integer :: i + character(len=*), parameter :: subname='(docn_datamode_multilev): ' + !------------------------------------------------------------------------------- + + rc = ESMF_SUCCESS + + do i = 1,size(So_omask) + if (So_omask(i) == 0.) then + So_t_depth(:,i) = 1.e30 + So_s_depth(:,i) = 1.e30 + end if + end do + + 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 From a6f0eb0be476d11692bc8c147c8ba64a557d1d32 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Tue, 23 Apr 2024 05:34:52 -0600 Subject: [PATCH 05/15] changes to get multilevel ocean input read --- docn/docn_datamode_multilev_mod.F90 | 103 +++++++++++++++++++++++----- docn/ocn_comp_nuopc.F90 | 8 +-- streams/dshr_strdata_mod.F90 | 9 +++ 3 files changed, 97 insertions(+), 23 deletions(-) diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 index fa02eb26..c43f6b35 100644 --- a/docn/docn_datamode_multilev_mod.F90 +++ b/docn/docn_datamode_multilev_mod.F90 @@ -2,12 +2,12 @@ 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 + 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_type + use dshr_strdata_mod , only : shr_strdata_get_stream_pointer, shr_strdata_type implicit none private ! except @@ -18,11 +18,23 @@ module docn_datamode_multilev_mod public :: docn_datamode_multilev_restart_read public :: docn_datamode_multilev_restart_write - ! export fields - real(r8), pointer :: So_omask(:) => null() ! real ocean fraction sent to mediator + ! 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) = (/ & + 60. , 120. , 180. , 240. , 300. , 360. , & + 420. , 480. , 540. , 600. , 660. , 720. , & + 780. , 840. , 900. , 960. , 1020., 1080., & + 1140., 1200., 1260., 1320., 1380., 1440., & + 1500., 1560., 1620., 1680., 1740., 1800. /) + ! constants character(*) , parameter :: nullstr = 'null' character(*) , parameter :: rpfile = 'rpointer.ocn' @@ -33,13 +45,12 @@ module docn_datamode_multilev_mod contains !=============================================================================== - subroutine docn_datamode_multilev_advertise(exportState, fldsexport, flds_scalar_name, nlev, rc) + 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(in) :: nlev integer , intent(out) :: rc ! local variables @@ -51,8 +62,8 @@ subroutine docn_datamode_multilev_advertise(exportState, fldsexport, flds_scalar ! 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) - call dshr_fldList_add(fldsExport, 'So_s_depth', ungridded_lbound=1, ungridded_ubound=nlev) + 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)) @@ -65,12 +76,13 @@ subroutine docn_datamode_multilev_advertise(exportState, fldsexport, flds_scalar end subroutine docn_datamode_multilev_advertise !=============================================================================== - subroutine docn_datamode_multilev_init_pointers(exportState, ocn_fraction, rc) + subroutine docn_datamode_multilev_init_pointers(exportState, sdat, ocn_fraction, rc) ! input/output variables - type(ESMF_State) , intent(inout) :: exportState - real(r8) , intent(in) :: ocn_fraction(:) - integer , intent(out) :: rc + 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): ' @@ -78,7 +90,15 @@ subroutine docn_datamode_multilev_init_pointers(exportState, ocn_fraction, rc) 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) @@ -96,25 +116,70 @@ subroutine docn_datamode_multilev_init_pointers(exportState, ocn_fraction, rc) end subroutine docn_datamode_multilev_init_pointers !=============================================================================== - subroutine docn_datamode_multilev_advance(rc) + subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) ! input/output variables - integer, intent(out) :: rc + type(shr_strdata_type) , intent(in) :: sdat + integer , intent(in) :: logunit + logical , intent(in) :: mainproc + integer , intent(out) :: rc ! local variables - integer :: i + 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 - do i = 1,size(So_omask) - if (So_omask(i) == 0.) then - So_t_depth(:,i) = 1.e30 - So_s_depth(:,i) = 1.e30 + ! 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) = 1.e30 + So_s_depth(ko,i) = 1.e30 + else + if (stream_So_t_depth(ki+1,i) == shr_const_spval) then + So_t_depth(ko,i) = stream_So_t_depth(ki,i) + So_s_depth(ko,i) = stream_So_s_depth(ki,i) + 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 + + 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 !=============================================================================== diff --git a/docn/ocn_comp_nuopc.F90 b/docn/ocn_comp_nuopc.F90 index 3a9c12a8..ced8cb86 100644 --- a/docn/ocn_comp_nuopc.F90 +++ b/docn/ocn_comp_nuopc.F90 @@ -184,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 @@ -198,7 +199,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: bcasttmp(4) real(r8) :: rtmp(1) type(ESMF_VM) :: vm - integer :: nlev = 60 !DEBUG - remove this and put into namelist character(len=*),parameter :: subname=trim(module_name)//':(InitializeAdvertise) ' character(*) ,parameter :: F00 = "('(" // trim(module_name) // ") ',8a)" character(*) ,parameter :: F01 = "('(" // trim(module_name) // ") ',a,2x,i8)" @@ -331,7 +331,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) 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, nlev, rc) + call docn_datamode_multilev_advertise(exportState, fldsExport, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -561,7 +561,7 @@ subroutine docn_comp_run(importState, exportState, clock, target_ymd, target_tod 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, model_frac, rc) + call docn_datamode_multilev_init_pointers(exportState, sdat, model_frac, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end select @@ -621,7 +621,7 @@ subroutine docn_comp_run(importState, exportState, clock, target_ymd, target_tod call docn_datamode_cplhist_advance(rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return case('multilev') - call docn_datamode_multilev_advance(rc=rc) + call docn_datamode_multilev_advance(sdat, logunit, mainproc, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end select diff --git a/streams/dshr_strdata_mod.F90 b/streams/dshr_strdata_mod.F90 index e9159f69..5830f5e9 100644 --- a/streams/dshr_strdata_mod.F90 +++ b/streams/dshr_strdata_mod.F90 @@ -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,6 +680,7 @@ 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 character(*), parameter :: subname = '(shr_strdata_set_stream_domain) ' ! ---------------------------------------------- @@ -698,10 +700,17 @@ 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) + ! DEBUG: input is in cm + sdat%pstrm(stream_index)%stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) / 100. + ! DEBUG 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 From 458f0df591a569ccbc5b1505ea8fed85d1013144 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Tue, 23 Apr 2024 05:53:49 -0600 Subject: [PATCH 06/15] added query for units of ocean depth layers --- docn/docn_datamode_multilev_mod.F90 | 3 ++- streams/dshr_strdata_mod.F90 | 18 ++++++++++++++---- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 index c43f6b35..bb528a58 100644 --- a/docn/docn_datamode_multilev_mod.F90 +++ b/docn/docn_datamode_multilev_mod.F90 @@ -142,7 +142,8 @@ subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) nlev_stream = sdat%pstrm(stream_index)%stream_nlev allocate(stream_vlevs(nlev_stream)) - stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) + ! TODO: for now hard-wired input in cm + stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) / 100. do ko = 1,nlev_export level_found = .false. diff --git a/streams/dshr_strdata_mod.F90 b/streams/dshr_strdata_mod.F90 index 5830f5e9..60e3d0dd 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 @@ -682,6 +682,8 @@ subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc) 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) ' ! ---------------------------------------------- @@ -703,9 +705,17 @@ subroutine shr_strdata_get_stream_nlev(sdat, stream_index, rc) 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) - ! DEBUG: input is in cm - sdat%pstrm(stream_index)%stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) / 100. - ! DEBUG + + ! 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 From 38ce4fbca00f964d75cbc9cfa03331cbdf8727bc Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 24 Apr 2024 09:37:48 -0600 Subject: [PATCH 07/15] updates needed for multi-level ocean input --- dglc/dglc_datamode_noevolve_mod.F90 | 11 ++++++++++- docn/docn_datamode_multilev_mod.F90 | 10 +++++----- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 95a6d846..08453147 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 @@ -42,8 +43,11 @@ module dglc_datamode_noevolve_mod 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' @@ -56,6 +60,8 @@ module dglc_datamode_noevolve_mod ! 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 = & @@ -115,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 diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 index bb528a58..6a1ec760 100644 --- a/docn/docn_datamode_multilev_mod.F90 +++ b/docn/docn_datamode_multilev_mod.F90 @@ -142,8 +142,7 @@ subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) nlev_stream = sdat%pstrm(stream_index)%stream_nlev allocate(stream_vlevs(nlev_stream)) - ! TODO: for now hard-wired input in cm - stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) / 100. + stream_vlevs(:) = sdat%pstrm(stream_index)%stream_vlevs(:) do ko = 1,nlev_export level_found = .false. @@ -157,15 +156,16 @@ subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) level_found = .true. do i = 1,size(So_omask) if (So_omask(i) == 0.) then - So_t_depth(ko,i) = 1.e30 - So_s_depth(ko,i) = 1.e30 + So_t_depth(ko,i) = shr_const_spval + So_s_depth(ko,i) = shr_const_spval else if (stream_So_t_depth(ki+1,i) == shr_const_spval) then - So_t_depth(ko,i) = stream_So_t_depth(ki,i) + 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) 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 From 2121f2b88109415554cfa68022faa1e35c60f8c9 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 24 Apr 2024 09:46:06 -0600 Subject: [PATCH 08/15] update comment --- docn/docn_datamode_multilev_mod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 index 6a1ec760..6a982866 100644 --- a/docn/docn_datamode_multilev_mod.F90 +++ b/docn/docn_datamode_multilev_mod.F90 @@ -159,6 +159,7 @@ subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) 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) == shr_const_spval) then 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) From 7e1096f4a1630210d3fca454d6f9a08563333381 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Mon, 29 Apr 2024 04:22:41 -0600 Subject: [PATCH 09/15] fixes to permit stream field with ungridded dimensions --- dshr/dshr_dfield_mod.F90 | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/dshr/dshr_dfield_mod.F90 b/dshr/dshr_dfield_mod.F90 index bfd5cdc1..2f3b33d6 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 @@ -442,6 +442,8 @@ subroutine dshr_dfield_copy(dfields, sdat, rc) integer :: nf integer :: fldbun_index integer :: stream_index + integer :: ungriddedUBound_output(1) + integer :: ungriddedCount !------------------------------------------------------------------------------- rc = ESMF_SUCCESS @@ -469,9 +471,18 @@ subroutine dshr_dfield_copy(dfields, sdat, rc) 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, fldptr2=data2d, rc=rc) + call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound_output, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - dfield%state_data2d(:,:) = data2d(:,:) + 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 From 1ba9b8f1716255b20469b2d0de93a563454d310f Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Mon, 29 Apr 2024 08:30:04 -0600 Subject: [PATCH 10/15] updated testlist so that DGLC%NOEVOLVE passes --- dglc/cime_config/testdefs/testlist_dglc.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dglc/cime_config/testdefs/testlist_dglc.xml b/dglc/cime_config/testdefs/testlist_dglc.xml index 194773fa..8383469c 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 @@ - + From dc6923e09a0959c7b3ec05ba0e07e2501d81e906 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 1 May 2024 14:29:38 -0600 Subject: [PATCH 11/15] fixes for vertical interpolation --- docn/docn_datamode_multilev_mod.F90 | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/docn/docn_datamode_multilev_mod.F90 b/docn/docn_datamode_multilev_mod.F90 index 6a982866..0c6ea681 100644 --- a/docn/docn_datamode_multilev_mod.F90 +++ b/docn/docn_datamode_multilev_mod.F90 @@ -29,10 +29,8 @@ module docn_datamode_multilev_mod integer, parameter :: nlev_export = 30 real(r8) :: vertical_levels(nlev_export) = (/ & - 60. , 120. , 180. , 240. , 300. , 360. , & - 420. , 480. , 540. , 600. , 660. , 720. , & - 780. , 840. , 900. , 960. , 1020., 1080., & - 1140., 1200., 1260., 1320., 1380., 1440., & + 60. , 120. , 180. , 240. , 300. , 360. , 420. , 480. , 540. , 600. , 660. , 720. , & + 780. , 840. , 900. , 960. , 1020., 1080., 1140., 1200., 1260., 1320., 1380., 1440., & 1500., 1560., 1620., 1680., 1740., 1800. /) ! constants @@ -160,9 +158,14 @@ subroutine docn_datamode_multilev_advance(sdat, logunit, mainproc, rc) 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) == shr_const_spval) then - 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) + 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 From 27726a1cab6aa0a3d863797fefb47e8d5a66671e Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 2 May 2024 19:34:21 +0200 Subject: [PATCH 12/15] update for blom output stream --- docn/cime_config/namelist_definition_docn.xml | 2 +- docn/cime_config/stream_definition_docn.xml | 45 ++++++++++++++++--- docn/docn_datamode_multilev_mod.F90 | 6 +-- 3 files changed, 44 insertions(+), 9 deletions(-) diff --git a/docn/cime_config/namelist_definition_docn.xml b/docn/cime_config/namelist_definition_docn.xml index 71109b1c..b2f78811 100644 --- a/docn/cime_config/namelist_definition_docn.xml +++ b/docn/cime_config/namelist_definition_docn.xml @@ -32,7 +32,7 @@ '' '' '' - sst_depth,salinity_depth + sst_salinity_depth_blom diff --git a/docn/cime_config/stream_definition_docn.xml b/docn/cime_config/stream_definition_docn.xml index 4144e0e5..4e4dd665 100644 --- a/docn/cime_config/stream_definition_docn.xml +++ b/docn/cime_config/stream_definition_docn.xml @@ -209,12 +209,12 @@ single - + - $DIN_LOC_ROOT/share/meshes/gx1v7_151008_ESMFmesh.nc + $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 + /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 @@ -242,10 +242,10 @@ - $DIN_LOC_ROOT/share/meshes/gx1v7_151008_ESMFmesh.nc + $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 + /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 @@ -271,4 +271,39 @@ 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 index 0c6ea681..897797f3 100644 --- a/docn/docn_datamode_multilev_mod.F90 +++ b/docn/docn_datamode_multilev_mod.F90 @@ -29,9 +29,9 @@ module docn_datamode_multilev_mod integer, parameter :: nlev_export = 30 real(r8) :: vertical_levels(nlev_export) = (/ & - 60. , 120. , 180. , 240. , 300. , 360. , 420. , 480. , 540. , 600. , 660. , 720. , & - 780. , 840. , 900. , 960. , 1020., 1080., 1140., 1200., 1260., 1320., 1380., 1440., & - 1500., 1560., 1620., 1680., 1740., 1800. /) + 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' From 017e6a992b61488347ff64dc616035506bd405c1 Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Tue, 14 May 2024 07:40:34 -0600 Subject: [PATCH 13/15] update workflow versions --- .github/workflows/extbuild.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/extbuild.yml b/.github/workflows/extbuild.yml index 52afb7b8..185f6d85 100644 --- a/.github/workflows/extbuild.yml +++ b/.github/workflows/extbuild.yml @@ -19,11 +19,11 @@ jobs: CPPFLAGS: "-I/usr/include -I/usr/local/include " LDFLAGS: "-L/usr/lib/x86_64-linux-gnu " # Versions of all dependencies can be updated here - these match tag names in the github repo - ESMF_VERSION: v8.5.0 - ParallelIO_VERSION: pio2_6_0 + ESMF_VERSION: v8.6.1 + ParallelIO_VERSION: pio2_6_2 steps: - id: checkout-CDEPS - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: submodules: recursive - id: load-env @@ -37,7 +37,7 @@ jobs: sudo apt-get install autotools-dev autoconf - name: Cache PARALLELIO id: cache-PARALLELIO - uses: actions/cache@v3 + uses: actions/cache@v4 with: path: ${GITHUB_WORKSPACE}/pio key: ${{ runner.os }}-${{ env.ParallelIO_VERSION }}-parallelio2 From 472b568fadfcfb744ea4cdaf805af7b635238d06 Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Tue, 14 May 2024 07:50:55 -0600 Subject: [PATCH 14/15] comment out unused vars --- dglc/dglc_datamode_noevolve_mod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 08453147..e4392105 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -46,8 +46,8 @@ module dglc_datamode_noevolve_mod 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(:) +! 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' From a8baf6644e46ceaeb47fbf03ef47c4a4e61b9753 Mon Sep 17 00:00:00 2001 From: Jim Edwards Date: Tue, 14 May 2024 07:54:57 -0600 Subject: [PATCH 15/15] remove more unused vars --- dglc/glc_comp_nuopc.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/dglc/glc_comp_nuopc.F90 b/dglc/glc_comp_nuopc.F90 index 4597c406..09c858ff 100644 --- a/dglc/glc_comp_nuopc.F90 +++ b/dglc/glc_comp_nuopc.F90 @@ -454,7 +454,6 @@ subroutine ModelAdvance(gcomp, rc) integer :: yr ! year integer :: mon ! month integer :: day ! day in month - integer :: tod ! seconds in day logical :: restart_write type(ESMF_Alarm) :: valid_alarm logical :: valid_inputs ! if true, inputs from mediator are valid @@ -680,7 +679,6 @@ subroutine ModelSetRunClock(gcomp, rc) integer :: stop_n ! Number until stop interval integer :: stop_ymd ! Stop date (YYYYMMDD) type(ESMF_ALARM) :: stop_alarm - character(len=128) :: name integer :: alarmcount character(len=CS) :: glc_avg_period ! averaging period in mediator type(ESMF_ALARM) :: valid_alarm ! model alarm