diff --git a/lanl_cice/drivers/cice/CICE_InitMod.F90 b/lanl_cice/drivers/cice/CICE_InitMod.F90 index 5b627bc..a84a3eb 100644 --- a/lanl_cice/drivers/cice/CICE_InitMod.F90 +++ b/lanl_cice/drivers/cice/CICE_InitMod.F90 @@ -172,12 +172,14 @@ subroutine init_restart use ice_brine, only: init_hbrine use ice_calendar, only: time, calendar use ice_domain, only: nblocks - use ice_domain_size, only: ncat + use ice_domain_size, only: ncat, max_ntrcr use ice_dyn_eap, only: read_restart_eap use ice_dyn_shared, only: kdyn use ice_firstyear, only: init_fy, restart_FY, read_restart_FY use ice_flux, only: sss + use ice_grid, only: tmask use ice_init, only: ice_ic + use ice_itd, only: aggregate use ice_lvl, only: init_lvl, restart_lvl, read_restart_lvl use ice_meltpond_cesm, only: init_meltponds_cesm, & restart_pond_cesm, read_restart_pond_cesm @@ -187,9 +189,12 @@ subroutine init_restart restart_pond_topo, read_restart_pond_topo use ice_restart_shared, only: runtype, restart use ice_restart_driver, only: restartfile, restartfile_v4 - use ice_state, only: tr_iage, tr_FY, tr_lvl, tr_pond_cesm, & - tr_pond_lvl, tr_pond_topo, tr_aero, trcrn, & - nt_iage, nt_FY, nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, tr_brine + use ice_state, only: aicen, vicen, vsnon, aice, trcr, & + vice, vsno, aice0, tr_iage, trcrn, trcr_depend, & + tr_brine, tr_aero, tr_pond_topo, tr_pond_lvl, tr_pond_cesm, & + tr_lvl, nt_alvl, nt_vlvl, tr_fy, nt_iage, nt_fy, & + nt_hpnd, nt_apnd, nt_ipnd + use ice_zbgc, only: init_bgc use ice_zbgc_shared, only: skl_bgc @@ -289,6 +294,29 @@ subroutine init_restart if (tr_brine) call init_hbrine ! brine height tracer if (skl_bgc) call init_bgc ! biogeochemistry + !----------------------------------------------------------------- + ! aggregate tracers + !----------------------------------------------------------------- + + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + + call aggregate (nx_block, ny_block, & + aicen(:,:,:,iblk), & + trcrn(:,:,:,:,iblk),& + vicen(:,:,:,iblk), & + vsnon(:,:,:,iblk), & + aice (:,:, iblk), & + trcr (:,:,:,iblk), & + vice (:,:, iblk), & + vsno (:,:, iblk), & + aice0(:,:, iblk), & + tmask(:,:, iblk), & + max_ntrcr, & + trcr_depend) + + enddo + !$OMP END PARALLEL DO end subroutine init_restart !======================================================================= diff --git a/lanl_cice/io_netcdf/ice_restart.F90 b/lanl_cice/io_netcdf/ice_restart.F90 index 69cdd31..0332b94 100644 --- a/lanl_cice/io_netcdf/ice_restart.F90 +++ b/lanl_cice/io_netcdf/ice_restart.F90 @@ -208,6 +208,7 @@ subroutine init_restart_write(filename_spec) call define_rest_field(ncid,'uvel',dims) call define_rest_field(ncid,'vvel',dims) + call define_rest_field(ncid,'coszen',dims) call define_rest_field(ncid,'scale_factor',dims) call define_rest_field(ncid,'swvdr',dims) @@ -395,6 +396,7 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3, & use ice_domain_size, only: max_blocks, ncat use ice_fileunits, only: nu_diag use ice_read_write, only: ice_read, ice_read_nc + use ice_communicate, only: my_task, master_task integer (kind=int_kind), intent(in) :: & nu , & ! unit number (not used for netcdf) @@ -425,22 +427,37 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3, & varid, & ! variable id status ! status variable from netCDF routine + logical (kind=log_kind) :: & + restart_ext2 ! temporary value to allow restart_ext for writing + ! but not reading if runtyp = initial + real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: & work2 ! input array (real, 8-byte) + ! set temporary variable + restart_ext2 = restart_ext + if(runtype .eq. 'initial' .and. restart_ext)restart_ext2 = .false. + if (my_task == master_task)then + if(restart_ext2)then + write(nu_diag,*)'NOTE: will read ghost cells in restart' + else + write(nu_diag,*)'NOTE: will not read ghost cells in restart' + end if + end if + if (restart_format == 'nc') then if (present(field_loc)) then if (ndim3 == ncat) then - if (restart_ext) then + if (restart_ext2) then call ice_read_nc(ncid,1,vname,work,diag, & - field_loc=field_loc,field_type=field_type,restart_ext=restart_ext) + field_loc=field_loc,field_type=field_type,restart_ext=restart_ext2) else call ice_read_nc(ncid,1,vname,work,diag,field_loc,field_type) endif elseif (ndim3 == 1) then - if (restart_ext) then + if (restart_ext2) then call ice_read_nc(ncid,1,vname,work2,diag, & - field_loc=field_loc,field_type=field_type,restart_ext=restart_ext) + field_loc=field_loc,field_type=field_type,restart_ext=restart_ext2) else call ice_read_nc(ncid,1,vname,work2,diag,field_loc,field_type) endif @@ -450,14 +467,14 @@ subroutine read_restart_field(nu,nrec,work,atype,vname,ndim3, & endif else if (ndim3 == ncat) then - if (restart_ext) then - call ice_read_nc(ncid, 1, vname, work, diag, restart_ext=restart_ext) + if (restart_ext2) then + call ice_read_nc(ncid, 1, vname, work, diag, restart_ext=restart_ext2) else call ice_read_nc(ncid, 1, vname, work, diag) endif elseif (ndim3 == 1) then - if (restart_ext) then - call ice_read_nc(ncid, 1, vname, work2, diag, restart_ext=restart_ext) + if (restart_ext2) then + call ice_read_nc(ncid, 1, vname, work2, diag, restart_ext=restart_ext2) else call ice_read_nc(ncid, 1, vname, work2, diag) endif diff --git a/lanl_cice/source/ice_flux.F90 b/lanl_cice/source/ice_flux.F90 index 2aab4e3..4bbb3a1 100644 --- a/lanl_cice/source/ice_flux.F90 +++ b/lanl_cice/source/ice_flux.F90 @@ -190,7 +190,8 @@ module ice_flux albice , & ! bare ice albedo albsno , & ! snow albedo albpnd , & ! melt pond albedo - apeff_ai ! effective pond area used for radiation calculation + apeff_ai , & ! effective pond area used for radiation calculation + snowfrac ! snow fraction used in radiation real (kind=dbl_kind), & dimension(nx_block,ny_block,max_blocks,max_nstrm), target, public :: & @@ -210,6 +211,9 @@ module ice_flux fswthruidr, & ! nir dir shortwave penetrating to ocean (W/m^2) fswthruidf ! nir dif shortwave penetrating to ocean (W/m^2) + real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: & + fswthrun_ai ! per-category fswthru * ai (W/m^2) + real (kind=dbl_kind), & dimension (nx_block,ny_block,max_aero,max_blocks), target, public :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) diff --git a/lanl_cice/source/ice_restart_driver.F90 b/lanl_cice/source/ice_restart_driver.F90 index 098a333..8ac691b 100644 --- a/lanl_cice/source/ice_restart_driver.F90 +++ b/lanl_cice/source/ice_restart_driver.F90 @@ -53,7 +53,7 @@ subroutine dumpfile(filename_spec) use ice_domain_size, only: nilyr, nslyr, ncat, max_blocks use ice_fileunits, only: nu_diag, nu_rst_pointer, nu_dump use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, & - strocnxT, strocnyT, sst, frzmlt, iceumask, & + strocnxT, strocnyT, sst, frzmlt, iceumask, coszen, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 @@ -125,6 +125,7 @@ subroutine dumpfile(filename_spec) !----------------------------------------------------------------- ! radiation fields !----------------------------------------------------------------- + call write_restart_field(nu_dump,0,coszen,'ruf8','coszen',1,diag) call write_restart_field(nu_dump,0,scale_factor,'ruf8','scale_factor',1,diag) call write_restart_field(nu_dump,0,swvdr,'ruf8','swvdr',1,diag) @@ -200,7 +201,7 @@ subroutine restartfile (ice_ic) max_ntrcr, max_blocks use ice_fileunits, only: nu_diag, nu_rst_pointer, nu_restart use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, & - strocnxT, strocnyT, sst, frzmlt, iceumask, & + strocnxT, strocnyT, sst, frzmlt, iceumask, coszen, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 @@ -327,6 +328,8 @@ subroutine restartfile (ice_ic) if (my_task == master_task) & write(nu_diag,*) 'radiation fields' + call read_restart_field(nu_restart,0,coszen,'ruf8', & + 'coszen',1,diag, field_loc_center, field_type_scalar) call read_restart_field(nu_restart,0,scale_factor,'ruf8', & 'scale_factor',1,diag, field_loc_center, field_type_scalar) call read_restart_field(nu_restart,0,swvdr,'ruf8', & diff --git a/lanl_cice/source/ice_shortwave.F90 b/lanl_cice/source/ice_shortwave.F90 index 9a71b89..e9b53b5 100644 --- a/lanl_cice/source/ice_shortwave.F90 +++ b/lanl_cice/source/ice_shortwave.F90 @@ -1,4 +1,4 @@ -! SVN:$Id$ +! SVN:$Id: ice_shortwave.F90 861 2014-10-21 16:44:30Z tcraig $ !======================================================================= ! ! The albedo and absorbed/transmitted flux parameterizations for @@ -51,7 +51,6 @@ module ice_shortwave use ice_communicate, only: my_task implicit none - save private public :: init_shortwave, run_dEdd, shortwave_ccsm3 @@ -108,13 +107,18 @@ module ice_shortwave public :: & fswpenln ! visible SW entering ice layers (W m-2) + real (kind=dbl_kind), dimension (nx_block,ny_block,ncat,max_blocks), & + public :: & + snowfracn ! Category snow fraction used in radiation + ! dEdd tuning parameters, set in namelist real (kind=dbl_kind), public :: & - R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo - R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo - R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo - dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C) - rsnw_mlt ! maximum melting snow grain radius (10^-6 m) + R_ice , & ! sea ice tuning parameter; +1 > 1sig increase in albedo + R_pnd , & ! ponded ice tuning parameter; +1 > 1sig increase in albedo + R_snw , & ! snow tuning parameter; +1 > ~.01 change in broadband albedo + dT_mlt, & ! change in temp for non-melt to melt snow grain radius change (C) + rsnw_mlt, & ! maximum melting snow grain radius (10^-6 m) + kalg ! algae absorption coefficient for 0.5 m thick layer real (kind=dbl_kind), parameter, public :: & hi_ssl = 0.050_dbl_kind, & ! ice surface scattering layer thickness (m) @@ -125,7 +129,7 @@ module ice_shortwave hp0 = 0.200_dbl_kind ! pond depth below which transition to bare ice real (kind=dbl_kind) :: & - exp_min ! minimum exponential value + exp_min ! minimum exponential value !======================================================================= @@ -142,7 +146,8 @@ subroutine init_shortwave use ice_flux, only: alvdf, alidf, alvdr, alidr, & alvdr_ai, alidr_ai, alvdf_ai, alidf_ai, & swvdr, swvdf, swidr, swidf, & - albice, albsno, albpnd, apeff_ai, albcnt, coszen, fsnow + albice, albsno, albpnd, albcnt, coszen, fsnow, & + apeff_ai, snowfrac use ice_orbital, only: init_orbit use ice_state, only: aicen, vicen, vsnon, trcrn, nt_Tsfc use ice_blocks, only: block, get_block @@ -161,36 +166,56 @@ subroutine init_shortwave ilo,ihi,jlo,jhi, & ! beginning and end of physical domain n ! thickness category index - real (kind=dbl_kind) :: cszn ! counter for history averaging + real (kind=dbl_kind) :: cszn,netsw ! counter for history averaging type (block) :: & this_block ! block information for current block - !$OMP PARALLEL DO PRIVATE(iblk,i,j) - do iblk=1,nblocks - do j = 1, ny_block - do i = 1, nx_block - alvdf(i,j,iblk) = c0 - alidf(i,j,iblk) = c0 - alvdr(i,j,iblk) = c0 - alidr(i,j,iblk) = c0 - alvdr_ai(i,j,iblk) = c0 - alidr_ai(i,j,iblk) = c0 - alvdf_ai(i,j,iblk) = c0 - alidf_ai(i,j,iblk) = c0 - enddo - enddo - enddo + !$OMP PARALLEL DO PRIVATE(iblk,i,j,n) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + alvdf(i,j,iblk) = c0 + alidf(i,j,iblk) = c0 + alvdr(i,j,iblk) = c0 + alidr(i,j,iblk) = c0 + alvdr_ai(i,j,iblk) = c0 + alidr_ai(i,j,iblk) = c0 + alvdf_ai(i,j,iblk) = c0 + alidf_ai(i,j,iblk) = c0 + enddo + enddo + + ! Initialize + do n = 1, ncat + do j = 1, ny_block + do i = 1, nx_block + alvdrn(i,j,n,iblk) = c0 + alidrn(i,j,n,iblk) = c0 + alvdfn(i,j,n,iblk) = c0 + alidfn(i,j,n,iblk) = c0 + fswsfcn(i,j,n,iblk) = c0 + fswintn(i,j,n,iblk) = c0 + fswthrun(i,j,n,iblk) = c0 + enddo ! i + enddo ! j + enddo ! ncat + + fswpenln(:,:,:,:,iblk) = c0 + Iswabsn(:,:,:,:,iblk) = c0 + Sswabsn(:,:,:,:,iblk) = c0 + + enddo ! iblk !$OMP END PARALLEL DO if (trim(shortwave) == 'dEdd') then ! delta Eddington -#ifndef CCSMCOUPLED +#ifndef CESMCOUPLED ! These come from the driver in the coupled model. call init_orbit ! initialize orbital parameters #endif !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) - do iblk=1,nblocks + do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -213,11 +238,13 @@ subroutine init_shortwave fswthrun(:,:,:,iblk), & fswthrunvdr(:,:,:,iblk), fswthrunvdf(:,:,:,iblk), & fswthrunidr(:,:,:,iblk), fswthrunidf(:,:,:,iblk), & - fswpenln(:,:,:,:,iblk), & + fswpenln(:,:,:,:,iblk), & Sswabsn(:,:,:,:,iblk), Iswabsn(:,:,:,:,iblk), & albicen(:,:,:,iblk), albsnon(:,:,:,iblk), & albpndn(:,:,:,iblk), apeffn(:,:,:,iblk), & - dhsn(:,:,:,iblk), ffracn(:,:,:,iblk)) + snowfracn(:,:,:,iblk), & + dhsn(:,:,:,iblk), ffracn(:,:,:,iblk), & + initonly = .true. ) enddo !$OMP END PARALLEL DO @@ -225,7 +252,7 @@ subroutine init_shortwave else ! basic (ccsm3) shortwave !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block) - do iblk=1,nblocks + do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -261,8 +288,8 @@ subroutine init_shortwave !----------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblk,i,j,n,ilo,ihi,jlo,jhi,this_block, & - !$OMP ij,icells,cszn) - do iblk=1,nblocks + !$OMP ij,icells,cszn,netsw,indxi,indxj) + do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo ihi = this_block%ihi @@ -295,7 +322,8 @@ subroutine init_shortwave alidr(i,j,iblk) = alidr(i,j,iblk) & + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk) - if (coszen(i,j,iblk) > puny) then ! sun above horizon + netsw = swvdr(i,j,iblk)+swidr(i,j,iblk)+swvdf(i,j,iblk)+swidf(i,j,iblk) + if (netsw > puny) then ! sun above horizon albice(i,j,iblk) = albice(i,j,iblk) & + albicen(i,j,n,iblk)*aicen(i,j,n,iblk) albsno(i,j,iblk) = albsno(i,j,iblk) & @@ -306,6 +334,8 @@ subroutine init_shortwave apeff_ai(i,j,iblk) = apeff_ai(i,j,iblk) & + apeffn(i,j,n,iblk)*aicen(i,j,n,iblk) + snowfrac(i,j,iblk) = snowfrac(i,j,iblk) & + + snowfracn(i,j,n,iblk)*aicen(i,j,n,iblk) enddo enddo ! ncat @@ -320,13 +350,6 @@ subroutine init_shortwave alidf_ai (i,j,iblk) = alidf (i,j,iblk) alvdr_ai (i,j,iblk) = alvdr (i,j,iblk) alidr_ai (i,j,iblk) = alidr (i,j,iblk) - - ! for history averaging - cszn = c0 - if (coszen(i,j,iblk) > puny) cszn = c1 - do n = 1, nstreams - albcnt(i,j,iblk,n) = albcnt(i,j,iblk,n) + cszn - enddo enddo enddo @@ -372,21 +395,21 @@ subroutine shortwave_ccsm3 (nx_block, ny_block, & swidf ! sw down, near IR, diffuse (W/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), & - intent(out) :: & - alvdrn , & ! visible, direct, avg (fraction) - alidrn , & ! near-ir, direct, avg (fraction) - alvdfn , & ! visible, diffuse, avg (fraction) - alidfn , & ! near-ir, diffuse, avg (fraction) - fswsfc , & ! SW absorbed at ice/snow surface (W m-2) - fswint , & ! SW absorbed in ice interior, below surface (W m-2) - fswthru , & ! SW through ice to ocean (W m-2) - fswthruvdr , & ! visible direct SW through ice to ocean (W m-2) - fswthruvdf , & ! visible diffuse SW through ice to ocean (W m-2) - albin , & ! bare ice albedo - albsn ! snow albedo + intent(inout) :: & + alvdrn , & ! visible, direct, avg (fraction) + alidrn , & ! near-ir, direct, avg (fraction) + alvdfn , & ! visible, diffuse, avg (fraction) + alidfn , & ! near-ir, diffuse, avg (fraction) + fswsfc , & ! SW absorbed at ice/snow surface (W m-2) + fswint , & ! SW absorbed in ice interior, below surface (W m-2) + fswthru , & ! SW through ice to ocean (W m-2) + fswthruvdr , & ! visible direct SW through ice to ocean (W m-2) + fswthruvdf , & ! visible diffuse SW through ice to ocean (W m-2) + albin , & ! bare ice albedo + albsn ! snow albedo real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr+1,ncat), & - intent(out) :: & + intent(inout) :: & fswpenl ! SW entering ice layers (W m-2) real (kind=dbl_kind), dimension (nx_block,ny_block), & @@ -394,11 +417,11 @@ subroutine shortwave_ccsm3 (nx_block, ny_block, & coszen ! cosine(zenith angle) real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr,ncat), & - intent(out) :: & + intent(inout) :: & Iswabs ! SW absorbed in particular layer (W m-2) real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr,ncat), & - intent(out) :: & + intent(inout) :: & Sswabs ! SW absorbed in particular layer (W m-2) ! local variables @@ -903,8 +926,8 @@ subroutine absorbed_solar (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block) :: & fswpen , & ! SW penetrating beneath surface (W m-2) - fswpenvdr , & ! SW vis dir penetrating beneath surface (W m-2) - fswpenvdf , & ! SW vis dif penetrating beneath surface (W m-2) + fswpenvdr , & ! SW vis dir penetrating beneath surface (W m-2) + fswpenvdf , & ! SW vis dif penetrating beneath surface (W m-2) trantop , & ! transmitted frac of penetrating SW at layer top tranbot ! transmitted frac of penetrating SW at layer bot @@ -929,8 +952,6 @@ subroutine absorbed_solar (nx_block, ny_block, & fswthruvdr(i,j) = c0 fswthruvdf(i,j) = c0 fswpen (i,j) = c0 - fswpenvdr (i,j) = c0 - fswpenvdf (i,j) = c0 trantop(i,j) = c0 tranbot(i,j) = c0 enddo @@ -971,7 +992,6 @@ subroutine absorbed_solar (nx_block, ny_block, & + (c1-alidfns(i,j))*asnow ) swabs = swabsv + swabsi - if (hs <= hs_min) swabs = c0 fswpenvdr(i,j) = swvdr(i,j) * (c1-alvdrni(i,j)) * (c1-asnow) * i0vis fswpenvdf(i,j) = swvdf(i,j) * (c1-alvdfni(i,j)) * (c1-asnow) * i0vis @@ -1087,16 +1107,18 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & fswthrun, & fswthrunvdr, fswthrunvdf, & fswthrunidr, fswthrunidf, & - fswpenln, & + fswpenln, & Sswabsn, Iswabsn, & albicen, albsnon, & albpndn, apeffn, & - dhsn, ffracn) + snowfracn, & + dhsn, ffracn, & + initonly ) use ice_calendar, only: dt use ice_meltpond_cesm, only: hs0 use ice_meltpond_topo, only: hp1 - use ice_meltpond_lvl, only: hs1, pndaspect, snowinfil + use ice_meltpond_lvl, only: hs1, pndaspect use ice_orbital, only: compute_coszen use ice_state, only: ntrcr, nt_Tsfc, nt_alvl, nt_apnd, nt_hpnd, nt_ipnd, & tr_pond_cesm, tr_pond_lvl, tr_pond_topo @@ -1132,7 +1154,7 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & real(kind=dbl_kind), dimension(nx_block,ny_block), intent(out) :: & coszen ! cosine solar zenith angle, < 0 for sun below horizon - real(kind=dbl_kind), dimension(nx_block,ny_block,ncat), intent(out) :: & + real(kind=dbl_kind), dimension(nx_block,ny_block,ncat), intent(inout) :: & alvdrn, & ! visible direct albedo (fraction) alvdfn, & ! near-ir direct albedo (fraction) alidrn, & ! visible diffuse albedo (fraction) @@ -1147,17 +1169,21 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & albicen, & ! albedo bare ice albsnon, & ! albedo snow albpndn, & ! albedo pond - apeffn ! effective pond area used for radiation calculation + apeffn, & ! effective pond area used for radiation calculation + snowfracn ! Snow fraction used in radiation - real(kind=dbl_kind), dimension(nx_block,ny_block,nslyr,ncat), intent(out) :: & + real(kind=dbl_kind), dimension(nx_block,ny_block,nslyr,ncat), intent(inout) :: & Sswabsn ! SW radiation absorbed in snow layers (W m-2) - real(kind=dbl_kind), dimension(nx_block,ny_block,nilyr,ncat), intent(out) :: & + real(kind=dbl_kind), dimension(nx_block,ny_block,nilyr,ncat), intent(inout) :: & Iswabsn ! SW radiation absorbed in ice layers (W m-2) - real(kind=dbl_kind), dimension(nx_block,ny_block,nilyr+1,ncat), intent(out) :: & + real(kind=dbl_kind), dimension(nx_block,ny_block,nilyr+1,ncat), intent(inout) :: & fswpenln ! visible SW entering ice layers (W m-2) + logical (kind=log_kind), optional :: & + initonly ! flag to indicate init only, default is false + ! local temporary variables integer (kind=int_kind) :: & @@ -1195,9 +1221,17 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & spn , & ! snow depth on refrozen pond (m) tmp ! 0 or 1 + logical (kind=log_kind) :: & + linitonly ! local initonly value + real (kind=dbl_kind), parameter :: & argmax = c10 ! maximum argument of exponential + linitonly = .false. + if (present(initonly)) then + linitonly = initonly + endif + exp_min = exp(-argmax) ! identify ice-ocean cells @@ -1243,9 +1277,11 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & trcrn(:,:,nt_Tsfc,n), fsn, hsn, & rhosnwn, rsnwn) + apeffn(:,:,n) = c0 ! for history + snowfracn(:,:,n) = c0 ! for history + ! set pond properties if (tr_pond_cesm) then - apeffn(:,:,n) = c0 ! for history do ij = 1, icells i = indxi(ij) j = indxj(ij) @@ -1260,27 +1296,31 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & fpn(i,j) = (c1 - asnow) * fpn(i,j) hpn(i,j) = pndaspect * fpn(i,j) endif + + ! Zero out fraction of thin ponds for radiation only + if (hpn(i,j) < hpmin) fpn(i,j) = c0 + apeffn(i,j,n) = fpn(i,j) ! for history enddo elseif (tr_pond_lvl) then - apeffn(:,:,n) = c0 ! for history do ij = 1, icells i = indxi(ij) j = indxj(ij) fpn(i,j) = c0 ! fraction of ice covered in pond hpn(i,j) = c0 ! pond depth over fpn + ! refrozen pond lid thickness avg over ice ! allow snow to cover pond ice ipn = trcrn(i,j,nt_alvl,n) * trcrn(i,j,nt_apnd,n) & * trcrn(i,j,nt_ipnd,n) dhs = dhsn(i,j,n) ! snow depth difference, sea ice - pond - if (ipn > puny .and. & + if (.not. linitonly .and. ipn > puny .and. & dhs < puny .and. fsnow(i,j)*dt > hs_min) & dhs = hsn(i,j) - fsnow(i,j)*dt ! initialize dhs>0 spn = hsn(i,j) - dhs ! snow depth on pond ice - if (ipn*spn < puny) dhs = c0 + if (.not. linitonly .and. ipn*spn < puny) dhs = c0 dhsn(i,j,n) = dhs ! save: constant until reset to 0 ! not using ipn assumes that lid ice is perfectly clear @@ -1303,7 +1343,7 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & ! infiltrate snow hp = hpn(i,j) - if (snowinfil .and. hp > puny) then + if (hp > puny) then hs = hsn(i,j) rp = rhofresh*hp/(rhofresh*hp + rhos*hs) if (rp < p15) then @@ -1318,22 +1358,24 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & hpn(i,j) = hp * tmp fpn(i,j) = fpn(i,j) * tmp endif - fsn(i,j) = min(fsn(i,j), c1-fpn(i,j)) + endif ! hp > puny + + ! Zero out fraction of thin ponds for radiation only + if (hpn(i,j) < hpmin) fpn(i,j) = c0 - endif ! snowinfil + fsn(i,j) = min(fsn(i,j), c1-fpn(i,j)) ! endif ! masking by lid ice apeffn(i,j,n) = fpn(i,j) ! for history enddo ! ij elseif (tr_pond_topo) then - apeffn(:,:,n) = c0 ! for history do ij = 1, icells i = indxi(ij) j = indxj(ij) ! Lid effective if thicker than hp1 if (trcrn(i,j,nt_apnd,n)*aicen(i,j,n) > puny .and. & - trcrn(i,j,nt_ipnd,n) < hp1) then + trcrn(i,j,nt_ipnd,n) < hp1) then fpn(i,j) = trcrn(i,j,nt_apnd,n) else fpn(i,j) = c0 @@ -1344,6 +1386,14 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & fpn(i,j) = c0 hpn(i,j) = c0 endif + + ! Zero out fraction of thin ponds for radiation only + if (hpn(i,j) < hpmin) fpn(i,j) = c0 + + ! If ponds are present snow fraction reduced to + ! non-ponded part dEdd scheme + fsn(i,j) = min(fsn(i,j), c1-fpn(i,j)) + apeffn(i,j,n) = fpn(i,j) enddo else @@ -1353,11 +1403,13 @@ subroutine run_dEdd(ilo,ihi,jlo,jhi, & trcrn(:,:,nt_Tsfc,n), & fsn, fpn, & hpn) - apeffn(:,:,n) = fpn(i,j) ! for history + apeffn(:,:,n) = fpn(:,:) ! for history fpn = c0 hpn = c0 endif + snowfracn(:,:,n) = fsn(:,:) ! for history + call shortwave_dEdd(nx_block, ny_block, & ntrcr, icells, & indxi, indxj, & @@ -1452,7 +1504,6 @@ subroutine shortwave_dEdd (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), & intent(in) :: & - coszen , & ! cosine of solar zenith angle aice , & ! concentration of ice vice , & ! volume of ice hs , & ! snow depth @@ -1477,7 +1528,8 @@ subroutine shortwave_dEdd (nx_block, ny_block, & swidf ! sw down, near IR, diffuse (W/m^2) real (kind=dbl_kind), dimension (nx_block,ny_block), & - intent(out) :: & + intent(inout) :: & + coszen , & ! cosine of solar zenith angle alvdr , & ! visible, direct, albedo (fraction) alvdf , & ! visible, diffuse, albedo (fraction) alidr , & ! near-ir, direct, albedo (fraction) @@ -1491,15 +1543,15 @@ subroutine shortwave_dEdd (nx_block, ny_block, & fswthruidf ! SW through snow/bare ice/ponded ice into ocean (W m-2) real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr+1), & - intent(out) :: & + intent(inout) :: & fswpenl ! visible SW entering ice layers (W m-2) real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr), & - intent(out) :: & + intent(inout) :: & Sswabs ! SW absorbed in snow layer (W m-2) real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr), & - intent(out) :: & + intent(inout) :: & Iswabs ! SW absorbed in ice layer (W m-2) real (kind=dbl_kind), dimension (nx_block,ny_block), & @@ -1557,6 +1609,8 @@ subroutine shortwave_dEdd (nx_block, ny_block, & aidrl , & ! near-ir, direct, albedo (fraction) aidfl ! near-ir, diffuse, albedo (fraction) + real (kind=dbl_kind) :: netsw + !----------------------------------------------------------------------- do j = 1, ny_block @@ -1608,12 +1662,13 @@ subroutine shortwave_dEdd (nx_block, ny_block, & i = indxi(ij) j = indxj(ij) vsno = hs(i,j) * aice(i,j) - if (coszen(i,j) > puny) then ! sun above horizon + netsw = swvdr(i,j)+swidr(i,j)+swvdf(i,j)+swidf(i,j) + if (netsw > puny) then ! sun above horizon aero_mp(i,j,na ) = trcr(i,j,nt_aero-1+na )*vsno aero_mp(i,j,na+1) = trcr(i,j,nt_aero-1+na+1)*vsno aero_mp(i,j,na+2) = trcr(i,j,nt_aero-1+na+2)*vice(i,j) aero_mp(i,j,na+3) = trcr(i,j,nt_aero-1+na+3)*vice(i,j) - endif ! aice > 0 and coszen > 0 + endif ! aice > 0 and netsw > 0 enddo ! ij enddo ! na endif ! if aerosols @@ -1630,7 +1685,9 @@ subroutine shortwave_dEdd (nx_block, ny_block, & i = indxi(ij) j = indxj(ij) ! sea ice points with sun above horizon - if (coszen(i,j) > puny) then + netsw = swvdr(i,j)+swidr(i,j)+swvdf(i,j)+swidf(i,j) + if (netsw > puny) then + coszen(i,j) = max(puny,coszen(i,j)) ! evaluate sea ice thickness and fraction hi(i,j) = vice(i,j) / aice(i,j) fi(i,j) = c1 - fs(i,j) - fp(i,j) @@ -1641,7 +1698,7 @@ subroutine shortwave_dEdd (nx_block, ny_block, & indxj_DE(icells_DE) = j ! bare ice endif ! fi > 0 - endif ! coszen > 0 + endif ! netsw > 0 enddo ! ij ! calculate bare sea ice @@ -1685,7 +1742,9 @@ subroutine shortwave_dEdd (nx_block, ny_block, & i = indxi(ij) j = indxj(ij) ! sea ice points with sun above horizon - if (coszen(i,j) > puny) then + netsw = swvdr(i,j)+swidr(i,j)+swvdf(i,j)+swidf(i,j) + if (netsw > puny) then + coszen(i,j) = max(puny,coszen(i,j)) ! snow-covered sea ice points if(fs(i,j) > c0) then icells_DE = icells_DE + 1 @@ -1693,7 +1752,7 @@ subroutine shortwave_dEdd (nx_block, ny_block, & indxj_DE(icells_DE) = j ! snow-covered ice endif ! fs > 0 - endif ! coszen > 0 + endif ! netsw > 0 enddo ! ij ! calculate snow covered sea ice @@ -1738,16 +1797,18 @@ subroutine shortwave_dEdd (nx_block, ny_block, & j = indxj(ij) hi(i,j) = c0 ! sea ice points with sun above horizon - if (coszen(i,j) > puny) then + netsw = swvdr(i,j)+swidr(i,j)+swvdf(i,j)+swidf(i,j) + if (netsw > puny) then + coszen(i,j) = max(puny,coszen(i,j)) hi(i,j) = vice(i,j) / aice(i,j) ! if non-zero pond fraction and sufficient pond depth - if( fp(i,j) > puny .and. hp(i,j) > hpmin ) then + if( fp(i,j) > c0 ) then icells_DE = icells_DE + 1 indxi_DE(icells_DE) = i indxj_DE(icells_DE) = j ! ponded ice endif - endif ! coszen > puny + endif ! netsw > puny enddo ! ij ! calculate ponded ice @@ -1782,6 +1843,19 @@ subroutine shortwave_dEdd (nx_block, ny_block, & + awtvdf*avdfl(i,j) + awtidf*aidfl(i,j) enddo +! If no incoming shortwave, set albedos to 1. + do j = 1, ny_block + do i = 1, nx_block + netsw = swvdr(i,j)+swidr(i,j)+swvdf(i,j)+swidf(i,j) + if (netsw <= puny) then + alvdr(i,j) = c1 + alvdf(i,j) = c1 + alidr(i,j) = c1 + alidf(i,j) = c1 + endif + enddo + enddo + dbug = .false. if (dbug .and. print_points) then do n = 1, npnt @@ -1816,7 +1890,7 @@ subroutine shortwave_dEdd (nx_block, ny_block, & alvdr(i,j),alvdf(i,j) write(nu_diag,*) ' alidr alidf = ', & alidr(i,j),alidf(i,j) - write(nu_diag,*) ' fswsfc fswint fswthru = ', & + write(nu_diag,*) ' fswsfc fswint fswthru = ', & fswsfc(i,j),fswint(i,j),fswthru(i,j) write(nu_diag,*) ' fswthruvdr fswthruvdf fswthruidr fswthruidf = ', & fswthruvdr(i,j),fswthruvdf(i,j), & @@ -1867,7 +1941,7 @@ subroutine compute_dEdd & fswthru, & fswthruvdr, fswthruvdf, & fswthruidr, fswthruidf, & - Sswabs, Iswabs, fswpenl) + Sswabs, Iswabs, fswpenl) use ice_therm_shared, only: heat_capacity use ice_state, only: tr_aero @@ -2370,9 +2444,6 @@ subroutine compute_dEdd & rhoi = 917.0_dbl_kind,& ! pure ice mass density (kg/m3) fr_max = 1.00_dbl_kind, & ! snow grain adjustment factor max fr_min = 0.80_dbl_kind, & ! snow grain adjustment factor min - ! algae absorption coefficient for 0.5 m thick layer - kalg = 0.60_dbl_kind, & ! for 0.5 m path of 75 mg Chl a / m2 -!cesm turns kalg off kalg = 0.00_dbl_kind, & ! for 0.5 m path of 75 mg Chl a / m2 ! tuning parameters ! ice and pond scat coeff fractional change for +- one-sigma in albedo fp_ice = 0.15_dbl_kind, & ! ice fraction of scat coeff for + stn dev in alb @@ -3073,7 +3144,6 @@ subroutine compute_dEdd & i = indxi_DE(ij) j = indxj_DE(ij) Iswabs(i,j,k) = Iswabs(i,j,k) + Iabs(ij,k)*fi(i,j) - ! bgc layer fswpenl(i,j,k) = fswpenl(i,j,k) + fthrul(ij,k)* fi(i,j) if (k == nilyr) then @@ -3360,17 +3430,6 @@ subroutine solution_dEdd & !----------------------------------------------------------------------- - ! initialize all layer apparent optical properties to 0 - do k = 0, klev - rdir (k) = c0 - rdif_a(k) = c0 - rdif_b(k) = c0 - tdir (k) = c0 - tdif_a(k) = c0 - tdif_b(k) = c0 - trnlay(k) = c0 - enddo - ! initialize all output to 0 do ij = 1, icells_DE i = indxi_DE(ij) @@ -3422,8 +3481,17 @@ subroutine solution_dEdd & j = indxj_DE(ij) ! begin main level loop - do k=0,klev + do k = 0, klev + ! initialize all layer apparent optical properties to 0 + rdir (k) = c0 + rdif_a(k) = c0 + rdif_b(k) = c0 + tdir (k) = c0 + tdif_a(k) = c0 + tdif_b(k) = c0 + trnlay(k) = c0 + ! compute next layer Delta-eddington solution only if total transmission ! of radiation to the interface just above the layer exceeds trmin. diff --git a/lanl_cice/source/ice_step_mod.F90 b/lanl_cice/source/ice_step_mod.F90 index cbe1680..a7bbfeb 100644 --- a/lanl_cice/source/ice_step_mod.F90 +++ b/lanl_cice/source/ice_step_mod.F90 @@ -1347,7 +1347,8 @@ subroutine step_radiation (dt, iblk) Sswabsn, Iswabsn, shortwave, & albicen, albsnon, albpndn, & alvdrn, alidrn, alvdfn, alidfn, & - run_dedd, shortwave_ccsm3, apeffn + run_dedd, shortwave_ccsm3, apeffn, & + snowfracn use ice_state, only: aicen, vicen, vsnon, trcrn, nt_Tsfc use ice_timers, only: ice_timer_start, ice_timer_stop, timer_sw @@ -1417,6 +1418,7 @@ subroutine step_radiation (dt, iblk) Sswabsn(:,:,:,:,iblk), Iswabsn(:,:,:,:,iblk), & albicen(:,:,:,iblk), albsnon(:,:,:,iblk), & albpndn(:,:,:,iblk), apeffn(:,:,:,iblk), & + snowfracn(:,:,:,iblk), & dhsn(:,:,:,iblk), ffracn(:,:,:,iblk)) else ! .not. dEdd diff --git a/lanl_cice_cap/cice_cap.F90 b/lanl_cice_cap/cice_cap.F90 index 2eead12..4567b78 100644 --- a/lanl_cice_cap/cice_cap.F90 +++ b/lanl_cice_cap/cice_cap.F90 @@ -26,7 +26,7 @@ module cice_cap_mod use ice_domain_size, only: max_blocks, nx_global, ny_global use ice_domain, only: nblocks, blocks_ice, halo_info, distrb_info use ice_distribution, only: ice_distributiongetblockloc - use ice_constants, only: Tffresh, rad_to_deg + use ice_constants, only: c0, Tffresh, rad_to_deg, depressT use ice_calendar, only: dt use ice_flux use ice_grid, only: TLAT, TLON, ULAT, ULON, hm, tarea, ANGLET, ANGLE, & @@ -902,6 +902,8 @@ subroutine ModelAdvance_slow(gcomp, rc) real(ESMF_KIND_R8), pointer :: dataPtr_vice(:,:,:) real(ESMF_KIND_R8), pointer :: dataPtr_vsno(:,:,:) #endif + character(240) :: import_timestr, export_timestr + character(240) :: fname character(240) :: msgString character(len=*),parameter :: subname='(cice_cap:ModelAdvance_slow)' @@ -958,10 +960,13 @@ subroutine ModelAdvance_slow(gcomp, rc) file=__FILE__)) & return ! bail out - if(write_diagnostics) then - import_slice = import_slice + 1 + call ESMF_TimeGet(currTime, timestring=import_timestr, rc=rc) + call ESMF_TimeGet(currTime+timestep, timestring=export_timestr, rc=rc) + if(write_diagnostics) then call state_diagnose(importState, 'cice_import', rc) + + fname = 'field_ice_import_'//trim(import_timestr)//'.nc' do i = 1,fldsToice_num fldname = fldsToice(i)%shortname call ESMF_StateGet(importState, itemName=trim(fldname), itemType=itemType, rc=rc) @@ -976,8 +981,8 @@ subroutine ModelAdvance_slow(gcomp, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return #ifdef CMEPS - call ESMF_FieldWrite(lfield, fileName='field_ice_import_'//trim(fldname)//'.nc', & - timeslice=import_slice, overwrite=overwrite_timeslice, rc=rc) + call ESMF_FieldWrite(lfield, fileName=trim(fname), & + timeslice=1, overwrite=overwrite_timeslice, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -998,8 +1003,8 @@ subroutine ModelAdvance_slow(gcomp, rc) fldptr2d(:,:) = fldptr(:,:,1) - call ESMF_FieldWrite(lfield2d, fileName='field_ice_import_'//trim(fldname)//'.nc', & - timeslice=import_slice, overwrite=overwrite_timeslice, rc=rc) + call ESMF_FieldWrite(lfield2d, fileName=trim(fname), & + timeslice=1, overwrite=overwrite_timeslice, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -1055,9 +1060,14 @@ subroutine ModelAdvance_slow(gcomp, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) return call State_getFldPtr(importState,'inst_height_lowest',dataPtr_zlvl,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) return - call State_getFldPtr(importState,'air_density_height_lowest',dataPtr_rhoabot,rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) return + ! for CMEPS, the air density is no longer sent from the mediator. For NEMS, this field can still + ! be advertised and sent (so that no changes are required to NEMS), but this field is now calculated + ! below with a calculation of air density equivalent to the mediator calculation + !call State_getFldPtr(importState,'air_density_height_lowest',dataPtr_rhoabot,rc=rc) + !if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) return + rhoa = 0.0_ESMF_KIND_R8 + potT = 0.0_ESMF_KIND_R8 do iblk = 1,nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -1071,9 +1081,11 @@ subroutine ModelAdvance_slow(gcomp, rc) i1 = i - ilo + 1 j1 = j - jlo + 1 #ifdef CMEPS - rhoa (i,j,iblk) = dataPtr_rhoabot(i1,j1) ! import directly from mediator - if(dataPtr_pbot(i1,j1) .gt. 0.0) & - potT (i,j,iblk) = dataPtr_Tbot (i1,j1) * (100000./dataPtr_pbot(i1,j1))**0.286 ! Potential temperature (K) + if (dataPtr_Tbot(i1,j1) /= 0._ESMF_KIND_R8) rhoa(i,j,iblk) = dataPtr_pbot(i1,j1)/& + (287.058_ESMF_KIND_R8*(1._ESMF_KIND_R8+0.608_ESMF_KIND_R8*dataPtr_qbot(i1,j1))*dataPtr_Tbot(i1,j1)) + !rhoa (i,j,iblk) = dataPtr_rhoabot(i1,j1) ! import directly from mediator + if(dataPtr_pbot(i1,j1) .gt. 0.0_ESMF_KIND_R8) & + potT (i,j,iblk) = dataPtr_Tbot (i1,j1) * (100000._ESMF_KIND_R8/dataPtr_pbot(i1,j1))**0.286_ESMF_KIND_R8 ! Potential temperature (K) Tair (i,j,iblk) = dataPtr_Tbot (i1,j1) ! near surface temp, maybe lowest level (K) Qa (i,j,iblk) = dataPtr_qbot (i1,j1) ! near surface humidity, maybe lowest level (kg/kg) zlvl (i,j,iblk) = dataPtr_zlvl (i1,j1) ! height of the lowest level (m) @@ -1086,7 +1098,7 @@ subroutine ModelAdvance_slow(gcomp, rc) frain (i,j,iblk) = dataPtr_lprec (i1,j1) ! flux of rain (liquid only) fsnow (i,j,iblk) = dataPtr_fprec (i1,j1) ! flux of frozen precip sss (i,j,iblk) = dataPtr_sss (i1,j1) ! sea surface salinity (maybe for mushy layer) - sst (i,j,iblk) = dataPtr_sst (i1,j1) - 273.15 ! sea surface temp (may not be needed?) + sst (i,j,iblk) = dataPtr_sst (i1,j1) - Tffresh ! sea surface temp (may not be needed?) frzmlt (i,j,iblk) = dataPtr_fmpot (i1,j1) ! ! --- rotate these vectors from east/north to i/j --- uocn (i,j,iblk) = dataPtr_ocncz (i1,j1) @@ -1096,9 +1108,11 @@ subroutine ModelAdvance_slow(gcomp, rc) ss_tltx(i,j,iblk) = dataPtr_sssz (i1,j1) ss_tlty(i,j,iblk) = dataPtr_sssm (i1,j1) #else - rhoa (i,j,iblk) = dataPtr_rhoabot(i1,j1,iblk) ! import directly from mediator - if(dataPtr_pbot(i1,j1,iblk) .gt. 0.0) & - potT (i,j,iblk) = dataPtr_Tbot (i1,j1,iblk) * (100000./dataPtr_pbot(i1,j1,iblk))**0.286 ! Potential temperature (K) + if (dataPtr_Tbot(i1,j1,iblk) /= 0._ESMF_KIND_R8) rhoa(i,j,iblk) = dataPtr_pbot(i1,j1,iblk)/& + (287.058_ESMF_KIND_R8*(1._ESMF_KIND_R8+0.608_ESMF_KIND_R8*dataPtr_qbot(i1,j1,iblk))*dataPtr_Tbot(i1,j1,iblk)) + !rhoa (i,j,iblk) = dataPtr_rhoabot(i1,j1,iblk) ! import directly from mediator + if(dataPtr_pbot(i1,j1,iblk) .gt. 0.0_ESMF_KIND_R8) & + potT (i,j,iblk) = dataPtr_Tbot (i1,j1,iblk) * (100000._ESMF_KIND_R8/dataPtr_pbot(i1,j1,iblk))**0.286_ESMF_KIND_R8 ! Potential temperature (K) Tair (i,j,iblk) = dataPtr_Tbot (i1,j1,iblk) ! near surface temp, maybe lowest level (K) Qa (i,j,iblk) = dataPtr_qbot (i1,j1,iblk) ! near surface humidity, maybe lowest level (kg/kg) zlvl (i,j,iblk) = dataPtr_zlvl (i1,j1,iblk) ! height of the lowest level (m) @@ -1111,7 +1125,7 @@ subroutine ModelAdvance_slow(gcomp, rc) frain (i,j,iblk) = dataPtr_lprec (i1,j1,iblk) ! flux of rain (liquid only) fsnow (i,j,iblk) = dataPtr_fprec (i1,j1,iblk) ! flux of frozen precip sss (i,j,iblk) = dataPtr_sss (i1,j1,iblk) ! sea surface salinity (maybe for mushy layer) - sst (i,j,iblk) = dataPtr_sst (i1,j1,iblk) - 273.15 ! sea surface temp (may not be needed?) + sst (i,j,iblk) = dataPtr_sst (i1,j1,iblk) - Tffresh ! sea surface temp (may not be needed?) frzmlt (i,j,iblk) = dataPtr_fmpot (i1,j1,iblk) ! ! --- rotate these vectors from east/north to i/j --- uocn (i,j,iblk) = dataPtr_ocncz (i1,j1,iblk) @@ -1312,8 +1326,7 @@ subroutine ModelAdvance_slow(gcomp, rc) if(write_diagnostics) then call state_diagnose(exportState, 'cice_export', rc) - export_slice = export_slice + 1 - + fname = 'field_ice_export_'//trim(export_timestr)//'.nc' do i = 1,fldsFrIce_num fldname = fldsFrIce(i)%shortname call ESMF_StateGet(exportState, itemName=trim(fldname), itemType=itemType, rc=rc) @@ -1328,8 +1341,8 @@ subroutine ModelAdvance_slow(gcomp, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return #ifdef CMEPS - call ESMF_FieldWrite(lfield, fileName='field_ice_export_'//trim(fldname)//'.nc', & - timeslice=export_slice, overwrite=overwrite_timeslice, rc=rc) + call ESMF_FieldWrite(lfield, fileName=trim(fname), & + timeslice=1, overwrite=overwrite_timeslice, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -1350,8 +1363,8 @@ subroutine ModelAdvance_slow(gcomp, rc) fldptr2d(:,:) = fldptr(:,:,1) - call ESMF_FieldWrite(lfield2d, fileName='field_ice_export_'//trim(fldname)//'.nc', & - timeslice=export_slice, overwrite=overwrite_timeslice,rc=rc) + call ESMF_FieldWrite(lfield2d, fileName=trim(fname), & + timeslice=1, overwrite=overwrite_timeslice,rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -1372,24 +1385,27 @@ subroutine ModelAdvance_slow(gcomp, rc) ! Dump out all the cice internal fields to cross-examine with those connected with mediator ! This will help to determine roughly which fields can be hooked into cice - call dumpCICEInternal(ice_grid_i, import_slice, "inst_zonal_wind_height10m", "will provide", strax) - call dumpCICEInternal(ice_grid_i, import_slice, "inst_merid_wind_height10m", "will provide", stray) - call dumpCICEInternal(ice_grid_i, import_slice, "inst_pres_height_surface" , "will provide", zlvl) - call dumpCICEInternal(ice_grid_i, import_slice, "ocn_current_zonal", "will provide", uocn) - call dumpCICEInternal(ice_grid_i, import_slice, "ocn_current_merid", "will provide", vocn) - call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_slope_zonal", "will provide", ss_tltx) - call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_slope_merid", "will provide", ss_tlty) - call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_salinity", "will provide", sss) - call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_temperature", "will provide", sst) + !import_slice = import_slice + 1 + !call dumpCICEInternal(ice_grid_i, import_slice, "air_density_height_lowest", "will provide", rhoa) + !call dumpCICEInternal(ice_grid_i, import_slice, "inst_zonal_wind_height10m", "will provide", strax) + !call dumpCICEInternal(ice_grid_i, import_slice, "inst_merid_wind_height10m", "will provide", stray) + !call dumpCICEInternal(ice_grid_i, import_slice, "inst_pres_height_surface" , "will provide", zlvl) + !call dumpCICEInternal(ice_grid_i, import_slice, "ocn_current_zonal", "will provide", uocn) + !call dumpCICEInternal(ice_grid_i, import_slice, "ocn_current_merid", "will provide", vocn) + !call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_slope_zonal", "will provide", ss_tltx) + !call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_slope_merid", "will provide", ss_tlty) + !call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_salinity", "will provide", sss) + !call dumpCICEInternal(ice_grid_i, import_slice, "sea_surface_temperature", "will provide", sst) !--------- export fields from Sea Ice ------------- - call dumpCICEInternal(ice_grid_i, export_slice, "ice_fraction" , "will provide", aice) - call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_air_ice_zonal" , "will provide", strairxT) - call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_air_ice_merid" , "will provide", strairyT) - call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_ocn_ice_zonal" , "will provide", strocnxT) - call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_ocn_ice_merid" , "will provide", strocnyT) - call dumpCICEInternal(ice_grid_i, export_slice, "mean_sw_pen_to_ocn" , "will provide", fswthru) + !export_slice = export_slice + 1 + !call dumpCICEInternal(ice_grid_i, export_slice, "ice_fraction" , "will provide", aice) + !call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_air_ice_zonal" , "will provide", strairxT) + !call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_air_ice_merid" , "will provide", strairyT) + !call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_ocn_ice_zonal" , "will provide", strocnxT) + !call dumpCICEInternal(ice_grid_i, export_slice, "stress_on_ocn_ice_merid" , "will provide", strocnyT) + !call dumpCICEInternal(ice_grid_i, export_slice, "mean_sw_pen_to_ocn" , "will provide", fswthru) if(profile_memory) call ESMF_VMLogMemInfo("Leaving CICE Model_ADVANCE: ") end subroutine ModelAdvance_slow