Skip to content

Commit

Permalink
Update CICE for latest Consortium master (NOAA-EMC#38)
Browse files Browse the repository at this point in the history 
    * Implement advanced snow physics in icepack and CICE
    * Fix time-stamping of CICE history files
    * Fix CICE history file precision
  • Loading branch information
@DeniseWorthen
DeniseWorthen committed May 10, 2024
1 parent bf1dfc5 commit 955ef8d
Show file tree
Hide file tree
Showing 55 changed files with 5,638 additions and 3,637 deletions.
143 changes: 115 additions & 28 deletions cicecore/cicedyn/analysis/ice_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ module ice_diagnostics
use ice_communicate, only: my_task, master_task
use ice_constants, only: c0, c1
use ice_calendar, only: istep1
use ice_domain_size, only: nslyr
use ice_fileunits, only: nu_diag
use ice_fileunits, only: flush_fileunit
use ice_exit, only: abort_ice
Expand All @@ -39,7 +40,7 @@ module ice_diagnostics
print_global ! if true, print global data

integer (kind=int_kind), public :: &
debug_model_step = 999999999 ! begin printing at istep1=debug_model_step
debug_model_step = 0 ! begin printing at istep1=debug_model_step

integer (kind=int_kind), parameter, public :: &
npnt = 2 ! total number of points to be printed
Expand Down Expand Up @@ -73,6 +74,12 @@ module ice_diagnostics
integer (kind=int_kind), dimension(npnt), public :: &
piloc, pjloc, pbloc, pmloc ! location of diagnostic points

integer (kind=int_kind), public :: &
debug_model_i = -1, & ! location of debug_model point, local i index
debug_model_j = -1, & ! location of debug_model point, local j index
debug_model_iblk = -1, & ! location of debug_model point, local block number
debug_model_task = -1 ! location of debug_model point, local task number

! for hemispheric water and heat budgets
real (kind=dbl_kind) :: &
totmn , & ! total ice/snow water mass (nh)
Expand Down Expand Up @@ -136,15 +143,19 @@ subroutine runtime_diags (dt)
i, j, k, n, iblk, nc, &
ktherm, &
nt_tsfc, nt_aero, nt_fbri, nt_apnd, nt_hpnd, nt_fsd, &
nt_isosno, nt_isoice
nt_isosno, nt_isoice, nt_rsnw, nt_rhos, nt_smice, nt_smliq

logical (kind=log_kind) :: &
tr_pond_topo, tr_brine, tr_iso, tr_aero, calc_Tsfc, tr_fsd
tr_pond_topo, tr_brine, tr_iso, tr_aero, calc_Tsfc, tr_fsd, &
tr_snow, snwgrain

real (kind=dbl_kind) :: &
rhow, rhos, rhoi, puny, awtvdr, awtidr, awtvdf, awtidf, &
rhofresh, lfresh, lvap, ice_ref_salinity, Tffresh

character (len=char_len) :: &
snwredist

! hemispheric state quantities
real (kind=dbl_kind) :: &
umaxn, hmaxn, shmaxn, arean, snwmxn, extentn, shmaxnt, &
Expand Down Expand Up @@ -184,7 +195,8 @@ subroutine runtime_diags (dt)
pTsfc, pevap, pfswabs, pflwout, pflat, pfsens, &
pfsurf, pfcondtop, psst, psss, pTf, hiavg, hsavg, hbravg, &
pfhocn, psalt, fsdavg, &
pmeltt, pmeltb, pmeltl, psnoice, pdsnow, pfrazil, pcongel
pmeltt, pmeltb, pmeltl, psnoice, pdsnow, pfrazil, pcongel, &
prsnwavg, prhosavg, psmicetot, psmliqtot, psmtot

real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
work1, work2
Expand All @@ -193,15 +205,19 @@ subroutine runtime_diags (dt)

call icepack_query_parameters(ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc)
call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_aero_out=tr_aero, &
tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso, &
tr_snow_out=tr_snow)
call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, nt_Tsfc_out=nt_Tsfc, &
nt_aero_out=nt_aero, nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, &
nt_fsd_out=nt_fsd,nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
nt_fsd_out=nt_fsd,nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
nt_rsnw_out=nt_rsnw, nt_rhos_out=nt_rhos, &
nt_smice_out=nt_smice, nt_smliq_out=nt_smliq)
call icepack_query_parameters(Tffresh_out=Tffresh, rhos_out=rhos, &
rhow_out=rhow, rhoi_out=rhoi, puny_out=puny, &
awtvdr_out=awtvdr, awtidr_out=awtidr, awtvdf_out=awtvdf, awtidf_out=awtidf, &
rhofresh_out=rhofresh, lfresh_out=lfresh, lvap_out=lvap, &
ice_ref_salinity_out=ice_ref_salinity)
ice_ref_salinity_out=ice_ref_salinity,snwredist_out=snwredist, &
snwgrain_out=snwgrain)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)
Expand Down Expand Up @@ -819,6 +835,27 @@ subroutine runtime_diags (dt)
enddo
endif
endif
if (tr_snow) then ! snow tracer quantities
prsnwavg (n) = c0 ! avg snow grain radius
prhosavg (n) = c0 ! avg snow density
psmicetot(n) = c0 ! total mass of ice in snow (kg/m2)
psmliqtot(n) = c0 ! total mass of liquid in snow (kg/m2)
psmtot (n) = c0 ! total mass of snow volume (kg/m2)
if (vsno(i,j,iblk) > c0) then
do k = 1, nslyr
prsnwavg (n) = prsnwavg (n) + trcr(i,j,nt_rsnw +k-1,iblk) ! snow grain radius
prhosavg (n) = prhosavg (n) + trcr(i,j,nt_rhos +k-1,iblk) ! compacted snow density
psmicetot(n) = psmicetot(n) + trcr(i,j,nt_smice+k-1,iblk) * vsno(i,j,iblk)
psmliqtot(n) = psmliqtot(n) + trcr(i,j,nt_smliq+k-1,iblk) * vsno(i,j,iblk)
end do
endif
psmtot (n) = rhos * vsno(i,j,iblk) ! mass of ice in standard density snow
prsnwavg (n) = prsnwavg (n) / real(nslyr,kind=dbl_kind) ! snow grain radius
prhosavg (n) = prhosavg (n) / real(nslyr,kind=dbl_kind) ! compacted snow density
psmicetot(n) = psmicetot(n) / real(nslyr,kind=dbl_kind) ! mass of ice in snow
psmliqtot(n) = psmliqtot(n) / real(nslyr,kind=dbl_kind) ! mass of liquid in snow
end if
psalt(n) = c0
if (vice(i,j,iblk) /= c0) psalt(n) = work2(i,j,iblk)/vice(i,j,iblk)
pTsfc(n) = trcr(i,j,nt_Tsfc,iblk) ! ice/snow sfc temperature
pevap(n) = evap(i,j,iblk)*dt/rhoi ! sublimation/condensation
Expand Down Expand Up @@ -870,6 +907,11 @@ subroutine runtime_diags (dt)
call broadcast_scalar(pmeltl (n), pmloc(n))
call broadcast_scalar(psnoice (n), pmloc(n))
call broadcast_scalar(pdsnow (n), pmloc(n))
call broadcast_scalar(psmtot (n), pmloc(n))
call broadcast_scalar(prsnwavg (n), pmloc(n))
call broadcast_scalar(prhosavg (n), pmloc(n))
call broadcast_scalar(psmicetot(n), pmloc(n))
call broadcast_scalar(psmliqtot(n), pmloc(n))
call broadcast_scalar(pfrazil (n), pmloc(n))
call broadcast_scalar(pcongel (n), pmloc(n))
call broadcast_scalar(pdhi (n), pmloc(n))
Expand Down Expand Up @@ -1053,6 +1095,26 @@ subroutine runtime_diags (dt)
write(nu_diag,900) 'effective dhi (m) = ',pdhi(1),pdhi(2)
write(nu_diag,900) 'effective dhs (m) = ',pdhs(1),pdhs(2)
write(nu_diag,900) 'intnl enrgy chng(W/m^2)= ',pde (1),pde (2)

if (tr_snow) then
if (trim(snwredist) /= 'none') then
write(nu_diag,900) 'avg snow density(kg/m3)= ',prhosavg(1) &
,prhosavg(2)
endif
if (snwgrain) then
write(nu_diag,900) 'avg snow grain radius = ',prsnwavg(1) &
,prsnwavg(2)
write(nu_diag,900) 'mass ice in snow(kg/m2)= ',psmicetot(1) &
,psmicetot(2)
write(nu_diag,900) 'mass liq in snow(kg/m2)= ',psmliqtot(1) &
,psmliqtot(2)
write(nu_diag,900) 'mass std snow (kg/m2)= ',psmtot(1) &
,psmtot(2)
write(nu_diag,900) 'max ice+liq (kg/m2)= ',rhow * hsavg(1) &
,rhow * hsavg(2)
endif
endif

write(nu_diag,*) '----------ocn----------'
write(nu_diag,900) 'sst (C) = ',psst(1),psst(2)
write(nu_diag,900) 'sss (ppt) = ',psss(1),psss(2)
Expand Down Expand Up @@ -1432,9 +1494,9 @@ subroutine init_diags
write(nu_diag,*) ' Find indices of diagnostic points '
endif

piloc(:) = 0
pjloc(:) = 0
pbloc(:) = 0
piloc(:) = -1
pjloc(:) = -1
pbloc(:) = -1
pmloc(:) = -999
plat(:) = -999._dbl_kind
plon(:) = -999._dbl_kind
Expand Down Expand Up @@ -1535,16 +1597,29 @@ subroutine debug_ice(iblk, plabeld)
integer (kind=int_kind) :: i, j, m
character(len=*), parameter :: subname='(debug_ice)'

! tcraig, do this only on one point, the first point
! do m = 1, npnt
m = 1
if (istep1 >= debug_model_step .and. &
iblk == pbloc(m) .and. my_task == pmloc(m)) then
i = piloc(m)
j = pjloc(m)
call print_state(plabeld,i,j,iblk)
if (istep1 >= debug_model_step) then

! set debug point to 1st global point if not set as local values
if (debug_model_i < 0 .and. debug_model_j < 0 .and. &
debug_model_iblk < 0 .and. debug_model_task < 0) then
debug_model_i = piloc(1)
debug_model_j = pjloc(1)
debug_model_task = pmloc(1)
debug_model_iblk = pbloc(1)
endif

! if debug point is messed up, abort
if (debug_model_i < 0 .or. debug_model_j < 0 .or. &
debug_model_iblk < 0 .or. debug_model_task < 0) then
call abort_ice (subname//'ERROR: debug_model_[i,j,iblk,mytask] not set correctly')
endif
! enddo

! write out debug info
if (debug_model_iblk == iblk .and. debug_model_task == my_task) then
call print_state(plabeld,debug_model_i,debug_model_j,debug_model_iblk)
endif

endif

end subroutine debug_ice

Expand Down Expand Up @@ -1573,23 +1648,25 @@ subroutine print_state(plabel,i,j,iblk)

real (kind=dbl_kind) :: &
eidebug, esdebug, &
qi, qs, Tsnow, &
qi, qs, Tsnow, si, &
rad_to_deg, puny, rhoi, lfresh, rhos, cp_ice

integer (kind=int_kind) :: n, k, nt_Tsfc, nt_qice, nt_qsno, nt_fsd, &
nt_isosno, nt_isoice
nt_isosno, nt_isoice, nt_sice, nt_smice, nt_smliq

logical (kind=log_kind) :: tr_fsd, tr_iso
logical (kind=log_kind) :: tr_fsd, tr_iso, tr_snow

type (block) :: &
this_block ! block information for current block

character(len=*), parameter :: subname = '(print_state)'

call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso, &
tr_snow_out=tr_snow)
call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, &
nt_qsno_out=nt_qsno, nt_fsd_out=nt_fsd, &
nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, &
nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
nt_smice_out=nt_smice, nt_smliq_out=nt_smliq)
call icepack_query_parameters( &
rad_to_deg_out=rad_to_deg, puny_out=puny, rhoi_out=rhoi, lfresh_out=lfresh, &
rhos_out=rhos, cp_ice_out=cp_ice)
Expand Down Expand Up @@ -1619,8 +1696,11 @@ subroutine print_state(plabel,i,j,iblk)
endif
write(nu_diag,*) 'Tsfcn',trcrn(i,j,nt_Tsfc,n,iblk)
if (tr_fsd) write(nu_diag,*) 'afsdn',trcrn(i,j,nt_fsd,n,iblk) ! fsd cat 1
! if (tr_iso) write(nu_diag,*) 'isosno',trcrn(i,j,nt_isosno,n,iblk) ! isotopes in snow
! if (tr_iso) write(nu_diag,*) 'isoice',trcrn(i,j,nt_isoice,n,iblk) ! isotopes in ice
! layer 1 diagnostics
! if (tr_iso) write(nu_diag,*) 'isosno',trcrn(i,j,nt_isosno,n,iblk) ! isotopes in snow
! if (tr_iso) write(nu_diag,*) 'isoice',trcrn(i,j,nt_isoice,n,iblk) ! isotopes in ice
! if (tr_snow) write(nu_diag,*) 'smice', trcrn(i,j,nt_smice, n,iblk) ! ice mass in snow
! if (tr_snow) write(nu_diag,*) 'smliq', trcrn(i,j,nt_smliq, n,iblk) ! liquid mass in snow
write(nu_diag,*) ' '

! dynamics (transport and/or ridging) causes the floe size distribution to become non-normal
Expand All @@ -1633,7 +1713,6 @@ subroutine print_state(plabel,i,j,iblk)

enddo ! n


eidebug = c0
do n = 1,ncat
do k = 1,nilyr
Expand Down Expand Up @@ -1666,6 +1745,14 @@ subroutine print_state(plabel,i,j,iblk)
write(nu_diag,*) 'qsnow(i,j)',esdebug
write(nu_diag,*) ' '

do n = 1,ncat
do k = 1,nilyr
si = trcrn(i,j,nt_sice+k-1,n,iblk)
write(nu_diag,*) 'sice, cat ',n,' layer ',k, si
enddo
enddo
write(nu_diag,*) ' '

write(nu_diag,*) 'uvel(i,j)',uvel(i,j,iblk)
write(nu_diag,*) 'vvel(i,j)',vvel(i,j,iblk)

Expand Down
Loading

0 comments on commit 955ef8d

Please sign in to comment.