Skip to content

Commit

Permalink
Merge cgridDEV branch including C grid implementation and other fixes (…
Browse files Browse the repository at this point in the history 
…CICE-Consortium#715)

Merge cgridDEV branch to main

Co-authored-by: Jean-Francois Lemieux <[email protected]>
Co-authored-by: David A. Bailey <[email protected]>
Co-authored-by: Elizabeth Hunke <[email protected]>
Co-authored-by: daveh150 <[email protected]>
Co-authored-by: TRasmussen <[email protected]>
Co-authored-by: Philippe Blain <[email protected]>
  • Loading branch information
@apcraig @JFLemieux73
@dabail10 @eclare108213 @daveh150 @TillRasmussen @phil-blain
7 people authored May 10, 2022
1 parent 8d7314f commit 078aab4
Show file tree
Hide file tree
Showing 135 changed files with 14,907 additions and 5,105 deletions.
1 change: 1 addition & 0 deletions .github/workflows/test-cice.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ jobs:
run: |
sudo xcode-select -r
sudo xcode-select -s /Library/Developer/CommandLineTools
sudo ln -s /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/* /usr/local/include/
echo "xcrun --show-sdk-path: $(xcrun --show-sdk-path)"
echo "xcode-select -p: $(xcode-select -p)"
- name: system info
Expand Down
172 changes: 159 additions & 13 deletions cicecore/cicedynB/analysis/ice_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ subroutine runtime_diags (dt)
alvdr_init, alvdf_init, alidr_init, alidf_init
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_grid, only: lmask_n, lmask_s, tarean, tareas, grid_ice
use ice_state ! everything
! tcraig, this is likely to cause circular dependency because ice_prescribed_mod is high level routine
#ifdef CESMCOUPLED
Expand Down Expand Up @@ -201,6 +201,17 @@ subroutine runtime_diags (dt)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
work1, work2

real (kind=dbl_kind), parameter :: &
maxval_spval = -0.9_dbl_kind*HUGE(0.0_dbl_kind) ! spval to detect
! undefined values returned from global_maxval. if global_maxval
! is applied to a region that does not exist (for instance
! southern hemisphere in box cases), global_maxval
! returns -HUGE which we want to avoid writing. The
! return value is checked against maxval_spval before writing.

! real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
! uvelT, vvelT

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

call icepack_query_parameters(ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc)
Expand Down Expand Up @@ -228,6 +239,8 @@ subroutine runtime_diags (dt)
! hemispheric quantities

! total ice area
arean = c0
areas = c0
arean = global_sum(aice, distrb_info, field_loc_center, tarean)
areas = global_sum(aice, distrb_info, field_loc_center, tareas)
arean = arean * m2_to_km2
Expand All @@ -244,6 +257,8 @@ subroutine runtime_diags (dt)
enddo
enddo
!$OMP END PARALLEL DO
extentn = c0
extents = c0
extentn = global_sum(work1, distrb_info, field_loc_center, &
tarean)
extents = global_sum(work1, distrb_info, field_loc_center, &
Expand All @@ -252,10 +267,14 @@ subroutine runtime_diags (dt)
extents = extents * m2_to_km2

! total ice volume
shmaxn = c0
shmaxs = c0
shmaxn = global_sum(vice, distrb_info, field_loc_center, tarean)
shmaxs = global_sum(vice, distrb_info, field_loc_center, tareas)

! total snow volume
snwmxn = c0
snwmxs = c0
snwmxn = global_sum(vsno, distrb_info, field_loc_center, tarean)
snwmxs = global_sum(vsno, distrb_info, field_loc_center, tareas)

Expand Down Expand Up @@ -293,7 +312,25 @@ subroutine runtime_diags (dt)
enddo
enddo
enddo
!$OMP END PARALLEL DO
! Eventually do energy diagnostic on T points.
! if (grid_ice == 'CD') then
! !$OMP PARALLEL DO PRIVATE(iblk,i,j)
! do iblk = 1, nblocks
! do j = 1, ny_block
! do i = 1, nx_block
! call grid_average_X2Y('E2TS',uvelE,uvelT)
! call grid_average_X2Y('N2TS',vvelN,vvelT)
! work1(i,j,iblk) = p5 &
! * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &
! * (uvelT(i,j,iblk)*uvelT(i,j,iblk) &
! + vvelT(i,j,iblk)*vvelT(i,j,iblk))
! enddo
! enddo
! enddo
! endif
! !$OMP END PARALLEL DO
ketotn = c0
ketots = c0
ketotn = global_sum(work1, distrb_info, field_loc_center, tarean)
ketots = global_sum(work1, distrb_info, field_loc_center, tareas)

Expand Down Expand Up @@ -370,23 +407,57 @@ subroutine runtime_diags (dt)
endif

! maximum ice volume (= mean thickness including open water)
hmaxn = c0
hmaxs = c0
hmaxn = global_maxval(vice, distrb_info, lmask_n)
hmaxs = global_maxval(vice, distrb_info, lmask_s)
if (hmaxn < maxval_spval) hmaxn = c0
if (hmaxs < maxval_spval) hmaxs = c0

! maximum ice speed
!$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) = sqrt(uvel(i,j,iblk)**2 &
+ vvel(i,j,iblk)**2)
if (grid_ice == 'CD') then
!$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) = max(sqrt(uvelE(i,j,iblk)**2 &
+ vvelE(i,j,iblk)**2), &
sqrt(uvelN(i,j,iblk)**2 &
+ vvelN(i,j,iblk)**2))
enddo
enddo
enddo
!$OMP END PARALLEL DO
elseif (grid_ice == 'C') then
!$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) = sqrt(uvelE(i,j,iblk)**2 &
+ vvelN(i,j,iblk)**2)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
else
!$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) = sqrt(uvel(i,j,iblk)**2 &
+ vvel(i,j,iblk)**2)
enddo
enddo
enddo
!$OMP END PARALLEL DO
endif

umaxn = c0
umaxs = c0
umaxn = global_maxval(work1, distrb_info, lmask_n)
umaxs = global_maxval(work1, distrb_info, lmask_s)
if (umaxn < maxval_spval) umaxn = c0
if (umaxs < maxval_spval) umaxs = c0

! Write warning message if ice speed is too big
! (Ice speeds of ~1 m/s or more usually indicate instability)
Expand Down Expand Up @@ -427,8 +498,12 @@ subroutine runtime_diags (dt)

! maximum ice strength

pmaxn = c0
pmaxs = c0
pmaxn = global_maxval(strength, distrb_info, lmask_n)
pmaxs = global_maxval(strength, distrb_info, lmask_s)
if (pmaxn < maxval_spval) pmaxn = c0
if (pmaxs < maxval_spval) pmaxs = c0

pmaxn = pmaxn / c1000 ! convert to kN/m
pmaxs = pmaxs / c1000
Expand All @@ -437,7 +512,9 @@ subroutine runtime_diags (dt)

! total ice/snow internal energy
call total_energy (work1)


etotn = c0
etots = c0
etotn = global_sum(work1, distrb_info, &
field_loc_center, tarean)
etots = global_sum(work1, distrb_info, &
Expand All @@ -452,6 +529,8 @@ subroutine runtime_diags (dt)

! evaporation

evpn = c0
evps = c0
evpn = global_sum_prod(evap, aice, distrb_info, &
field_loc_center, tarean)
evps = global_sum_prod(evap, aice, distrb_info, &
Expand All @@ -470,6 +549,8 @@ subroutine runtime_diags (dt)
endif

! salt flux
sfsaltn = c0
sfsalts = c0
sfsaltn = global_sum(fsalt_ai, distrb_info, &
field_loc_center, tarean)
sfsalts = global_sum(fsalt_ai, distrb_info, &
Expand All @@ -478,6 +559,8 @@ subroutine runtime_diags (dt)
sfsalts = sfsalts*dt

! fresh water flux
sfreshn = c0
sfreshs = c0
sfreshn = global_sum(fresh_ai, distrb_info, &
field_loc_center, tarean)
sfreshs = global_sum(fresh_ai, distrb_info, &
Expand All @@ -499,6 +582,8 @@ subroutine runtime_diags (dt)

! ocean heat
! Note: fswthru not included because it does not heat ice
fhocnn = c0
fhocns = c0
fhocnn = global_sum(fhocn_ai, distrb_info, &
field_loc_center, tarean)
fhocns = global_sum(fhocn_ai, distrb_info, &
Expand Down Expand Up @@ -548,6 +633,8 @@ subroutine runtime_diags (dt)

endif ! calc_Tsfc

fhatmn = c0
fhatms = c0
fhatmn = global_sum(work1, distrb_info, &
field_loc_center, tarean)
fhatms = global_sum(work1, distrb_info, &
Expand All @@ -564,6 +651,8 @@ subroutine runtime_diags (dt)
enddo
!$OMP END PARALLEL DO

fswnetn = c0
fswnets = c0
fswnetn = global_sum(work1, distrb_info, &
field_loc_center, tarean)
fswnets = global_sum(work1, distrb_info, &
Expand All @@ -582,6 +671,8 @@ subroutine runtime_diags (dt)
enddo
!$OMP END PARALLEL DO

fswdnn = c0
fswdns = c0
fswdnn = global_sum(work1, distrb_info, &
field_loc_center, tarean)
fswdns = global_sum(work1, distrb_info, &
Expand All @@ -597,12 +688,17 @@ subroutine runtime_diags (dt)
enddo
enddo
!$OMP END PARALLEL DO

fhfrzn = c0
fhfrzs = c0
fhfrzn = global_sum(work1, distrb_info, &
field_loc_center, tarean)
fhfrzs = global_sum(work1, distrb_info, &
field_loc_center, tareas)

! rain
rnn = c0
rns = c0
rnn = global_sum_prod(frain, aice_init, distrb_info, &
field_loc_center, tarean)
rns = global_sum_prod(frain, aice_init, distrb_info, &
Expand All @@ -611,6 +707,8 @@ subroutine runtime_diags (dt)
rns = rns*dt

! snow
snn = c0
sns = c0
snn = global_sum_prod(fsnow, aice_init, distrb_info, &
field_loc_center, tarean)
sns = global_sum_prod(fsnow, aice_init, distrb_info, &
Expand All @@ -623,6 +721,8 @@ subroutine runtime_diags (dt)
work1(:,:,:) = frazil(:,:,:)*rhoi/dt
if (ktherm == 2 .and. .not.update_ocn_f) &
work1(:,:,:) = (frazil(:,:,:)-frazil_diag(:,:,:))*rhoi/dt
frzn = c0
frzs = c0
frzn = global_sum(work1, distrb_info, &
field_loc_center, tarean)
frzs = global_sum(work1, distrb_info, &
Expand Down Expand Up @@ -706,6 +806,16 @@ subroutine runtime_diags (dt)

! isotopes
if (tr_iso) then
fisoan = c0
fisoas = c0
fisoon = c0
fisoos = c0
isototn = c0
isotots = c0
isomx1n = c0
isomx1s = c0
isorn = c0
isors = c0
do n = 1, n_iso
fisoan(n) = global_sum_prod(fiso_atm(:,:,n,:), aice_init, &
distrb_info, field_loc_center, tarean)
Expand Down Expand Up @@ -738,13 +848,25 @@ subroutine runtime_diags (dt)
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)
if (isomx1n(n) < maxval_spval) isomx1n(n) = c0
if (isomx1s(n) < maxval_spval) isomx1s(n) = c0
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
faeran = c0
faeras = c0
faeron = c0
faeros = c0
aerototn = c0
aerotots = c0
aeromx1n = c0
aeromx1s = c0
aerrn = c0
aerrs = c0
do n = 1, n_aero
faeran(n) = global_sum_prod(faero_atm(:,:,n,:), aice_init, &
distrb_info, field_loc_center, tarean)
Expand Down Expand Up @@ -776,6 +898,8 @@ subroutine runtime_diags (dt)
aerotots(n) = global_sum(work1, distrb_info, field_loc_center, tareas)
aeromx1n(n) = global_maxval(work1, distrb_info, lmask_n)
aeromx1s(n) = global_maxval(work1, distrb_info, lmask_s)
if (aeromx1n(n) < maxval_spval) aeromx1n(n) = c0
if (aeromx1s(n) < maxval_spval) aeromx1s(n) = c0

aerrn(n) = (totaeron(n)-aerototn(n)+faeran(n)-faeron(n)) &
/ (aerototn(n) + c1)
Expand Down Expand Up @@ -1629,10 +1753,12 @@ end subroutine debug_ice

subroutine print_state(plabel,i,j,iblk)

use ice_grid, only: grid_ice
use ice_blocks, only: block, get_block
use ice_domain, only: blocks_ice
use ice_domain_size, only: ncat, nilyr, nslyr, nfsd
use ice_state, only: aice0, aicen, vicen, vsnon, uvel, vvel, trcrn
use ice_state, only: aice0, aicen, vicen, vsnon, uvel, vvel, &
uvelE, vvelE, uvelN, vvelN, trcrn
use ice_flux, only: uatm, vatm, potT, Tair, Qa, flw, frain, fsnow, &
fsens, flat, evap, flwout, swvdr, swvdf, swidr, swidf, rhoa, &
frzmlt, sst, sss, Tf, Tref, Qref, Uref, uocn, vocn, strtltx, strtlty
Expand Down Expand Up @@ -1754,6 +1880,15 @@ subroutine print_state(plabel,i,j,iblk)

write(nu_diag,*) 'uvel(i,j)',uvel(i,j,iblk)
write(nu_diag,*) 'vvel(i,j)',vvel(i,j,iblk)
if (grid_ice == 'C') then
write(nu_diag,*) 'uvelE(i,j)',uvelE(i,j,iblk)
write(nu_diag,*) 'uvelN(i,j)',uvelN(i,j,iblk)
elseif (grid_ice == 'CD') then
write(nu_diag,*) 'uvelE(i,j)',uvelE(i,j,iblk)
write(nu_diag,*) 'vvelE(i,j)',vvelE(i,j,iblk)
write(nu_diag,*) 'uvelN(i,j)',uvelN(i,j,iblk)
write(nu_diag,*) 'vvelN(i,j)',vvelN(i,j,iblk)
endif

write(nu_diag,*) ' '
write(nu_diag,*) 'atm states and fluxes'
Expand Down Expand Up @@ -1801,10 +1936,12 @@ end subroutine print_state

subroutine print_points_state(plabel,ilabel)

use ice_grid, only: grid_ice
use ice_blocks, only: block, get_block
use ice_domain, only: blocks_ice
use ice_domain_size, only: ncat, nilyr, nslyr
use ice_state, only: aice0, aicen, vicen, vsnon, uvel, vvel, trcrn
use ice_state, only: aice0, aicen, vicen, vsnon, uvel, vvel, &
uvelE, vvelE, uvelE, vvelE, trcrn
use ice_flux, only: uatm, vatm, potT, Tair, Qa, flw, frain, fsnow, &
fsens, flat, evap, flwout, swvdr, swvdf, swidr, swidf, rhoa, &
frzmlt, sst, sss, Tf, Tref, Qref, Uref, uocn, vocn, strtltx, strtlty
Expand Down Expand Up @@ -1896,6 +2033,15 @@ subroutine print_points_state(plabel,ilabel)

write(nu_diag,*) trim(llabel),'uvel=',uvel(i,j,iblk)
write(nu_diag,*) trim(llabel),'vvel=',vvel(i,j,iblk)
if (grid_ice == 'C') then
write(nu_diag,*) trim(llabel),'uvelE=',uvelE(i,j,iblk)
write(nu_diag,*) trim(llabel),'vvelN=',vvelN(i,j,iblk)
elseif (grid_ice == 'CD') then
write(nu_diag,*) trim(llabel),'uvelE=',uvelE(i,j,iblk)
write(nu_diag,*) trim(llabel),'vvelE=',vvelE(i,j,iblk)
write(nu_diag,*) trim(llabel),'uvelN=',uvelN(i,j,iblk)
write(nu_diag,*) trim(llabel),'vvelN=',vvelN(i,j,iblk)
endif

write(nu_diag,*) ' '
write(nu_diag,*) 'atm states and fluxes'
Expand Down
Loading

0 comments on commit 078aab4

Please sign in to comment.