Skip to content

Commit

Permalink
Isotopes for CICE (#423)
Browse files Browse the repository at this point in the history 
* initial changes to cice driver for isotopes

* add isotope fluxes

* forcing and restarts for isotopes

* interim commit for isotopes

* more isotopes code including ice_step

* diagnostics and history for isotopes

* update both RunMods

* isotope testing option

* fix diags and ice_in; add isotope tests

* changing badger acct

* fix isotope indexing

* initialize trcrn(isotopes)

* fix isotope restarts

* Add isotope documentation

* updating copyright to 2020

* update icepack version

* update non standalone drivers for isotope changes

Co-authored-by: apcraig <[email protected]>
Co-authored-by: David Bailey <[email protected]>
  • Loading branch information
@eclare108213 @apcraig @dabail10
3 people authored Apr 8, 2020
1 parent 9ac1863 commit 1bf6e8d
Show file tree
Hide file tree
Showing 36 changed files with 858 additions and 1,584 deletions.
Binary file modified LICENSE.pdf
View file
Binary file not shown.
124 changes: 106 additions & 18 deletions cicecore/cicedynB/analysis/ice_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ module ice_diagnostics
use ice_fileunits, only: flush_fileunit
use ice_exit, only: abort_ice
use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
use icepack_intfc, only: icepack_max_aero
use icepack_intfc, only: icepack_max_aero, icepack_max_iso
use icepack_intfc, only: icepack_query_parameters
use icepack_intfc, only: icepack_query_tracer_flags
use icepack_intfc, only: icepack_query_tracer_indices
Expand Down Expand Up @@ -79,6 +79,10 @@ module ice_diagnostics
toten , & ! total ice/snow energy (J)
totes ! total ice/snow energy (J)

real (kind=dbl_kind), dimension(icepack_max_iso) :: &
totison , & ! total isotope mass
totisos ! total isotope mass

real (kind=dbl_kind), dimension(icepack_max_aero) :: &
totaeron , & ! total aerosol mass
totaeros ! total aerosol mass
Expand All @@ -89,8 +93,8 @@ module ice_diagnostics
integer (kind=int_kind), parameter, public :: &
check_step = 999999999, & ! begin printing at istep1=check_step
iblkp = 1, & ! block number
ip = 2, & ! i index
jp = 11, & ! j index
ip = 72, & ! i index
jp = 11, & ! j index
mtask = 0 ! my_task

!=======================================================================
Expand All @@ -113,15 +117,15 @@ subroutine runtime_diags (dt)
use ice_constants, only: c1, c1000, c2, p001, p5, &
field_loc_center, m2_to_km2
use ice_domain, only: distrb_info, nblocks
use ice_domain_size, only: ncat, n_aero, max_blocks, nfsd
use ice_domain_size, only: ncat, n_iso, n_aero, max_blocks, nfsd
use ice_flux, only: alvdr, alidr, alvdf, alidf, evap, fsnow, frazil, &
fswabs, fswthru, flw, flwout, fsens, fsurf, flat, frzmlt_init, frain, fpond, &
fhocn_ai, fsalt_ai, fresh_ai, frazil_diag, &
update_ocn_f, Tair, Qa, fsw, fcondtop, meltt, meltb, meltl, snoice, &
dsnow, congel, sst, sss, Tf, fhocn, &
swvdr, swvdf, swidr, swidf, &
alvdr_init, alvdf_init, alidr_init, alidf_init
use ice_flux_bgc, only: faero_atm, faero_ocn
use ice_flux_bgc, only: faero_atm, faero_ocn, fiso_atm, fiso_ocn
use ice_global_reductions, only: global_sum, global_sum_prod, global_maxval
use ice_grid, only: lmask_n, lmask_s, tarean, tareas
use ice_state ! everything
Expand All @@ -138,10 +142,11 @@ subroutine runtime_diags (dt)
integer (kind=int_kind) :: &
i, j, k, n, iblk, nc, &
ktherm, &
nt_tsfc, nt_aero, nt_fbri, nt_apnd, nt_hpnd, nt_fsd
nt_tsfc, nt_aero, nt_fbri, nt_apnd, nt_hpnd, nt_fsd, &
nt_isosno, nt_isoice

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

real (kind=dbl_kind) :: &
rhow, rhos, rhoi, puny, awtvdr, awtidr, awtvdf, awtidf, &
Expand All @@ -166,6 +171,13 @@ subroutine runtime_diags (dt)
delein, werrn, herrn, msltn, delmsltn, serrn, &
deleis, werrs, herrs, mslts, delmslts, serrs

! isotope diagnostics
real (kind=dbl_kind), dimension(icepack_max_aero) :: &
fisoan, fisoon, isorn, &
fisoas, fisoos, isors, &
isomx1n, isomx1s, &
isototn, isotots

! aerosol diagnostics
real (kind=dbl_kind), dimension(icepack_max_aero) :: &
faeran, faeron, aerrn, &
Expand All @@ -188,10 +200,10 @@ 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_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
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_fsd_out=nt_fsd,nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
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, &
Expand Down Expand Up @@ -683,6 +695,45 @@ subroutine runtime_diags (dt)
serrn = (sfsaltn + delmsltn) / (msltn + c1)
serrs = (sfsalts + delmslts) / (mslts + c1)

! isotopes
if (tr_iso) then
do n = 1, n_iso
fisoan(n) = global_sum_prod(fiso_atm(:,:,n,:), aice_init, &
distrb_info, field_loc_center, tarean)
fisoas(n) = global_sum_prod(fiso_atm(:,:,n,:), aice_init, &
distrb_info, field_loc_center, tareas)
fisoan(n) = fisoan(n)*dt
fisoas(n) = fisoas(n)*dt
fisoon(n) = global_sum_prod(fiso_ocn(:,:,n,:), aice, &
distrb_info, field_loc_center, tarean)
fisoos(n) = global_sum_prod(fiso_ocn(:,:,n,:), aice, &
distrb_info, field_loc_center, tareas)
fisoon(n) = fisoon(n)*dt
fisoos(n) = fisoos(n)*dt

!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = c0
do k = 1, n_iso
work1(i,j,iblk) = work1(i,j,iblk) &
+ vsno(i,j,iblk)*trcr(i,j,nt_isosno+k-1,iblk) &
+ vice(i,j,iblk)*trcr(i,j,nt_isoice+k-1,iblk)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
isototn(n) = global_sum(work1, distrb_info, field_loc_center, tarean)
isotots(n) = global_sum(work1, distrb_info, field_loc_center, tareas)
isomx1n(n) = global_maxval(work1, distrb_info, lmask_n)
isomx1s(n) = global_maxval(work1, distrb_info, lmask_s)
isorn(n) = (totison(n)-isototn(n)+fisoan(n)-fisoon(n))/(isototn(n)+c1)
isors(n) = (totisos(n)-isotots(n)+fisoas(n)-fisoos(n))/(isotots(n)+c1)
enddo ! n_iso
endif ! tr_iso

! aerosols
if (tr_aero) then
do n = 1, n_aero
Expand Down Expand Up @@ -917,6 +968,17 @@ subroutine runtime_diags (dt)
write(nu_diag,901) 'arwt salt flx error = ',serrn,serrs

write(nu_diag,*) '----------------------------'
if (tr_iso) then
do n = 1, n_iso
write(nu_diag,*) ' isotope ',n
write(nu_diag,901) 'fiso_atm (kg/m2) = ', fisoan(n), fisoas(n)
write(nu_diag,901) 'fiso_ocn (kg/m2) = ', fisoon(n), fisoos(n)
write(nu_diag,901) 'total iso (kg/m2) = ', isototn(n), isotots(n)
write(nu_diag,901) 'iso error = ', isorn(n), isors(n)
write(nu_diag,901) 'maximum iso (kg/m2) = ', isomx1n(n),isomx1s(n)
enddo
write(nu_diag,*) '----------------------------'
endif ! tr_iso
if (tr_aero) then
do n = 1, n_aero
write(nu_diag,*) ' aerosol ',n
Expand Down Expand Up @@ -1030,16 +1092,16 @@ subroutine init_mass_diags
use ice_blocks, only: nx_block, ny_block
use ice_constants, only: field_loc_center
use ice_domain, only: distrb_info, nblocks
use ice_domain_size, only: n_aero, ncat, max_blocks
use ice_domain_size, only: n_iso, n_aero, ncat, max_blocks
use ice_global_reductions, only: global_sum
use ice_grid, only: tareas, tarean
use ice_state, only: aicen, vice, vsno, trcrn, trcr

integer (kind=int_kind) :: n, i, j, iblk, &
nt_hpnd, nt_apnd, nt_aero
integer (kind=int_kind) :: n, i, j, k, iblk, &
nt_hpnd, nt_apnd, nt_aero, nt_isosno, nt_isoice

logical (kind=log_kind) :: &
tr_aero, tr_pond_topo
tr_iso, tr_aero, tr_pond_topo

real (kind=dbl_kind) :: &
shmaxn, snwmxn, shmaxs, snwmxs, totpn, totps, &
Expand All @@ -1051,7 +1113,8 @@ subroutine init_mass_diags
character(len=*), parameter :: subname = '(init_mass_diags)'

call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_pond_topo_out=tr_pond_topo)
call icepack_query_tracer_indices( &
call icepack_query_tracer_flags(tr_iso_out=tr_iso)
call icepack_query_tracer_indices(nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
nt_hpnd_out=nt_hpnd, nt_apnd_out=nt_apnd, nt_aero_out=nt_aero)
call icepack_query_parameters( &
rhoi_out=rhoi, rhos_out=rhos, rhofresh_out=rhofresh)
Expand Down Expand Up @@ -1094,6 +1157,27 @@ subroutine init_mass_diags
enddo ! npnt
endif ! print_points

if (tr_iso) then
do n=1,n_iso
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = c0
do k = 1, n_iso
work1(i,j,iblk) = work1(i,j,iblk) &
+ vsno(i,j,iblk)*trcr(i,j,nt_isosno+k-1,iblk) &
+ vice(i,j,iblk)*trcr(i,j,nt_isoice+k-1,iblk)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
totison(n)= global_sum(work1, distrb_info, field_loc_center, tarean)
totisos(n)= global_sum(work1, distrb_info, field_loc_center, tareas)
enddo
endif

if (tr_aero) then
do n=1,n_aero
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
Expand Down Expand Up @@ -1480,18 +1564,20 @@ subroutine print_state(plabel,i,j,iblk)
qi, qs, Tsnow, &
rad_to_deg, puny, rhoi, lfresh, rhos, cp_ice

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

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

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)
call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso)
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_qsno_out=nt_qsno, nt_fsd_out=nt_fsd, &
nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
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 @@ -1521,6 +1607,8 @@ 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
write(nu_diag,*) ' '

! dynamics (transport and/or ridging) causes the floe size distribution to become non-normal
Expand Down
Loading

0 comments on commit 1bf6e8d

Please sign in to comment.