From aa1e066ad182c92ad94c83065d7964ae5ce3657c Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Mon, 7 Nov 2022 11:25:53 -0800 Subject: [PATCH 1/2] Update KE diagnostic (#784) * Update KE diagnostic - Average U and V fields to T grid for KE calculation for B, C, and CD grid, raw U and V were combined with T grid variables in prior implementation. Also - Update max speed calculation for C grid, calc speed at both E and N points - Refactor highfreq option u,v passed into icepack for windstress calculation - This is bit-for-bit - Create local variables uvelT_icep and vvelT_icep to store highfreq info - Compute [uvelT,vvelT]_icep at start of therm calc instead of end of dynamics calc - No longer need to recalc [uvelT,vvelT]_icep on restart - Update some indentation in ice_diagnostics.F90 * Add grid average NE2TA, NE2UA --- .../cicedynB/analysis/ice_diagnostics.F90 | 156 +++++++++--------- cicecore/cicedynB/general/ice_state.F90 | 4 - cicecore/cicedynB/general/ice_step_mod.F90 | 57 ++++--- cicecore/cicedynB/infrastructure/ice_grid.F90 | 134 ++++++++++----- .../infrastructure/ice_restart_driver.F90 | 5 +- 5 files changed, 207 insertions(+), 149 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_diagnostics.F90 b/cicecore/cicedynB/analysis/ice_diagnostics.F90 index 83eb840d6..8879d6632 100644 --- a/cicecore/cicedynB/analysis/ice_diagnostics.F90 +++ b/cicecore/cicedynB/analysis/ice_diagnostics.F90 @@ -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, grid_ice + use ice_grid, only: lmask_n, lmask_s, tarean, tareas, grid_ice, grid_average_X2Y use ice_state ! everything ! tcraig, this is likely to cause circular dependency because ice_prescribed_mod is high level routine #ifdef CESMCOUPLED @@ -199,7 +199,8 @@ subroutine runtime_diags (dt) prsnwavg, prhosavg, psmicetot, psmliqtot, psmtot real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & - work1, work2 + uvelT, vvelT, & ! u,v on T points + work1, work2 ! temporary real (kind=dbl_kind), parameter :: & maxval_spval = -0.9_dbl_kind*HUGE(0.0_dbl_kind) ! spval to detect @@ -209,9 +210,6 @@ subroutine runtime_diags (dt) ! 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) @@ -284,51 +282,45 @@ subroutine runtime_diags (dt) if (tr_pond_topo) then !$OMP PARALLEL DO PRIVATE(iblk,i,j,n) do iblk = 1, nblocks - do j = 1, ny_block - do i = 1, nx_block - work1(i,j,iblk) = c0 - do n = 1, ncat - work1(i,j,iblk) = work1(i,j,iblk) & - + aicen(i,j,n,iblk) & - * trcrn(i,j,nt_apnd,n,iblk) & - * trcrn(i,j,nt_hpnd,n,iblk) + do j = 1, ny_block + do i = 1, nx_block + work1(i,j,iblk) = c0 + do n = 1, ncat + work1(i,j,iblk) = work1(i,j,iblk) & + + aicen(i,j,n,iblk) & + * trcrn(i,j,nt_apnd,n,iblk) & + * trcrn(i,j,nt_hpnd,n,iblk) + enddo + enddo enddo - enddo - enddo enddo !$OMP END PARALLEL DO ptotn = global_sum(work1, distrb_info, field_loc_center, tarean) ptots = global_sum(work1, distrb_info, field_loc_center, tareas) endif - ! total ice-snow kinetic energy + ! total ice-snow kinetic energy, on T points. + if (grid_ice == 'B') then + call grid_average_X2Y('A',uvel ,'U',uvelT,'T') + call grid_average_X2Y('A',vvel ,'U',vvelT,'T') + elseif (grid_ice == 'C') then + call grid_average_X2Y('A',uvelE,'E',uvelT,'T') + call grid_average_X2Y('A',vvelN,'N',vvelT,'T') + elseif (grid_ice == 'CD') then + call grid_average_X2Y('A',uvelE,'E',uvelN,'N',uvelT,'T') + call grid_average_X2Y('A',vvelE,'E',vvelN,'N',vvelT,'T') + endif + !$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) = p5 & - * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) & - * (uvel(i,j,iblk)**2 + vvel(i,j,iblk)**2) + work1(i,j,iblk) = p5 * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) & + * (uvelT(i,j,iblk)**2 + vvelT(i,j,iblk)**2) enddo enddo enddo - ! 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 + !$OMP END PARALLEL DO ketotn = c0 ketots = c0 ketotn = global_sum(work1, distrb_info, field_loc_center, tarean) @@ -420,21 +412,22 @@ subroutine runtime_diags (dt) 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)) + 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 + ! map uvelE to N and vvelN to E then compute max on E and N + call grid_average_X2Y('A',uvelE,'E',work1,'N') ! work1 =~ uvelN + call grid_average_X2Y('A',vvelN,'N',work2,'E') ! work2 =~ vvelE !$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) + work1(i,j,iblk) = max(sqrt(uvelE(i,j,iblk)**2 + work2(i,j,iblk)**2), & + sqrt(work1(i,j,iblk)**2 + vvelN(i,j,iblk)**2)) enddo enddo enddo @@ -444,8 +437,7 @@ subroutine runtime_diags (dt) 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) + work1(i,j,iblk) = sqrt(uvel(i,j,iblk)**2 + vvel(i,j,iblk)**2) enddo enddo enddo @@ -466,31 +458,31 @@ subroutine runtime_diags (dt) if (umaxn > umax_stab) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks - do j = 1, ny_block - do i = 1, nx_block - if (abs(work1(i,j,iblk) - umaxn) < puny) then - write(nu_diag,*) ' ' - write(nu_diag,*) 'Warning, large ice speed' - write(nu_diag,*) 'my_task, iblk, i, j, umaxn:', & - my_task, iblk, i, j, umaxn - endif - enddo - enddo + do j = 1, ny_block + do i = 1, nx_block + if (abs(work1(i,j,iblk) - umaxn) < puny) then + write(nu_diag,*) ' ' + write(nu_diag,*) 'Warning, large ice speed' + write(nu_diag,*) 'my_task, iblk, i, j, umaxn:', & + my_task, iblk, i, j, umaxn + endif + enddo + enddo enddo !$OMP END PARALLEL DO elseif (umaxs > umax_stab) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks - do j = 1, ny_block - do i = 1, nx_block - if (abs(work1(i,j,iblk) - umaxs) < puny) then - write(nu_diag,*) ' ' - write(nu_diag,*) 'Warning, large ice speed' - write(nu_diag,*) 'my_task, iblk, i, j, umaxs:', & - my_task, iblk, i, j, umaxs - endif - enddo - enddo + do j = 1, ny_block + do i = 1, nx_block + if (abs(work1(i,j,iblk) - umaxs) < puny) then + write(nu_diag,*) ' ' + write(nu_diag,*) 'Warning, large ice speed' + write(nu_diag,*) 'my_task, iblk, i, j, umaxs:', & + my_task, iblk, i, j, umaxs + endif + enddo + enddo enddo !$OMP END PARALLEL DO endif ! umax @@ -1357,14 +1349,14 @@ subroutine init_mass_diags do n=1,n_aero !$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) = trcr(i,j,nt_aero +4*(n-1),iblk)*vsno(i,j,iblk) & - + trcr(i,j,nt_aero+1+4*(n-1),iblk)*vsno(i,j,iblk) & - + trcr(i,j,nt_aero+2+4*(n-1),iblk)*vice(i,j,iblk) & - + trcr(i,j,nt_aero+3+4*(n-1),iblk)*vice(i,j,iblk) - enddo - enddo + do j = 1, ny_block + do i = 1, nx_block + work1(i,j,iblk) = trcr(i,j,nt_aero +4*(n-1),iblk)*vsno(i,j,iblk) & + + trcr(i,j,nt_aero+1+4*(n-1),iblk)*vsno(i,j,iblk) & + + trcr(i,j,nt_aero+2+4*(n-1),iblk)*vice(i,j,iblk) & + + trcr(i,j,nt_aero+3+4*(n-1),iblk)*vice(i,j,iblk) + enddo + enddo enddo !$OMP END PARALLEL DO totaeron(n)= global_sum(work1, distrb_info, field_loc_center, tarean) @@ -1377,17 +1369,17 @@ subroutine init_mass_diags totps = c0 !$OMP PARALLEL DO PRIVATE(iblk,i,j,n) do iblk = 1, nblocks - do j = 1, ny_block - do i = 1, nx_block - work1(i,j,iblk) = c0 - do n = 1, ncat - work1(i,j,iblk) = work1(i,j,iblk) & - + aicen(i,j,n,iblk) & - * trcrn(i,j,nt_apnd,n,iblk) & - * trcrn(i,j,nt_hpnd,n,iblk) + do j = 1, ny_block + do i = 1, nx_block + work1(i,j,iblk) = c0 + do n = 1, ncat + work1(i,j,iblk) = work1(i,j,iblk) & + + aicen(i,j,n,iblk) & + * trcrn(i,j,nt_apnd,n,iblk) & + * trcrn(i,j,nt_hpnd,n,iblk) + enddo + enddo enddo - enddo - enddo enddo !$OMP END PARALLEL DO totpn = global_sum(work1, distrb_info, field_loc_center, tarean) diff --git a/cicecore/cicedynB/general/ice_state.F90 b/cicecore/cicedynB/general/ice_state.F90 index 10e0aabf8..7b20879bc 100644 --- a/cicecore/cicedynB/general/ice_state.F90 +++ b/cicecore/cicedynB/general/ice_state.F90 @@ -110,8 +110,6 @@ module ice_state public :: & uvel , & ! x-component of velocity on U grid (m/s) vvel , & ! y-component of velocity on U grid (m/s) - uvelT , & ! x-component of velocity on T grid (m/s) - vvelT , & ! y-component of velocity on T grid (m/s) uvelE , & ! x-component of velocity on E grid (m/s) vvelE , & ! y-component of velocity on E grid (m/s) uvelN , & ! x-component of velocity on N grid (m/s) @@ -159,8 +157,6 @@ subroutine alloc_state aice0 (nx_block,ny_block,max_blocks) , & ! concentration of open water uvel (nx_block,ny_block,max_blocks) , & ! x-component of velocity on U grid (m/s) vvel (nx_block,ny_block,max_blocks) , & ! y-component of velocity on U grid (m/s) - uvelT (nx_block,ny_block,max_blocks) , & ! x-component of velocity on T grid (m/s) - vvelT (nx_block,ny_block,max_blocks) , & ! y-component of velocity on T grid (m/s) uvelE (nx_block,ny_block,max_blocks) , & ! x-component of velocity on E grid (m/s) vvelE (nx_block,ny_block,max_blocks) , & ! y-component of velocity on E grid (m/s) uvelN (nx_block,ny_block,max_blocks) , & ! x-component of velocity on N grid (m/s) diff --git a/cicecore/cicedynB/general/ice_step_mod.F90 b/cicecore/cicedynB/general/ice_step_mod.F90 index 39f10ffdf..5742bb2b9 100644 --- a/cicecore/cicedynB/general/ice_step_mod.F90 +++ b/cicecore/cicedynB/general/ice_step_mod.F90 @@ -46,6 +46,10 @@ module ice_step_mod step_snow, prep_radiation, step_radiation, ocean_mixed_layer, & update_state, biogeochemistry, step_dyn_wave, step_prep + real (kind=dbl_kind), dimension (:,:,:), allocatable :: & + uvelT_icep, & ! uvel for wind stress computation in icepack + vvelT_icep ! vvel for wind stress computation in icepack + !======================================================================= contains @@ -80,6 +84,13 @@ subroutine step_prep use ice_flux, only: uatm, vatm, uatmT, vatmT use ice_grid, only: grid_atm_dynu, grid_atm_dynv, grid_average_X2Y + use ice_state, only: uvel, vvel + + logical (kind=log_kind) :: & + highfreq ! highfreq flag + + logical (kind=log_kind), save :: & + first_call = .true. ! first call flag character(len=*), parameter :: subname = '(step_prep)' @@ -92,6 +103,26 @@ subroutine step_prep call grid_average_X2Y('S',uatm,grid_atm_dynu,uatmT,'T') call grid_average_X2Y('S',vatm,grid_atm_dynv,vatmT,'T') + !----------------------------------------------------------------- + ! Compute uvelT_icep, vvelT_icep + !----------------------------------------------------------------- + + if (first_call) then + allocate(uvelT_icep(nx_block,ny_block,max_blocks)) + allocate(vvelT_icep(nx_block,ny_block,max_blocks)) + uvelT_icep = c0 + vvelT_icep = c0 + endif + + call icepack_query_parameters(highfreq_out=highfreq) + + if (highfreq) then + call grid_average_X2Y('A', uvel, 'U', uvelT_icep, 'T') + call grid_average_X2Y('A', vvel, 'U', vvelT_icep, 'T') + endif + + first_call = .false. + end subroutine step_prep !======================================================================= @@ -209,7 +240,7 @@ subroutine step_therm1 (dt, iblk) Qa_iso, Qref_iso, fiso_evap, HDO_ocn, H2_16O_ocn, H2_18O_ocn use ice_grid, only: lmask_n, lmask_s, tmask use ice_state, only: aice, aicen, aicen_init, vicen_init, & - vice, vicen, vsno, vsnon, trcrn, uvelT, vvelT, vsnon_init + vice, vicen, vsno, vsnon, trcrn, vsnon_init #ifdef CICE_IN_NEMO use ice_state, only: aice_init #endif @@ -251,8 +282,6 @@ subroutine step_therm1 (dt, iblk) tr_pond_lvl, tr_pond_topo, calc_Tsfc, highfreq, tr_snow real (kind=dbl_kind) :: & - uvelTij, & ! cell-centered velocity, x component (m/s) - vvelTij, & ! cell-centered velocity, y component (m/s) puny ! a very small number real (kind=dbl_kind), dimension(n_aero,2,ncat) :: & @@ -336,14 +365,6 @@ subroutine step_therm1 (dt, iblk) do j = jlo, jhi do i = ilo, ihi - if (highfreq) then ! include ice velocity in calculation of wind stress - uvelTij = uvelT(i,j,iblk) - vvelTij = vvelT(i,j,iblk) - else - uvelTij = c0 - vvelTij = c0 - endif ! highfreq - if (tr_snow) then do n = 1, ncat do k = 1, nslyr @@ -389,8 +410,8 @@ subroutine step_therm1 (dt, iblk) vicen = vicen (i,j,:,iblk), & vsno = vsno (i,j, iblk), & vsnon = vsnon (i,j,:,iblk), & - uvel = uvelTij , & - vvel = vvelTij , & + uvel = uvelT_icep (i,j, iblk), & + vvel = vvelT_icep (i,j, iblk), & Tsfc = trcrn (i,j,nt_Tsfc,:,iblk), & zqsn = trcrn (i,j,nt_qsno:nt_qsno+nslyr-1,:,iblk), & zqin = trcrn (i,j,nt_qice:nt_qice+nilyr-1,:,iblk), & @@ -939,7 +960,7 @@ subroutine step_dyn_horiz (dt) use ice_flux, only: strocnxU, strocnyU, strocnxT_iavg, strocnyT_iavg use ice_flux, only: init_history_dyn use ice_grid, only: grid_average_X2Y - use ice_state, only: aiU, uvel, vvel, uvelT, vvelT + use ice_state, only: aiU use ice_transport_driver, only: advection, transport_upwind, transport_remap real (kind=dbl_kind), intent(in) :: & @@ -971,14 +992,6 @@ subroutine step_dyn_horiz (dt) if (kdyn == 2) call eap (dt) if (kdyn == 3) call implicit_solver (dt) - !----------------------------------------------------------------- - ! Compute uvelT, vvelT - ! only needed for highfreq, but compute anyway - !----------------------------------------------------------------- - - call grid_average_X2Y('A', uvel, 'U', uvelT, 'T') - call grid_average_X2Y('A', vvel, 'U', vvelT, 'T') - !----------------------------------------------------------------- ! Compute strocnxT_iavg, strocnyT_iavg for thermo and coupling !----------------------------------------------------------------- diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index b7082bf93..dfccdd413 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -2680,13 +2680,23 @@ subroutine grid_average_X2Y_NEversion(type,work1a,grid1a,work1b,grid1b,work2,gri ! state masked case('NE2US') - call grid_average_X2YS_2('NE2US',work1a,narea,npm,work1b,earea,epm,work2) + call grid_average_X2Y_2('NE2US',work1a,narea,npm,work1b,earea,epm,work2) case('EN2US') - call grid_average_X2YS_2('NE2US',work1b,narea,npm,work1a,earea,epm,work2) + call grid_average_X2Y_2('NE2US',work1b,narea,npm,work1a,earea,epm,work2) case('NE2TS') - call grid_average_X2YS_2('NE2TS',work1a,narea,npm,work1b,earea,epm,work2) + call grid_average_X2Y_2('NE2TS',work1a,narea,npm,work1b,earea,epm,work2) case('EN2TS') - call grid_average_X2YS_2('NE2TS',work1b,narea,npm,work1a,earea,epm,work2) + call grid_average_X2Y_2('NE2TS',work1b,narea,npm,work1a,earea,epm,work2) + + ! state unmasked + case('NE2UA') + call grid_average_X2Y_2('NE2UA',work1a,narea,npm,work1b,earea,epm,work2) + case('EN2UA') + call grid_average_X2Y_2('NE2UA',work1b,narea,npm,work1a,earea,epm,work2) + case('NE2TA') + call grid_average_X2Y_2('NE2TA',work1a,narea,npm,work1b,earea,epm,work2) + case('EN2TA') + call grid_average_X2Y_2('NE2TA',work1b,narea,npm,work1a,earea,epm,work2) case default call abort_ice(subname//'ERROR: unknown X2Y '//trim(X2Y)) @@ -3580,36 +3590,6 @@ subroutine grid_average_X2YF(dir,work1,wght1,work2,wght2) end subroutine grid_average_X2YF -!======================================================================= -! Compute the minimum of adjacent values of a field at specific indices, -! depending on the grid location (U, E, N) -! - real(kind=dbl_kind) function grid_neighbor_min(field, i, j, grid_location) result(mini) - - real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & - field ! field defined at T point - - integer (kind=int_kind), intent(in) :: & - i, j - - character(len=*), intent(in) :: & - grid_location ! grid location at which to compute the minumum (U, E, N) - - character(len=*), parameter :: subname = '(grid_neighbor_min)' - - select case (trim(grid_location)) - case('U') - mini = min(field(i,j), field(i+1,j), field(i,j+1), field(i+1,j+1)) - case('E') - mini = min(field(i,j), field(i+1,j)) - case('N') - mini = min(field(i,j), field(i,j+1)) - case default - call abort_ice(subname // ' unknown grid_location: ' // grid_location) - end select - - end function grid_neighbor_min - !======================================================================= ! Shifts quantities from one grid to another ! State masked version, simple weighted averager @@ -3618,7 +3598,7 @@ end function grid_neighbor_min ! ! author: T. Craig - subroutine grid_average_X2YS_2(dir,work1a,wght1a,mask1a,work1b,wght1b,mask1b,work2) + subroutine grid_average_X2Y_2(dir,work1a,wght1a,mask1a,work1b,wght1b,mask1b,work2) use ice_constants, only: c0 @@ -3645,7 +3625,7 @@ subroutine grid_average_X2YS_2(dir,work1a,wght1a,mask1a,work1b,wght1b,mask1b,wor type (block) :: & this_block ! block information for current block - character(len=*), parameter :: subname = '(grid_average_X2YS_2)' + character(len=*), parameter :: subname = '(grid_average_X2Y_2)' work2(:,:,:) = c0 @@ -3701,11 +3681,91 @@ subroutine grid_average_X2YS_2(dir,work1a,wght1a,mask1a,work1b,wght1b,mask1b,wor enddo !$OMP END PARALLEL DO + case('NE2UA') + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp) + 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 + wtmp = (wght1a(i ,j ,iblk) & + + wght1a(i+1,j ,iblk) & + + wght1b(i ,j ,iblk) & + + wght1b(i ,j+1,iblk)) + if (wtmp /= c0) & + work2(i,j,iblk) = (work1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & + + work1a(i+1,j ,iblk)*wght1a(i+1,j ,iblk) & + + work1b(i ,j ,iblk)*wght1b(i ,j ,iblk) & + + work1b(i ,j+1,iblk)*wght1b(i ,j+1,iblk)) & + / wtmp + enddo + enddo + enddo + !$OMP END PARALLEL DO + + case('NE2TA') + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,wtmp) + 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 + wtmp = (wght1a(i ,j-1,iblk) & + + wght1a(i ,j ,iblk) & + + wght1b(i-1,j ,iblk) & + + wght1b(i ,j ,iblk)) + if (wtmp /= c0) & + work2(i,j,iblk) = (work1a(i ,j-1,iblk)*wght1a(i ,j-1,iblk) & + + work1a(i ,j ,iblk)*wght1a(i ,j ,iblk) & + + work1b(i-1,j ,iblk)*wght1b(i-1,j ,iblk) & + + work1b(i ,j ,iblk)*wght1b(i ,j ,iblk)) & + / wtmp + enddo + enddo + enddo + !$OMP END PARALLEL DO + case default call abort_ice(subname//'ERROR: unknown option '//trim(dir)) end select - end subroutine grid_average_X2YS_2 + end subroutine grid_average_X2Y_2 + +!======================================================================= +! Compute the minimum of adjacent values of a field at specific indices, +! depending on the grid location (U, E, N) +! + real(kind=dbl_kind) function grid_neighbor_min(field, i, j, grid_location) result(mini) + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + field ! field defined at T point + + integer (kind=int_kind), intent(in) :: & + i, j + + character(len=*), intent(in) :: & + grid_location ! grid location at which to compute the minumum (U, E, N) + + character(len=*), parameter :: subname = '(grid_neighbor_min)' + + select case (trim(grid_location)) + case('U') + mini = min(field(i,j), field(i+1,j), field(i,j+1), field(i+1,j+1)) + case('E') + mini = min(field(i,j), field(i+1,j)) + case('N') + mini = min(field(i,j), field(i,j+1)) + case default + call abort_ice(subname // ' unknown grid_location: ' // grid_location) + end select + + end function grid_neighbor_min !======================================================================= ! Compute the maximum of adjacent values of a field at specific indices, diff --git a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 index 2c7d3d63c..bd5a49eaf 100644 --- a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 +++ b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 @@ -289,7 +289,7 @@ subroutine restartfile (ice_ic) use ice_grid, only: tmask, grid_type, grid_ice, grid_average_X2Y use ice_state, only: trcr_depend, aice, vice, vsno, trcr, & aice0, aicen, vicen, vsnon, trcrn, aice_init, uvel, vvel, & - uvelE, vvelE, uvelN, vvelN, uvelT, vvelT, & + uvelE, vvelE, uvelN, vvelN, & trcr_base, nt_strata, n_trcr_strata character (*), optional :: ice_ic @@ -403,9 +403,6 @@ subroutine restartfile (ice_ic) 'vvelN',1,diag,field_loc_Nface, field_type_vector) endif - call grid_average_X2Y('A', uvel, 'U', uvelT, 'T') - call grid_average_X2Y('A', vvel, 'U', vvelT, 'T') - !----------------------------------------------------------------- ! radiation fields !----------------------------------------------------------------- From 251ca48aa8797ac75329ec569cec84e39a094b0a Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Mon, 7 Nov 2022 19:02:57 -0500 Subject: [PATCH 2/2] Add wave-ice coupling to nuopc/cmeps driver (#782) * merge latest master (#4) * Isotopes for CICE (#423) Co-authored-by: apcraig Co-authored-by: David Bailey Co-authored-by: Elizabeth Hunke * updated orbital calculations needed for cesm * fixed problems in updated orbital calculations needed for cesm * update CICE6 to support coupling with UFS * put in changes so that both ufsatm and cesm requirements for potential temperature and density are satisfied * Convergence on ustar for CICE. (#452) (#5) * Add atmiter_conv to CICE * Add documentation * trigger build the docs Co-authored-by: David A. Bailey * update icepack submodule * Revert "update icepack submodule" This reverts commit e70d1abcbeb4351195a2b81c6ce3f623c936426c. * update comp_ice.backend with temporary ice_timers fix * Fix threading problem in init_bgc * Fix additional OMP problems * changes for coldstart running * Move the forapps directory * remove cesmcoupled ifdefs * Fix logging issues for NUOPC * removal of many cpp-ifdefs * fix compile errors * fixes to get cesm working * fixed white space issue * Add restart_coszen namelist option * update icepack submodule * change Orion to orion in backend remove duplicate print lines from ice_transport_driver * add -link_mpi=dbg to debug flags (#8) * cice6 compile (#6) * enable debug build. fix to remove errors * fix an error in comp_ice.backend.libcice * change Orion to orion for machine identification * changes for consistency w/ current emc-cice5 (#13) Update to emc/develop fork to current CICE consortium Co-authored-by: David A. Bailey Co-authored-by: Tony Craig Co-authored-by: Elizabeth Hunke Co-authored-by: Mariana Vertenstein Co-authored-by: apcraig Co-authored-by: Philippe Blain * Fixcommit (#14) Align commit history between emc/develop and cice-consortium/master * Update CICE6 for integration to S2S * add wcoss_dell_p3 compiler macro * update to icepack w/ debug fix * replace SITE with MACHINE_ID * update compile scripts * Support TACC stampede (#19) * update icepack * add ice_dyn_vp module to CICE_InitMod * update gitmodules, update icepack * Update CICE to consortium master (#23) updates include: * deprecate upwind advection (CICE-Consortium#508) * add implicit VP solver (CICE-Consortium#491) * update icepack * switch icepack branches * update to icepack master but set abort flag in ITD routine to false * update icepack * Update CICE to latest Consortium master (#26) update CICE and Icepack * changes the criteria for aborting ice for thermo-conservation errors * updates the time manager * fixes two bugs in ice_therm_mushy * updates Icepack to Consortium master w/ flip of abort flag for troublesome IC cases * add cice changes for zlvs (#29) * update icepack and pointer * update icepack and revert gitmodules * Fix history features - Fix bug in history time axis when sec_init is not zero. - Fix issue with time_beg and time_end uninitialized values. - Add support for averaging with histfreq='1' by allowing histfreq_n to be any value in that case. Extend and clean up construct_filename for history files. More could be done, but wanted to preserve backwards compatibility. - Add new calendar_sec2hms to converts daily seconds to hh:mm:ss. Update the calchk calendar unit tester to check this method - Remove abort test in bcstchk, this was just causing problems in regression testing - Remove known problems documentation about problems writing when istep=1. This issue does not exist anymore with the updated time manager. - Add new tests with hist_avg = false. Add set_nml.histinst. * revert set_nml.histall * fix implementation error * update model log output in ice_init * Fix QC issues - Add netcdf ststus checks and aborts in ice_read_write.F90 - Check for end of file when reading records in ice_read_write.F90 for ice_read_nc methods - Update set_nml.qc to better specify the test, turn off leap years since we're cycling 2005 data - Add check in c ice.t-test.py to make sure there is at least 1825 files, 5 years of data - Add QC run to base_suite.ts to verify qc runs to completion and possibility to use those results directly for QC validation - Clean up error messages and some indentation in ice_read_write.F90 * Update testing - Add prod suite including 10 year gx1prod and qc test - Update unit test compare scripts * update documentation * reset calchk to 100000 years * update evp1d test * update icepack * update icepack * add memory profiling (#36) * add profile_memory calls to CICE cap * update icepack * fix rhoa when lowest_temp is 0.0 * provide default value for rhoa when imported temp_height_lowest (Tair) is 0.0 * resolves seg fault when frac_grid=false and do_ca=true * update icepack submodule * Update CICE for latest Consortium master (#38) * Implement advanced snow physics in icepack and CICE * Fix time-stamping of CICE history files * Fix CICE history file precision * Use CICE-Consortium/Icepack master (#40) * switch to icepack master at consortium * recreate cap update branch (#42) * add debug_model feature * add required variables and calls for tr_snow * remove 2 extraneous lines * remove two log print lines that were removed prior to merge of driver updates to consortium * duplicate gitmodule style for icepack * Update CICE to latest Consortium/main (#45) * Update CICE to Consortium/main (#48) Update OpenMP directives as needed including validation via new omp_suite. Fixed OpenMP in dynamics. Refactored eap puny/pi lookups to improve scalar performance Update Tsfc implementation to make sure land blocks don't set Tsfc to freezing temp Update for sea bed stress calculations * fix comment, fix env for orion and hera * replace save_init with step_prep in CICE_RunMod * fixes for cgrid repro * remove added haloupdates * baselines pass with these extra halo updates removed * change F->S for ocean velocities and tilts * fix debug failure when grid_ice=C * compiling in debug mode using -init=snan,arrays requires initialization of variables * respond to review comments * remove inserted whitespace for uvelE,N and vvelE,N * Add wave-cice coupling; update to Consortium main (#51) * add wave-ice fields * initialize aicen_init, which turns up as NaN in calc of floediam export * add call to icepack_init_wave to initialize wavefreq and dwavefreq * update to latest consortium main (PR 752) * add initializationsin ice_state * initialize vsnon/vsnon_init and vicen/vicen_init Co-authored-by: apcraig Co-authored-by: David Bailey Co-authored-by: Elizabeth Hunke Co-authored-by: Mariana Vertenstein Co-authored-by: Minsuk Ji <57227195+MinsukJi-NOAA@users.noreply.github.com> Co-authored-by: Tony Craig Co-authored-by: Philippe Blain --- cicecore/cicedynB/general/ice_state.F90 | 6 + cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 | 9 +- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 22 ++-- .../drivers/nuopc/cmeps/ice_import_export.F90 | 112 ++++++++++++++++-- 4 files changed, 129 insertions(+), 20 deletions(-) diff --git a/cicecore/cicedynB/general/ice_state.F90 b/cicecore/cicedynB/general/ice_state.F90 index 7b20879bc..862f0a8bc 100644 --- a/cicecore/cicedynB/general/ice_state.F90 +++ b/cicecore/cicedynB/general/ice_state.F90 @@ -188,6 +188,12 @@ subroutine alloc_state n_trcr_strata = 0 nt_strata = 0 trcr_base = c0 + aicen = c0 + aicen_init = c0 + vicen = c0 + vicen_init = c0 + vsnon = c0 + vsnon_init = c0 end subroutine alloc_state diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index f9b5116d0..3d5e5cc2a 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -78,7 +78,7 @@ subroutine cice_init2() use ice_calendar , only: dt, dt_dyn, istep, istep1, write_ic, init_calendar, calendar use ice_communicate , only: my_task, master_task use ice_diagnostics , only: init_diags - use ice_domain_size , only: ncat, nfsd + use ice_domain_size , only: ncat, nfsd, nfreq use ice_dyn_eap , only: init_eap, alloc_dyn_eap use ice_dyn_shared , only: kdyn, init_dyn use ice_dyn_vp , only: init_vp @@ -94,10 +94,12 @@ subroutine cice_init2() use ice_restoring , only: ice_HaloRestore_init use ice_timers , only: timer_total, init_ice_timers, ice_timer_start use ice_transport_driver , only: init_transport + use ice_arrays_column , only: wavefreq, dwavefreq logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec, tr_snow character(len=char_len) :: snw_aging_table + real(kind=dbl_kind), dimension(25) :: wave_spectrum_profile ! hardwire for now character(len=*), parameter :: subname = '(cice_init2)' !---------------------------------------------------- @@ -177,6 +179,11 @@ subroutine cice_init2() endif endif + if (wave_spec) then + call icepack_init_wave(nfreq=nfreq, & + wave_spectrum_profile=wave_spectrum_profile, wavefreq=wavefreq, dwavefreq=dwavefreq) + end if + ! Initialize shortwave components using swdn from previous timestep ! if restarting. These components will be scaled to current forcing ! in prep_radiation. diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 8920ea386..182308973 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -25,6 +25,7 @@ module ice_comp_nuopc use ice_calendar , only : force_restart_now, write_ic 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 use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long use ice_fileunits , only : nu_diag, nu_diag_set, inst_index, inst_name use ice_fileunits , only : inst_suffix, release_all_fileunits, flush_fileunit @@ -80,9 +81,6 @@ module ice_comp_nuopc character(len=*) , parameter :: orb_variable_year = 'variable_year' character(len=*) , parameter :: orb_fixed_parameters = 'fixed_parameters' - character(len=*),parameter :: shr_cal_noleap = 'NO_LEAP' - character(len=*),parameter :: shr_cal_gregorian = 'GREGORIAN' - type(ESMF_Mesh) :: ice_mesh integer :: nthrds ! Number of threads to use in this component @@ -216,7 +214,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) type(ESMF_Time) :: stopTime ! Stop time type(ESMF_Time) :: refTime ! Ref time type(ESMF_TimeInterval) :: timeStep ! Model timestep - type(ESMF_Calendar) :: esmf_calendar ! esmf calendar type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type integer :: start_ymd ! Start date (YYYYMMDD) integer :: start_tod ! start time of day (s) @@ -339,7 +336,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call get_component_instance(gcomp, inst_suffix, inst_index, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - inst_name = "ICE"//trim(inst_suffix) +! inst_name = "ICE"//trim(inst_suffix) + inst_name = "ICE" !---------------------------------------------------------------------------- ! start cice timers @@ -470,9 +468,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ESMF_TimeGet( currTime, calkindflag=esmf_caltype, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (esmf_caltype == ESMF_CALKIND_NOLEAP) then - calendar_type = shr_cal_noleap + calendar_type = ice_calendar_noleap else if (esmf_caltype == ESMF_CALKIND_GREGORIAN) then - calendar_type = shr_cal_gregorian + calendar_type = ice_calendar_gregorian else call abort_ice( subname//'ERROR:: bad calendar for ESMF' ) end if @@ -581,9 +579,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) end if call icepack_query_parameters( tfrz_option_out=tfrz_option) if (tfrz_option_driver /= tfrz_option) then - write(errmsg,'(a)') trim(subname)//'error: tfrz_option from driver '//trim(tfrz_option_driver)//& - ' must be the same as tfrz_option from cice namelist '//trim(tfrz_option) - call abort_ice(trim(errmsg)) + write(errmsg,'(a)') trim(subname)//'WARNING: tfrz_option from driver '//trim(tfrz_option_driver)//& + ' is overwriting tfrz_option from cice namelist '//trim(tfrz_option) + write(nu_diag,*) trim(errmsg) + call icepack_warnings_flush(nu_diag) + call icepack_init_parameters(tfrz_option_in=tfrz_option_driver) endif ! Flux convergence tolerance - always use the driver attribute value @@ -594,7 +594,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) read(cvalue,*) atmiter_conv_driver call icepack_query_parameters( atmiter_conv_out=atmiter_conv) if (atmiter_conv_driver /= atmiter_conv) then - write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'warning: atmiter_ from driver ',& + write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'WARNING: atmiter_ from driver ',& atmiter_conv_driver,' is overwritting atmiter_conv from cice namelist ',atmiter_conv write(nu_diag,*) trim(errmsg) call icepack_warnings_flush(nu_diag) diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index 5a6ce7572..d95a4d9b2 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -3,12 +3,13 @@ module ice_import_export use ESMF use NUOPC use NUOPC_Model - use ice_kinds_mod , only : int_kind, dbl_kind, char_len, log_kind + use ice_kinds_mod , only : int_kind, dbl_kind, char_len, char_len_long, log_kind use ice_constants , only : c0, c1, spval_dbl, radius use ice_constants , only : field_loc_center, field_type_scalar, field_type_vector use ice_blocks , only : block, get_block, nx_block, ny_block use ice_domain , only : nblocks, blocks_ice, halo_info, distrb_info use ice_domain_size , only : nx_global, ny_global, block_size_x, block_size_y, max_blocks, ncat + use ice_domain_size , only : nfreq, nfsd use ice_exit , only : abort_ice use ice_flux , only : strairxT, strairyT, strocnxT_iavg, strocnyT_iavg use ice_flux , only : alvdr, alidr, alvdf, alidf, Tref, Qref, Uref @@ -23,9 +24,10 @@ module ice_import_export use ice_flux , only : fsnow, uocn, vocn, sst, ss_tltx, ss_tlty, frzmlt use ice_flux , only : send_i2x_per_cat use ice_flux , only : sss, Tf, wind, fsw - use ice_state , only : vice, vsno, aice, aicen_init, trcr + use ice_arrays_column , only : floe_rad_c, wave_spectrum + use ice_state , only : vice, vsno, aice, aicen_init, trcr, trcrn use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm - use ice_grid , only : grid_type, grid_average_X2Y + use ice_grid , only : grid_type use ice_mesh_mod , only : ocn_gridcell_frac use ice_boundary , only : ice_HaloUpdate use ice_fileunits , only : nu_diag, flush_fileunit @@ -34,8 +36,10 @@ 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_parameters , only : puny, c2 use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf #ifdef CESMCOUPLED use shr_frz_mod , only : shr_frz_freezetemp @@ -112,6 +116,7 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam character(char_len) :: stdname character(char_len) :: cvalue logical :: flds_wiso ! use case + logical :: flds_wave ! use case logical :: isPresent, isSet character(len=*), parameter :: subname='(ice_import_export:ice_advertise_fields)' !------------------------------------------------------------------------------- @@ -148,6 +153,17 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam write(nu_diag,*)'flds_wiso = ',flds_wiso end if + flds_wave = .false. + call NUOPC_CompAttributeGet(gcomp, name='wav_coupling_to_cice', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) flds_wave + end if + if (my_task == master_task) then + write(nu_diag,*)'flds_wave = ',flds_wave + end if + !----------------- ! advertise import fields !----------------- @@ -192,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) @@ -225,6 +249,10 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam call fldlist_add(fldsFrIce_num, fldsFrIce, 'ice_fraction_n', & ungridded_lbound=1, ungridded_ubound=ncat) end if + if (flds_wave) then + call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_thick' ) + call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_floediam' ) + end if ! ice/atm fluxes computed by ice call fldlist_add(fldsFrIce_num, fldsFrIce, 'stress_on_air_ice_zonal' ) @@ -292,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 + integer :: i, j, iblk, n, k 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(:) @@ -403,11 +431,10 @@ subroutine ice_import( importState, rc ) ! local variables integer,parameter :: nflds=16 integer,parameter :: nfldv=6 - integer :: i, j, iblk, n + integer :: i, j, iblk, n, k 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 @@ -559,6 +586,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 @@ -845,7 +895,7 @@ subroutine ice_export( exportState, rc ) ! local variables type(block) :: this_block ! block information for current block - integer :: i, j, iblk, n ! incides + integer :: i, j, iblk, n, k ! indices integer :: n2 ! thickness category index integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain real (kind=dbl_kind) :: workx, worky ! tmps for converting grid @@ -859,7 +909,11 @@ subroutine ice_export( exportState, rc ) real (kind=dbl_kind) :: tauxo (nx_block,ny_block,max_blocks) ! ice/ocean stress real (kind=dbl_kind) :: tauyo (nx_block,ny_block,max_blocks) ! ice/ocean stress 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), allocatable :: tempfld(:,:,:) real (kind=dbl_kind), pointer :: dataptr_ifrac_n(:,:) real (kind=dbl_kind), pointer :: dataptr_swpen_n(:,:) @@ -877,6 +931,9 @@ subroutine ice_export( exportState, rc ) ! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & ! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + call icepack_query_tracer_indices(nt_fsd_out=nt_fsd) + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -890,8 +947,10 @@ subroutine ice_export( exportState, rc ) tauya(:,:,:) = c0 tauxo(:,:,:) = c0 tauyo(:,:,:) = c0 + floediam(:,:,:) = c0 + floethick(:,:,:) = c0 - !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky, this_block, ilo, ihi, jlo, jhi) + !$OMP PARALLEL DO PRIVATE(iblk,i,j,k,workx,worky, this_block, ilo, ihi, jlo, jhi) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -904,6 +963,27 @@ subroutine ice_export( exportState, rc ) ! ice fraction ailohi(i,j,iblk) = min(aice(i,j,iblk), c1) + if (tr_fsd) then + ! floe thickness (m) + if (aice(i,j,iblk) > puny) then + floethick(i,j,iblk) = vice(i,j,iblk) / aice(i,j,iblk) + else + floethick(i,j,iblk) = c0 + end if + + ! floe diameter (m) + workx = c0 + worky = c0 + do n = 1, ncat + do k = 1, nfsd + workx = workx + floe_rad_c(k) * aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk) + worky = worky + aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk) + end do + end do + if (worky > c0) workx = c2*workx / worky + 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 ???) @@ -1054,6 +1134,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 ! ------