From aecb0f6f16c2c8f5da9a0ef29024d3fb3550e94f Mon Sep 17 00:00:00 2001 From: David Bailey Date: Mon, 7 Nov 2022 18:27:28 -0700 Subject: [PATCH 1/3] CESM updates --- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 181 ++++++++++-------- .../drivers/nuopc/cmeps/ice_import_export.F90 | 119 ++++++------ cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 | 54 +++++- 3 files changed, 207 insertions(+), 147 deletions(-) diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 182308973..4d8210232 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -17,12 +17,13 @@ module ice_comp_nuopc use NUOPC_Model , only : NUOPC_ModelGet, SetVM use ice_constants , only : ice_init_constants, c0 use ice_shr_methods , only : chkerr, state_setscalar, state_getscalar, state_diagnose, alarmInit - use ice_shr_methods , only : set_component_logging, get_component_instance, state_flddebug + use ice_shr_methods , only : get_component_instance, state_flddebug + use ice_import_export , only : ice_import, ice_export, ice_advertise_fields, ice_realize_fields use ice_domain_size , only : nx_global, ny_global use ice_grid , only : grid_type, init_grid2 use ice_communicate , only : init_communicate, my_task, master_task, mpi_comm_ice - use ice_calendar , only : force_restart_now, write_ic + use ice_calendar , only : force_restart_now, write_ic, init_calendar use ice_calendar , only : idate, mday, mmonth, myear, year_init use ice_calendar , only : msec, dt, calendar, calendar_type, nextsw_cday, istep use ice_calendar , only : ice_calendar_noleap, ice_calendar_gregorian @@ -40,7 +41,10 @@ module ice_comp_nuopc #ifdef CESMCOUPLED use shr_const_mod use shr_orb_mod , only : shr_orb_decl, shr_orb_params, SHR_ORB_UNDEF_REAL, SHR_ORB_UNDEF_INT - use ice_scam , only : scmlat, scmlon, scol_mask, scol_frac, scol_ni, scol_nj + use ice_scam , only : scmlat, scmlon, scol_mask, scol_frac, scol_ni, scol_nj, scol_area + use nuopc_shr_methods , only : set_component_logging +#else + use ice_shr_methods , only : set_component_logging #endif use ice_timers use CICE_InitMod , only : cice_init1, cice_init2 @@ -87,6 +91,15 @@ module ice_comp_nuopc integer :: dbug = 0 logical :: profile_memory = .false. + logical :: mastertask + integer :: start_ymd ! Start date (YYYYMMDD) + integer :: start_tod ! start time of day (s) + integer :: curr_ymd ! Current date (YYYYMMDD) + integer :: curr_tod ! Current time of day (s) + integer :: stop_ymd ! stop date (YYYYMMDD) + integer :: stop_tod ! stop time of day (sec) + integer :: ref_ymd ! Reference date (YYYYMMDD) + integer :: ref_tod ! reference time of day (s) integer , parameter :: debug_import = 0 ! internal debug level integer , parameter :: debug_export = 0 ! internal debug level character(*), parameter :: modName = "(ice_comp_nuopc)" @@ -201,12 +214,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) real(kind=dbl_kind) :: atmiter_conv_driver integer (kind=int_kind) :: natmiter integer (kind=int_kind) :: natmiter_driver - character(len=char_len) :: tfrz_option_driver ! tfrz_option from driver attributes - character(len=char_len) :: tfrz_option ! tfrz_option from cice namelist - integer(int_kind) :: ktherm integer :: localPet integer :: npes - logical :: mastertask type(ESMF_VM) :: vm integer :: lmpicom ! local communicator type(ESMF_Time) :: currTime ! Current time @@ -215,14 +224,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) type(ESMF_Time) :: refTime ! Ref time type(ESMF_TimeInterval) :: timeStep ! Model timestep type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type - integer :: start_ymd ! Start date (YYYYMMDD) - integer :: start_tod ! start time of day (s) - integer :: curr_ymd ! Current date (YYYYMMDD) - integer :: curr_tod ! Current time of day (s) - integer :: stop_ymd ! stop date (YYYYMMDD) - integer :: stop_tod ! stop time of day (sec) - integer :: ref_ymd ! Reference date (YYYYMMDD) - integer :: ref_tod ! reference time of day (s) integer :: yy,mm,dd ! Temporaries for time query integer :: dtime ! time step integer :: shrlogunit ! original log unit @@ -232,12 +233,15 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: iblk, jblk ! indices integer :: ig, jg ! indices integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain + character(len=char_len_long) :: diag_filename = 'unset' character(len=char_len_long) :: logmsg character(len=char_len_long) :: single_column_lnd_domainfile real(dbl_kind) :: scol_lon real(dbl_kind) :: scol_lat real(dbl_kind) :: scol_spval + character(len=char_len) :: tfrz_option ! tfrz_option from cice namelist + character(len=char_len) :: tfrz_option_driver ! tfrz_option from cice namelist character(len=*), parameter :: subname=trim(modName)//':(InitializeAdvertise) ' !-------------------------------- @@ -541,12 +545,20 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call NUOPC_CompAttributeGet(gcomp, name='scol_nj', value=cvalue, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) scol_nj + call NUOPC_CompAttributeGet(gcomp, name='scol_area', value=cvalue, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + read(cvalue,*) scol_area call ice_mesh_create_scolumn(scmlon, scmlat, ice_mesh, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return scol_valid = (scol_mask == 1) if (.not. scol_valid) then + ! Read the cice namelist as part of the call to cice_init1 + ! Note that if single_column is true and scol_valid is not - will never get here + call t_startf ('cice_init1') + call cice_init1 + call t_stopf ('cice_init1') ! Advertise fields call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -559,13 +571,18 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) RETURN end if end if - ! Read the cice namelist as part of the call to cice_init1 ! Note that if single_column is true and scol_valid is not - will never get here - call t_startf ('cice_init1') call cice_init1 call t_stopf ('cice_init1') + + !----------------------------------------------------------------- + ! Advertise fields + !----------------------------------------------------------------- + call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! Form of ocean freezing temperature ! 'minus1p8' = -1.8 C @@ -675,6 +692,43 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) end if end if + call t_stopf ('cice_init_total') + + end subroutine InitializeAdvertise + + !=============================================================================== + + subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) + + ! Arguments + type(ESMF_GridComp) :: gcomp + type(ESMF_State) :: importState + type(ESMF_State) :: exportState + type(ESMF_Clock) :: clock + integer, intent(out) :: rc + + ! Local variables + integer :: n + integer :: fieldcount + type(ESMF_Field) :: lfield + character(len=char_len_long) :: cvalue + real(dbl_kind) :: scol_lon + real(dbl_kind) :: scol_lat + real(dbl_kind) :: scol_spval + real(dbl_kind), pointer :: fldptr1d(:) + real(dbl_kind), pointer :: fldptr2d(:,:) + integer :: rank + character(len=char_len) :: tfrz_option ! tfrz_option from cice namelist + integer(int_kind) :: ktherm + + character(len=char_len_long) :: single_column_lnd_domainfile + character(len=char_len_long) , pointer :: lfieldnamelist(:) => null() + character(len=*), parameter :: subname=trim(modName)//':(InitializeRealize) ' + !-------------------------------- + + rc = ESMF_SUCCESS + if (dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) + !---------------------------------------------------------------------------- ! Second cice initialization phase -after initializing grid info !---------------------------------------------------------------------------- @@ -685,29 +739,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call t_startf ('cice_init2') call cice_init2() call t_stopf ('cice_init2') - - !---------------------------------------------------------------------------- - ! reset shr logging to my log file - !---------------------------------------------------------------------------- - - call icepack_query_parameters(ktherm_out=ktherm) - call icepack_query_parameters(tfrz_option_out=tfrz_option) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & - file=__FILE__, line=__LINE__) - - ! Now write output to nu_diag - this must happen AFTER call to cice_init - if (mastertask) then - write(nu_diag,'(a,d21.14)') trim(subname)//' cice init nextsw_cday = ',nextsw_cday - write(nu_diag,'(a)') trim(subname)//' tfrz_option = '//trim(tfrz_option) - if (ktherm == 2 .and. trim(tfrz_option) /= 'mushy') then - write(nu_diag,*) trim(subname),' Warning: Using ktherm = 2 and tfrz_option = ', trim(tfrz_option) - endif - write(nu_diag,'(a )') trim(subname)//' inst_name = '//trim(inst_name) - write(nu_diag,'(a,i8 )') trim(subname)//' inst_index = ',inst_index - write(nu_diag,'(a )') trim(subname)//' inst_suffix = ',trim(inst_suffix) - endif - !--------------------------------------------------------------------------- ! use EClock to reset calendar information on initial start !--------------------------------------------------------------------------- @@ -758,6 +789,30 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) end if call calendar() ! update calendar info + + !---------------------------------------------------------------------------- + ! reset shr logging to my log file + !---------------------------------------------------------------------------- + + call icepack_query_parameters(ktherm_out=ktherm) + call icepack_query_parameters(tfrz_option_out=tfrz_option) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + ! Now write output to nu_diag - this must happen AFTER call to cice_init + if (mastertask) then + write(nu_diag,'(a,d21.14)') trim(subname)//' cice init nextsw_cday = ',nextsw_cday + write(nu_diag,'(a)') trim(subname)//' tfrz_option = '//trim(tfrz_option) + if (ktherm == 2 .and. trim(tfrz_option) /= 'mushy') then + write(nu_diag,*) trim(subname),' Warning: Using ktherm = 2 and tfrz_option = ', trim(tfrz_option) + endif + write(nu_diag,'(a )') trim(subname)//' inst_name = '//trim(inst_name) + write(nu_diag,'(a,i8 )') trim(subname)//' inst_index = ',inst_index + write(nu_diag,'(a )') trim(subname)//' inst_suffix = ',trim(inst_suffix) + endif + + if (write_ic) then call accum_hist(dt) ! write initial conditions end if @@ -769,50 +824,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ice_prescribed_init(clock, ice_mesh, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - !----------------------------------------------------------------- - ! Advertise fields - !----------------------------------------------------------------- - - ! NOTE: the advertise phase needs to be called after the ice - ! initialization since the number of ice categories is needed for - ! ice_fraction_n and mean_sw_pen_to_ocn_ifrac_n - call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - - call t_stopf ('cice_init_total') - - end subroutine InitializeAdvertise - - !=============================================================================== - - subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) - - ! Arguments - type(ESMF_GridComp) :: gcomp - type(ESMF_State) :: importState - type(ESMF_State) :: exportState - type(ESMF_Clock) :: clock - integer, intent(out) :: rc - - ! Local variables - integer :: n - integer :: fieldcount - type(ESMF_Field) :: lfield - character(len=char_len_long) :: cvalue - real(dbl_kind) :: scol_lon - real(dbl_kind) :: scol_lat - real(dbl_kind) :: scol_spval - real(dbl_kind), pointer :: fldptr1d(:) - real(dbl_kind), pointer :: fldptr2d(:,:) - integer :: rank - character(len=char_len_long) :: single_column_lnd_domainfile - character(len=char_len_long) , pointer :: lfieldnamelist(:) => null() - character(len=*), parameter :: subname=trim(modName)//':(InitializeRealize) ' - !-------------------------------- - - rc = ESMF_SUCCESS - if (dbug > 5) call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO) - #ifdef CESMCOUPLED ! if single column is not valid - set all export state fields to zero and return if (single_column .and. .not. scol_valid) then @@ -848,7 +859,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! *** RETURN HERE *** ! ******************* RETURN - else + else if(single_column) then write(nu_diag,'(a,3(f10.5,2x))')' (ice_comp_nuopc) single column mode lon/lat/frac is ',& scmlon,scmlat,scol_frac end if diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index d95a4d9b2..d1408d528 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -36,9 +36,9 @@ module ice_import_export use ice_shr_methods , only : chkerr, state_reset use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc , only : icepack_query_parameters, icepack_query_tracer_flags - use icepack_intfc , only : icepack_query_tracer_indices use icepack_intfc , only : icepack_liquidus_temperature use icepack_intfc , only : icepack_sea_freezing_temperature + use icepack_intfc , only : icepack_query_tracer_indices use icepack_parameters , only : puny, c2 use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf #ifdef CESMCOUPLED @@ -181,6 +181,13 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam if (flds_wiso) then call fldlist_add(fldsToIce_num, fldsToIce, 'So_roce_wiso', ungridded_lbound=1, ungridded_ubound=3) end if + ! the following are advertised but might not be connected if they are not advertised in the + ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific + ! from wave + if (flds_wave) then + call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, & + ungridded_ubound=25) + end if ! from atmosphere call fldlist_add(fldsToIce_num, fldsToIce, 'inst_height_lowest' ) @@ -208,14 +215,6 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam ! from atm - dry dust deposition fluxes (4 sizes) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_dstdry', ungridded_lbound=1, ungridded_ubound=4) - ! the following are advertised but might not be connected if they are not advertised in the - ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific - ! from wave - if (flds_wave) then - call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, & - ungridded_ubound=25) - end if - do n = 1,fldsToIce_num call NUOPC_Advertise(importState, standardName=fldsToIce(n)%stdname, & TransferOfferGeomObject='will provide', rc=rc) @@ -307,6 +306,7 @@ end subroutine ice_advertise_fields !============================================================================== subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc) + use ice_scam, only : single_column ! input/output variables type(ESMF_GridComp) :: gcomp @@ -320,7 +320,7 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc type(ESMF_State) :: exportState type(ESMF_Field) :: lfield integer :: numOwnedElements - integer :: i, j, iblk, n, k + integer :: i, j, iblk, n integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain type(block) :: this_block ! block information for current block real(dbl_kind), allocatable :: mesh_areas(:) @@ -361,10 +361,10 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc tag=subname//':CICE_Import',& mesh=mesh, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - #ifdef CESMCOUPLED ! Get mesh areas from second field - using second field since the ! first field is the scalar field + if (single_column) return call ESMF_MeshGet(mesh, numOwnedElements=numOwnedElements, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return @@ -435,6 +435,7 @@ subroutine ice_import( importState, rc ) integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain type(block) :: this_block ! block information for current block real (kind=dbl_kind),allocatable :: aflds(:,:,:,:) + real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: work real (kind=dbl_kind) :: workx, worky real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP real (kind=dbl_kind) :: Tffresh @@ -553,6 +554,29 @@ subroutine ice_import( importState, rc ) areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25) + if (State_FldChk(importState, 'Sw_elevation_spectrum')) then + if (nfreq /= 25) then + call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ") + end if + call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + do k = 1,nfreq + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n) + end do + end do + end do + end do + end if + ! perform a halo update if (.not.prescribed_ice) then @@ -586,29 +610,6 @@ subroutine ice_import( importState, rc ) end do !$OMP END PARALLEL DO - ! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25) - if (State_FldChk(importState, 'Sw_elevation_spectrum')) then - if (nfreq /= 25) then - call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ") - end if - call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - do k = 1,nfreq - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo; ihi = this_block%ihi - jlo = this_block%jlo; jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n) - end do - end do - end do - end do - end if - if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks @@ -889,6 +890,8 @@ end subroutine ice_import !=============================================================================== subroutine ice_export( exportState, rc ) + use ice_scam, only : single_column + ! input/output variables type(ESMF_State), intent(inout) :: exportState integer , intent(out) :: rc @@ -911,12 +914,13 @@ subroutine ice_export( exportState, rc ) real (kind=dbl_kind) :: ailohi(nx_block,ny_block,max_blocks) ! fractional ice area real (kind=dbl_kind) :: floediam(nx_block,ny_block,max_blocks) real (kind=dbl_kind) :: floethick(nx_block,ny_block,max_blocks) ! ice thickness - real (kind=dbl_kind) :: Tffresh logical (kind=log_kind) :: tr_fsd integer (kind=int_kind) :: nt_fsd + real (kind=dbl_kind) :: Tffresh real (kind=dbl_kind), allocatable :: tempfld(:,:,:) real (kind=dbl_kind), pointer :: dataptr_ifrac_n(:,:) real (kind=dbl_kind), pointer :: dataptr_swpen_n(:,:) + logical (kind=log_kind), save :: first_call = .true. character(len=*),parameter :: subname = 'ice_export' !----------------------------------------------------- @@ -963,6 +967,9 @@ subroutine ice_export( exportState, rc ) ! ice fraction ailohi(i,j,iblk) = min(aice(i,j,iblk), c1) + ! surface temperature + Tsrf(i,j,iblk) = Tffresh + trcr(i,j,1,iblk) !Kelvin (original ???) + if (tr_fsd) then ! floe thickness (m) if (aice(i,j,iblk) > puny) then @@ -984,9 +991,6 @@ subroutine ice_export( exportState, rc ) floediam(i,j,iblk) = MAX(c2*floe_rad_c(1),workx) endif - ! surface temperature - Tsrf(i,j,iblk) = Tffresh + trcr(i,j,1,iblk) !Kelvin (original ???) - ! wind stress (on POP T-grid: convert to lat-lon) workx = strairxT(i,j,iblk) ! N/m^2 worky = strairyT(i,j,iblk) ! N/m^2 @@ -1042,8 +1046,11 @@ subroutine ice_export( exportState, rc ) !--------------------------------- ! Zero out fields with tmask for proper coupler accumulation in ice free areas - call state_reset(exportState, c0, rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (first_call .or. .not.single_column) then + call state_reset(exportState, c0, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + first_call = .false. + endif ! Create a temporary field allocate(tempfld(nx_block,ny_block,nblocks)) @@ -1134,22 +1141,6 @@ subroutine ice_export( exportState, rc ) call state_setexport(exportState, 'Si_snowh' , input=tempfld , lmask=tmask, ifrac=ailohi, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! ------ - ! optional floe diameter and ice thickness to wave - ! ------ - - ! Sea ice thickness (m) - if (State_FldChk(exportState, 'Si_thick')) then - call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - ! Sea ice floe diameter (m) - if (State_FldChk(exportState, 'Si_floediam')) then - call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - ! ------ ! ice/atm fluxes computed by ice ! ------ @@ -1337,6 +1328,22 @@ subroutine ice_export( exportState, rc ) end do end if + ! ------ + ! optional floe diameter and ice thickness to wave + ! ------ + + ! Sea ice thickness (m) + if (State_FldChk(exportState, 'Si_thick')) then + call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Sea ice floe diameter (m) + if (State_FldChk(exportState, 'Si_floediam')) then + call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end subroutine ice_export !=============================================================================== diff --git a/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 index 0b1b9349c..a9b19df6b 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 @@ -427,7 +427,7 @@ subroutine ice_mesh_create_scolumn(scol_lon, scol_lat, ice_mesh, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Allocate module variable ocn_gridcell_frac - allocate(ocn_gridcell_frac(nx_block,ny_block,max_blocks)) + allocate(ocn_gridcell_frac(2,2,1)) ocn_gridcell_frac(:,:,:) = scol_frac end subroutine ice_mesh_create_scolumn @@ -560,7 +560,7 @@ subroutine ice_mesh_check(gcomp, ice_mesh, rc) ! Check CICE mesh use ice_constants, only : c1,c0,c360 - use ice_grid , only : tlon, tlat + use ice_grid , only : tlon, tlat, hm ! input/output parameters type(ESMF_GridComp) , intent(inout) :: gcomp @@ -569,7 +569,8 @@ subroutine ice_mesh_check(gcomp, ice_mesh, rc) ! local variables type(ESMF_DistGrid) :: distGrid - integer :: n,c,g,i,j,m ! indices + type(ESMF_Array) :: elemMaskArray + integer :: n,i,j ! indices integer :: iblk, jblk ! indices integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain type(block) :: this_block ! block information for current block @@ -578,11 +579,15 @@ subroutine ice_mesh_check(gcomp, ice_mesh, rc) real(dbl_kind), pointer :: ownedElemCoords(:) real(dbl_kind), pointer :: lat(:), latMesh(:) real(dbl_kind), pointer :: lon(:), lonMesh(:) + integer , pointer :: model_mask(:) real(dbl_kind) :: diff_lon real(dbl_kind) :: diff_lat real(dbl_kind) :: rad_to_deg real(dbl_kind) :: tmplon, eps_imesh logical :: isPresent, isSet + logical :: mask_error + integer :: mask_internal + integer :: mask_file character(len=char_len_long) :: cvalue character(len=char_len_long) :: logmsg character(len=*), parameter :: subname = ' ice_mesh_check: ' @@ -606,7 +611,7 @@ subroutine ice_mesh_check(gcomp, ice_mesh, rc) allocate(ownedElemCoords(spatialDim*numownedelements)) allocate(lonmesh(numOwnedElements)) allocate(latmesh(numOwnedElements)) - call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords) + call ESMF_MeshGet(ice_mesh, ownedElemCoords=ownedElemCoords, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return do n = 1,numOwnedElements lonMesh(n) = ownedElemCoords(2*n-1) @@ -650,8 +655,45 @@ subroutine ice_mesh_check(gcomp, ice_mesh, rc) enddo enddo -100 format('ERROR: CICE n, lonmesh, lon, diff_lon = ',i6,2(f21.13,3x),d21.5) -101 format('ERROR: CICE n, latmesh, lat, diff_lat = ',i6,2(f21.13,3x),d21.5) + ! obtain internally generated ice mask for error checks + allocate(model_mask(numOwnedElements)) + call ESMF_MeshGet(ice_mesh, elementdistGrid=distGrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + elemMaskArray = ESMF_ArrayCreate(distGrid, model_mask, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_MeshGet(ice_mesh, elemMaskArray=elemMaskArray, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + mask_error = .false. + n=0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + do j = jlo, jhi + jlo = this_block%jlo + jhi = this_block%jhi + do i = ilo, ihi + ilo = this_block%ilo + ihi = this_block%ihi + n = n+1 + mask_internal = nint(hm(i,j,iblk),kind=dbl_kind) + mask_file = model_mask(n) + if (mask_internal /= mask_file) then + write(6,102) n,mask_internal,mask_file + mask_error = .true. + end if + enddo !i + enddo !j + enddo !iblk + if (mask_error) then + call abort_ice(error_message=subname, file=__FILE__, line=__LINE__) + end if + + call ESMF_ArrayDestroy(elemMaskArray, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + +100 format('ERROR: CICE n, mesh_lon , lon, diff_lon = ',i8,2(f21.13,3x),d21.5) +101 format('ERROR: CICE n, mesh_lat , lat, diff_lat = ',i8,2(f21.13,3x),d21.5) +102 format('ERROR: CICE n, mesh_internal, mask_file = ',i8,2(i2,2x)) ! deallocate memory deallocate(ownedElemCoords) From b76342343e3dfdc8ab73ff28f66932734257e07e Mon Sep 17 00:00:00 2001 From: David Bailey Date: Tue, 8 Nov 2022 14:03:45 -0700 Subject: [PATCH 2/3] Rearrange some wave code --- .../drivers/nuopc/cmeps/ice_import_export.F90 | 93 ++++++++++--------- 1 file changed, 47 insertions(+), 46 deletions(-) diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index d1408d528..1b9f340a0 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -181,13 +181,6 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam if (flds_wiso) then call fldlist_add(fldsToIce_num, fldsToIce, 'So_roce_wiso', ungridded_lbound=1, ungridded_ubound=3) end if - ! the following are advertised but might not be connected if they are not advertised in the - ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific - ! from wave - if (flds_wave) then - call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, & - ungridded_ubound=25) - end if ! from atmosphere call fldlist_add(fldsToIce_num, fldsToIce, 'inst_height_lowest' ) @@ -215,6 +208,14 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam ! from atm - dry dust deposition fluxes (4 sizes) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_dstdry', ungridded_lbound=1, ungridded_ubound=4) + ! the following are advertised but might not be connected if they are not advertised in the + ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific + ! from wave + if (flds_wave) then + call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, & + ungridded_ubound=25) + end if + do n = 1,fldsToIce_num call NUOPC_Advertise(importState, standardName=fldsToIce(n)%stdname, & TransferOfferGeomObject='will provide', rc=rc) @@ -554,29 +555,6 @@ subroutine ice_import( importState, rc ) areacor=med2mod_areacor, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25) - if (State_FldChk(importState, 'Sw_elevation_spectrum')) then - if (nfreq /= 25) then - call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ") - end if - call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - do k = 1,nfreq - n = 0 - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo; ihi = this_block%ihi - jlo = this_block%jlo; jhi = this_block%jhi - do j = jlo, jhi - do i = ilo, ihi - n = n+1 - wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n) - end do - end do - end do - end do - end if - ! perform a halo update if (.not.prescribed_ice) then @@ -610,6 +588,29 @@ subroutine ice_import( importState, rc ) end do !$OMP END PARALLEL DO + ! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25) + if (State_FldChk(importState, 'Sw_elevation_spectrum')) then + if (nfreq /= 25) then + call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ") + end if + call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + do k = 1,nfreq + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n) + end do + end do + end do + end do + end if + if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks @@ -1141,6 +1142,22 @@ subroutine ice_export( exportState, rc ) call state_setexport(exportState, 'Si_snowh' , input=tempfld , lmask=tmask, ifrac=ailohi, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ------ + ! optional floe diameter and ice thickness to wave + ! ------ + + ! Sea ice thickness (m) + if (State_FldChk(exportState, 'Si_thick')) then + call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Sea ice floe diameter (m) + if (State_FldChk(exportState, 'Si_floediam')) then + call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + ! ------ ! ice/atm fluxes computed by ice ! ------ @@ -1328,22 +1345,6 @@ subroutine ice_export( exportState, rc ) end do end if - ! ------ - ! optional floe diameter and ice thickness to wave - ! ------ - - ! Sea ice thickness (m) - if (State_FldChk(exportState, 'Si_thick')) then - call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - - ! Sea ice floe diameter (m) - if (State_FldChk(exportState, 'Si_floediam')) then - call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - end if - end subroutine ice_export !=============================================================================== From 3edf6ea8e2dad17ac5d31f458b6bc9125eda81f7 Mon Sep 17 00:00:00 2001 From: David Bailey Date: Thu, 10 Nov 2022 12:21:02 -0700 Subject: [PATCH 3/3] Add updates for UFS --- cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 6 ++++++ cicecore/drivers/nuopc/cmeps/ice_import_export.F90 | 1 - 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 4d8210232..afdee5590 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -641,6 +641,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call cice_init1 call t_stopf ('cice_init1') + !----------------------------------------------------------------- + ! Advertise fields + !----------------------------------------------------------------- + call ice_advertise_fields(gcomp, importState, exportState, flds_scalar_name, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + #endif !---------------------------------------------------------------------------- diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index 1b9f340a0..e4db010de 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -436,7 +436,6 @@ subroutine ice_import( importState, rc ) integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain type(block) :: this_block ! block information for current block real (kind=dbl_kind),allocatable :: aflds(:,:,:,:) - real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: work real (kind=dbl_kind) :: workx, worky real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP real (kind=dbl_kind) :: Tffresh