diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index 919e7e2fd..388b848af 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -13,13 +13,14 @@ module icepack_atmo use icepack_kinds - use icepack_parameters, only: c0, c1, c2, c4, c5, c8, c10 - use icepack_parameters, only: c16, c20, p001, p01, p2, p4, p5, p75, puny - use icepack_parameters, only: senscoef, latncoef - use icepack_parameters, only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn - use icepack_parameters, only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit - use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow + use icepack_parameters, only: c0, c1, c2, c4, c5, c8, c10 + use icepack_parameters, only: c16, c20, p001, p01, p2, p4, p5, p75, puny + use icepack_parameters, only: senscoef, latncoef + use icepack_parameters, only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn + use icepack_parameters, only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit + use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow use icepack_parameters, only: atmbndy, calc_strair, formdrag + use icepack_parameters, only: icepack_chkoptargflag use icepack_tracers, only: n_iso use icepack_tracers, only: tr_iso use icepack_warnings, only: warnstr, icepack_warnings_add @@ -61,7 +62,6 @@ subroutine atmo_boundary_layer (sfctype, & Cdn_atm, & Cdn_atm_ratio_n, & Qa_iso, Qref_iso, & - iso_flag, & uvel, vvel, & Uref, zlvs ) @@ -103,13 +103,10 @@ subroutine atmo_boundary_layer (sfctype, & shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat - logical (kind=log_kind), intent(in), optional :: & - iso_flag ! flag to trigger iso calculations - - real (kind=dbl_kind), intent(in), optional, dimension(:) :: & + real (kind=dbl_kind), intent(in), dimension(:), optional :: & Qa_iso ! specific isotopic humidity (kg/kg) - real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & + real (kind=dbl_kind), intent(inout), dimension(:), optional :: & Qref_iso ! reference specific isotopic humidity (kg/kg) real (kind=dbl_kind), intent(in) :: & @@ -167,16 +164,8 @@ subroutine atmo_boundary_layer (sfctype, & real (kind=dbl_kind), parameter :: & zTrf = c2 ! reference height for air temp (m) - logical (kind=log_kind) :: & - l_iso_flag ! local flag to trigger iso calculations - character(len=*),parameter :: subname='(atmo_boundary_layer)' - l_iso_flag = .false. - if (present(iso_flag)) then - l_iso_flag = iso_flag - endif - al2 = log(zref/zTrf) !------------------------------------------------------------ @@ -389,21 +378,13 @@ subroutine atmo_boundary_layer (sfctype, & Uref = vmag * rd / rdn endif - if (l_iso_flag) then - if (present(Qref_iso) .and. present(Qa_iso)) then + if (tr_iso .and. sfctype(1:3)=='ice') then Qref_iso(:) = c0 - if (tr_iso) then - do n = 1, n_iso - ratio = c0 - if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2) - Qref_iso(n) = Qa_iso(n) - ratio*delq*fac - enddo - endif - else - call icepack_warnings_add(subname//' l_iso_flag true but optional arrays missing') - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - return - endif + do n = 1, n_iso + ratio = c0 + if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2) + Qref_iso(n) = Qa_iso(n) - ratio*delq*fac + enddo endif end subroutine atmo_boundary_layer @@ -895,18 +876,18 @@ subroutine icepack_atm_boundary(sfctype, & shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat - real (kind=dbl_kind), intent(in), optional, dimension(:) :: & + real (kind=dbl_kind), intent(in), dimension(:), optional :: & Qa_iso ! specific isotopic humidity (kg/kg) - real (kind=dbl_kind), intent(inout), optional, dimension(:) :: & + real (kind=dbl_kind), intent(inout), dimension(:), optional :: & Qref_iso ! reference specific isotopic humidity (kg/kg) - real (kind=dbl_kind), optional, intent(in) :: & + real (kind=dbl_kind), intent(in), optional :: & uvel , & ! x-direction ice speed (m/s) vvel , & ! y-direction ice speed (m/s) zlvs ! atm level height for scalars (if different than zlvl) (m) - real (kind=dbl_kind), optional, intent(out) :: & + real (kind=dbl_kind), intent(out), optional :: & Uref ! reference height wind speed (m/s) !autodocument_end @@ -916,14 +897,28 @@ subroutine icepack_atm_boundary(sfctype, & real (kind=dbl_kind) :: & l_uvel, l_vvel, l_Uref - real (kind=dbl_kind), dimension(:), allocatable :: & - l_Qa_iso, l_Qref_iso ! local copies of Qa_iso, Qref_iso - - logical (kind=log_kind) :: & - iso_flag ! flag to turn on iso calcs in other subroutines + logical (kind=log_kind), save :: & + first_call_ice = .true. ! first call flag character(len=*),parameter :: subname='(icepack_atm_boundary)' + !------------------------------------------------------------ + ! Check optional arguments + ! Need separate first_call flags for 'ice' and 'ocn' sfctype + !------------------------------------------------------------ + + if (sfctype == 'ice') then + if (icepack_chkoptargflag(first_call_ice)) then + if (tr_iso) then + if (.not.(present(Qa_iso).and.present(Qref_iso))) then + call icepack_warnings_add(subname//' error in fiso_ocn argument, tr_iso=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif + endif + endif + l_uvel = c0 l_vvel = c0 l_Uref = c0 @@ -933,19 +928,6 @@ subroutine icepack_atm_boundary(sfctype, & if (present(vvel)) then l_vvel = vvel endif - if (present(Qa_iso) .and. present(Qref_iso)) then - iso_flag = .true. - allocate(l_Qa_iso(size(Qa_iso,dim=1))) - allocate(l_Qref_iso(size(Qref_iso,dim=1))) - l_Qa_iso = Qa_iso - l_Qref_iso = Qref_iso - else - iso_flag = .false. - allocate(l_Qa_iso(1)) - allocate(l_Qref_iso(1)) - l_Qa_iso = c0 - l_Qref_iso = c0 - endif Cdn_atm_ratio_n = c1 @@ -972,11 +954,10 @@ subroutine icepack_atm_boundary(sfctype, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - iso_flag = iso_flag, & - Qa_iso=l_Qa_iso, & - Qref_iso=l_Qref_iso, & - uvel=l_uvel, vvel=l_vvel, & - Uref=l_Uref, zlvs=zlvs) + Qa_iso=Qa_iso, & + Qref_iso=Qref_iso, & + uvel=l_uvel, vvel=l_vvel,& + Uref=l_Uref, zlvs=zlvs ) if (icepack_warnings_aborted(subname)) return endif ! atmbndy @@ -984,12 +965,10 @@ subroutine icepack_atm_boundary(sfctype, & Uref = l_Uref endif - if (present(Qref_iso)) then - Qref_iso = l_Qref_iso + if (sfctype == 'ice') then + first_call_ice = .false. endif - deallocate(l_Qa_iso,l_Qref_iso) - end subroutine icepack_atm_boundary !======================================================================= diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90 index 50bbd72d4..54e33b8ad 100644 --- a/columnphysics/icepack_flux.F90 +++ b/columnphysics/icepack_flux.F90 @@ -89,10 +89,6 @@ subroutine merge_fluxes (aicen, & fsaltn , & ! salt flux to ocean (kg/m2/s) fhocnn , & ! actual ocn/ice heat flx (W/m**2) fswthrun, & ! sw radiation through ice bot (W/m**2) - fswthrun_vdr, & ! vis dir sw radiation through ice bot (W/m**2) - fswthrun_vdf, & ! vis dif sw radiation through ice bot (W/m**2) - fswthrun_idr, & ! nir dir sw radiation through ice bot (W/m**2) - fswthrun_idf, & ! nir dif sw radiation through ice bot (W/m**2) melttn , & ! top ice melt (m) meltbn , & ! bottom ice melt (m) meltsn , & ! snow melt (m) @@ -102,6 +98,10 @@ subroutine merge_fluxes (aicen, & snoicen ! snow-ice growth (m) real (kind=dbl_kind), optional, intent(in):: & + fswthrun_vdr, & ! vis dir sw radiation through ice bot (W/m**2) + fswthrun_vdf, & ! vis dif sw radiation through ice bot (W/m**2) + fswthrun_idr, & ! nir dir sw radiation through ice bot (W/m**2) + fswthrun_idf, & ! nir dif sw radiation through ice bot (W/m**2) Urefn ! air speed reference level (m/s) ! cumulative fluxes @@ -129,7 +129,6 @@ subroutine merge_fluxes (aicen, & meltb , & ! bottom ice melt (m) melts , & ! snow melt (m) meltsliq, & ! mass of snow melt (kg/m^2) - dsnow , & ! change in snow depth (m) congel , & ! congelation ice growth (m) snoice ! snow-ice growth (m) @@ -139,18 +138,19 @@ subroutine merge_fluxes (aicen, & fswthru_idr , & ! nir dir sw radiation through ice bot (W/m**2) fswthru_idf ! nir dif sw radiation through ice bot (W/m**2) - real (kind=dbl_kind), optional, intent(inout):: & + real (kind=dbl_kind), intent(inout), optional :: & + dsnow, & ! change in snow depth (m) Uref ! air speed reference level (m/s) - real (kind=dbl_kind), optional, dimension(:), intent(inout):: & - Qref_iso, & ! isotope air sp hum reference level (kg/kg) - fiso_ocn, & ! isotope fluxes to ocean (kg/m2/s) - fiso_evap ! isotope evaporation (kg/m2/s) + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + Qref_iso, & ! isotope air sp hum ref level (kg/kg) + fiso_ocn, & ! isotope fluxes to ocean (kg/m2/s) + fiso_evap ! isotope evaporation (kg/m2/s) - real (kind=dbl_kind), optional, dimension(:), intent(in):: & - Qrefn_iso, & ! isotope air sp hum reference level (kg/kg) - fiso_ocnn, & ! isotope fluxes to ocean (kg/m2/s) - fiso_evapn ! isotope evaporation (kg/m2/s) + real (kind=dbl_kind), dimension(:), intent(in), optional :: & + Qrefn_iso, & ! isotope air sp hum ref level (kg/kg) + fiso_ocnn, & ! isotope fluxes to ocean (kg/m2/s) + fiso_evapn ! isotope evaporation (kg/m2/s) character(len=*),parameter :: subname='(merge_fluxes)' @@ -220,7 +220,9 @@ subroutine merge_fluxes (aicen, & if (tr_snow) then meltsliq = meltsliq + meltsliqn * aicen endif - dsnow = dsnow + dsnown * aicen + if (present(dsnow)) then + dsnow = dsnow + dsnown * aicen + endif congel = congel + congeln * aicen snoice = snoice + snoicen * aicen diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90 index 495d0bde2..1db640556 100644 --- a/columnphysics/icepack_isotope.F90 +++ b/columnphysics/icepack_isotope.F90 @@ -83,21 +83,24 @@ subroutine update_isotope (dt, & aicen, & ! ice area aice_old, & ! beginning values vice_old, & - vsno_old, & - HDO_ocn, & ! - H2_16O_ocn, & ! - H2_18O_ocn ! - - real (kind=dbl_kind), dimension(:), intent(in) :: & - fiso_atm, & ! isotopic snowfall (kg/m^2/s of water) - Qref_iso ! isotope reference humidity + vsno_old real (kind=dbl_kind), dimension(:), intent(inout) :: & fiso_ocnn, & ! isotopic freshwater (kg/m^2/s) fiso_evapn ! evaporative water flux (kg/m^2/s) real (kind=dbl_kind), dimension(:), intent(inout) :: & - isosno, isoice ! mass of isotopes (kg) + isosno, & ! mass of isotopes (kg) + isoice + + real (kind=dbl_kind), dimension(:), intent(in) :: & + fiso_atm, & ! isotopic snowfall (kg/m^2/s of water) + Qref_iso ! isotope reference humidity + + real (kind=dbl_kind), intent(in) :: & + HDO_ocn, & ! ocean concentration of HDO (kg/kg) + H2_16O_ocn, & ! ocean concentration of H2_16O (kg/kg) + H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) ! local variables @@ -197,6 +200,7 @@ subroutine update_isotope (dt, & if (congel > c0) then do k = 1,n_iso + work = c0 if (k == 1) then alpha = isoice_alpha(congel/dt,'HDO',isotope_frac_method) work = alpha*HDO_ocn*rhoi*congel*aicen diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index b32954ca4..00f0f768f 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -833,7 +833,7 @@ subroutine cleanup_itd (dt, ntrcr, & real (kind=dbl_kind), dimension (:), intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) - real (kind=dbl_kind), dimension (:), intent(inout) :: & + real (kind=dbl_kind), dimension (:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) logical (kind=log_kind), intent(in), optional :: & diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 5ef805744..4dc87cb28 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -34,9 +34,10 @@ module icepack_mechred use icepack_kinds - use icepack_parameters, only: c0, c1, c2, c10, c25, Cf, Cp, Pstar, Cstar - use icepack_parameters, only: p05, p15, p25, p333, p5 - use icepack_parameters, only: puny, Lfresh, rhoi, rhos + use icepack_parameters, only: c0, c1, c2, c10, c25, Cf, Cp, Pstar, Cstar + use icepack_parameters, only: p05, p15, p25, p333, p5 + use icepack_parameters, only: puny, Lfresh, rhoi, rhos + use icepack_parameters, only: icepack_chkoptargflag use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg use icepack_parameters, only: conserv_check @@ -188,7 +189,7 @@ subroutine ridge_ice (dt, ndtd, & real (kind=dbl_kind), dimension(:), intent(inout), optional :: & faero_ocn ! aerosol flux to ocean (kg/m^2/s) - real (kind=dbl_kind), dimension(:), intent(inout) :: & + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) ! local variables @@ -1804,7 +1805,7 @@ subroutine icepack_step_ridge (dt, ndtd, & faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) flux_bio ! all bio fluxes to ocean - real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) real (kind=dbl_kind), dimension(:,:), intent(inout) :: & @@ -1825,17 +1826,32 @@ subroutine icepack_step_ridge (dt, ndtd, & real (kind=dbl_kind) :: & dtt ! thermo time step - real (kind=dbl_kind), dimension(:), allocatable :: & - l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s) - real (kind=dbl_kind) :: & l_closing ! local rate of closing due to divergence/shear (1/s) logical (kind=log_kind) :: & l_closing_flag ! flag if closing is passed + logical (kind=log_kind), save :: & + first_call = .true. ! first call flag + character(len=*),parameter :: subname='(icepack_step_ridge)' + !----------------------------------------------------------------- + ! Check optional arguments + !----------------------------------------------------------------- + + if (icepack_chkoptargflag(first_call)) then + if (tr_iso) then + if (.not.(present(fiso_ocn))) then + call icepack_warnings_add(subname//' error in fiso_ocn argument, tr_iso=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif + endif + + !----------------------------------------------------------------- ! Identify ice-ocean cells. ! Note: We can not limit the loop here using aice>puny because @@ -1843,15 +1859,6 @@ subroutine icepack_step_ridge (dt, ndtd, & ! it may be out of whack, which the ridging helps fix).-ECH !----------------------------------------------------------------- - if (present(fiso_ocn)) then - allocate(l_fiso_ocn(size(fiso_ocn))) - l_fiso_ocn = fiso_ocn - else - ! check tr_iso = true ??? - allocate(l_fiso_ocn(1)) - l_fiso_ocn = c0 - endif - if (present(closing)) then l_closing_flag = .true. l_closing = closing @@ -1879,7 +1886,7 @@ subroutine icepack_step_ridge (dt, ndtd, & dvirdgdt, opening, & fpond, & fresh, fhocn, & - faero_ocn, l_fiso_ocn, & + faero_ocn, fiso_ocn, & aparticn, krdgn, & aredistn, vredistn, & dardg1ndt, dardg2ndt, & @@ -1910,13 +1917,12 @@ subroutine icepack_step_ridge (dt, ndtd, & n_trcr_strata, nt_strata, & fpond, fresh, & fsalt, fhocn, & - faero_ocn, l_fiso_ocn, & + faero_ocn, fiso_ocn, & fzsal, & flux_bio) if (icepack_warnings_aborted(subname)) return - if (present(fiso_ocn)) fiso_ocn = l_fiso_ocn - deallocate(l_fiso_ocn) + first_call = .false. end subroutine icepack_step_ridge diff --git a/columnphysics/icepack_meltpond_lvl.F90 b/columnphysics/icepack_meltpond_lvl.F90 index 47ca49c98..5467db4b9 100644 --- a/columnphysics/icepack_meltpond_lvl.F90 +++ b/columnphysics/icepack_meltpond_lvl.F90 @@ -48,6 +48,7 @@ subroutine compute_ponds_lvl(dt, nilyr, & integer (kind=int_kind), intent(in) :: & nilyr, & ! number of ice layers ktherm ! type of thermodynamics (-1 none, 1 BL99, 2 mushy) + real (kind=dbl_kind), intent(in) :: & dt, & ! time step (s) hi_min, & ! minimum ice thickness allowed for thermo (m) @@ -56,8 +57,7 @@ subroutine compute_ponds_lvl(dt, nilyr, & character (len=char_len), intent(in) :: & frzpnd ! pond refreezing parameterization - real (kind=dbl_kind), & - intent(in) :: & + real (kind=dbl_kind), intent(in) :: & Tsfcn, & ! surface temperature (C) alvl, & ! fraction of level ice rfrac, & ! water fraction retained for melt ponds @@ -71,20 +71,17 @@ subroutine compute_ponds_lvl(dt, nilyr, & vsnon, & ! snow volume (m) meltsliqn ! liquid contribution to meltponds in dt (kg/m^2) - real (kind=dbl_kind), & - intent(inout) :: & + real (kind=dbl_kind), intent(inout) :: & apnd, hpnd, ipnd real (kind=dbl_kind), dimension (:), intent(in) :: & qicen, & ! ice layer enthalpy (J m-3) sicen ! salinity (ppt) - real (kind=dbl_kind), & - intent(in) :: & + real (kind=dbl_kind), intent(in) :: & dhs ! depth difference for snow on sea ice and pond ice - real (kind=dbl_kind), & - intent(out) :: & + real (kind=dbl_kind), intent(out) :: & ffrac ! fraction of fsurfn over pond used to melt ipond ! local temporary variables diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 0c979772e..70b147474 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -17,6 +17,7 @@ module icepack_parameters public :: icepack_query_parameters public :: icepack_write_parameters public :: icepack_recompute_constants + public :: icepack_chkoptargflag !----------------------------------------------------------------- ! control options @@ -1985,6 +1986,22 @@ end subroutine icepack_recompute_constants !======================================================================= + function icepack_chkoptargflag(first_call) result(chkoptargflag) + + logical(kind=log_kind), intent(in) :: first_call + + logical(kind=log_kind) :: chkoptargflag + + character(len=*),parameter :: subname='(icepack_chkoptargflag)' + + chkoptargflag = & + (argcheck == 'always' .or. (argcheck == 'first' .and. first_call)) + + end function icepack_chkoptargflag + +!======================================================================= + + end module icepack_parameters !======================================================================= diff --git a/columnphysics/icepack_snow.F90 b/columnphysics/icepack_snow.F90 index a2a481743..20b9a980a 100644 --- a/columnphysics/icepack_snow.F90 +++ b/columnphysics/icepack_snow.F90 @@ -919,7 +919,7 @@ subroutine snow_dry_metamorph (nslyr,nilyr, dt, rsnw, drsnw_dry, zqsn, & dzi, & ! ice layer thickness (m) dz ! dzs + dzi (m) - logical (kind=log_kind) :: & + logical (kind=log_kind), save :: & first_call = .true. ! first call flag character (char_len) :: & @@ -1194,7 +1194,6 @@ subroutine drain_snow (nslyr, vsnon, aicen, & else sliq = meltsliq ! computed in thickness_changes endif - meltsliq = meltsliq if (use_smliq_pnd) meltsliq = sliq end subroutine drain_snow diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 743d96d0f..e4af94abc 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -28,6 +28,7 @@ module icepack_therm_itd use icepack_parameters, only: kitd, ktherm use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin use icepack_parameters, only: saltflux_option + use icepack_parameters, only: icepack_chkoptargflag use icepack_tracers, only: ntrcr, nbtrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice @@ -1730,14 +1731,14 @@ subroutine add_new_ice (ncat, nilyr, & enddo endif - frazil_conc = c0 if (tr_iso .and. vtmp > puny) then do it=1,n_iso - if (it==1) & + frazil_conc = c0 + if (it==1) & frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn - if (it==2) & + if (it==2) & frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn - if (it==3) & + if (it==3) & frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn ! dilution and uptake in the ice @@ -1852,15 +1853,15 @@ subroutine add_new_ice (ncat, nilyr, & enddo endif - frazil_conc = c0 if (tr_iso) then do it=1,n_iso - if (it==1) & - frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn - if (it==2) & - frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn - if (it==3) & - frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn + frazil_conc = c0 + if (it==1) & + frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn + if (it==2) & + frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn + if (it==3) & + frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn trcrn(nt_isoice+it-1,1) & = (trcrn(nt_isoice+it-1,1)*vice1) & @@ -2093,10 +2094,10 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & yday ! day of year ! water isotopes - real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) - real (kind=dbl_kind), optional, intent(in) :: & + real (kind=dbl_kind), intent(in), optional :: & HDO_ocn , & ! ocean concentration of HDO (kg/kg) H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) H2_18O_ocn ! ocean concentration of H2_18O (kg/kg) @@ -2104,14 +2105,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & ! local variables - ! water isotopes - real (kind=dbl_kind), dimension(:), allocatable :: & - l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s) - - real (kind=dbl_kind) :: & - l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg) - l_H2_16O_ocn , & ! local ocean concentration of H2_16O (kg/kg) - l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg) + logical (kind=log_kind), save :: & + first_call = .true. ! first call flag character(len=*),parameter :: subname='(icepack_step_therm2)' @@ -2119,19 +2114,18 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & ! Check optional arguments and set local values !----------------------------------------------------------------- - if (present(fiso_ocn)) then - allocate(l_fiso_ocn(size(fiso_ocn))) - l_fiso_ocn(:) = fiso_ocn(:) - else - allocate(l_fiso_ocn(1)) - l_fiso_ocn(:) = c0 + if (icepack_chkoptargflag(first_call)) then + if (tr_iso) then + if (.not.(present(fiso_ocn) .and. & + present(HDO_ocn) .and. & + present(H2_16O_ocn) .and. & + present(H2_18O_ocn))) then + call icepack_warnings_add(subname//' error in iso arguments, tr_iso=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif endif - l_HDO_ocn = c0 - l_H2_16O_ocn = c0 - l_H2_18O_ocn = c0 - if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn - if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn - if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn !----------------------------------------------------------------- ! Let rain drain through to the ocean. @@ -2208,9 +2202,9 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & - frazil_diag, l_fiso_ocn, & - l_HDO_ocn, l_H2_16O_ocn, & - l_H2_18O_ocn, & + frazil_diag, fiso_ocn, & + HDO_ocn, H2_16O_ocn, & + H2_18O_ocn, & wave_sig_ht, & wave_spectrum, & wavefreq, dwavefreq, & @@ -2228,7 +2222,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & n_aero, fpond, & fresh, fsalt, & fhocn, faero_ocn, & - l_fiso_ocn, & + fiso_ocn, & rside, meltl, & fside, sss, & aicen, vicen, & @@ -2280,14 +2274,11 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & n_trcr_strata, nt_strata, & fpond, fresh, & fsalt, fhocn, & - faero_ocn, l_fiso_ocn, & + faero_ocn, fiso_ocn, & fzsal, flux_bio) if (icepack_warnings_aborted(subname)) return - if (present(fiso_ocn)) then - fiso_ocn(:) = l_fiso_ocn(:) - endif - deallocate(l_fiso_ocn) + first_call = .false. end subroutine icepack_step_therm2 diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 618fa0d1c..ff9b775f3 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -28,6 +28,7 @@ module icepack_therm_vertical use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist use icepack_parameters, only: saltflux_option + use icepack_parameters, only: icepack_chkoptargflag use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso use icepack_tracers, only: tr_pond_lvl, tr_pond_topo @@ -190,8 +191,7 @@ subroutine thermo_vertical (nilyr, nslyr, & fhocnn ! net heat flux to ocean (W/m^2) ! diagnostic fields - real (kind=dbl_kind), & - intent(inout):: & + real (kind=dbl_kind), intent(inout):: & Tsnice , & ! snow ice interface temperature (deg C) meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) @@ -2171,15 +2171,15 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & zlvs) integer (kind=int_kind), intent(in) :: & - ncat , & ! number of thickness categories - nilyr , & ! number of ice layers - nslyr ! number of snow layers + ncat , & ! number of thickness categories + nilyr , & ! number of ice layers + nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & dt , & ! time step - uvel , & ! x-component of velocity (m/s) - vvel , & ! y-component of velocity (m/s) - strax , & ! wind stress components (N/m^2) + uvel , & ! x-component of velocity (m/s) + vvel , & ! y-component of velocity (m/s) + strax , & ! wind stress components (N/m^2) stray , & ! yday ! day of year @@ -2192,42 +2192,42 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration - vice , & ! volume per unit area of ice (m) - vsno , & ! volume per unit area of snow (m) + vice , & ! volume per unit area of ice (m) + vsno , & ! volume per unit area of snow (m) zlvl , & ! atm level height for momentum (and scalars if zlvs is not present) (m) - uatm , & ! wind velocity components (m/s) - vatm , & - wind , & ! wind speed (m/s) - potT , & ! air potential temperature (K) - Tair , & ! air temperature (K) - Qa , & ! specific humidity (kg/kg) - rhoa , & ! air density (kg/m^3) - frain , & ! rainfall rate (kg/m^2 s) - fsnow , & ! snowfall rate (kg/m^2 s) - fpond , & ! fresh water flux to ponds (kg/m^2/s) - fresh , & ! fresh water flux to ocean (kg/m^2/s) - fsalt , & ! salt flux to ocean (kg/m^2/s) - fhocn , & ! net heat flux to ocean (W/m^2) - fswthru , & ! shortwave penetrating to ocean (W/m^2) + uatm , & ! wind velocity components (m/s) + vatm , & ! (m/s) + wind , & ! wind speed (m/s) + potT , & ! air potential temperature (K) + Tair , & ! air temperature (K) + Qa , & ! specific humidity (kg/kg) + rhoa , & ! air density (kg/m^3) + frain , & ! rainfall rate (kg/m^2 s) + fsnow , & ! snowfall rate (kg/m^2 s) + fpond , & ! fresh water flux to ponds (kg/m^2/s) + fresh , & ! fresh water flux to ocean (kg/m^2/s) + fsalt , & ! salt flux to ocean (kg/m^2/s) + fhocn , & ! net heat flux to ocean (W/m^2) + fswthru , & ! shortwave penetrating to ocean (W/m^2) fsurf , & ! net surface heat flux (excluding fcondtop)(W/m^2) fcondtop , & ! top surface conductive flux (W/m^2) fcondbot , & ! bottom surface conductive flux (W/m^2) - fsens , & ! sensible heat flux (W/m^2) - flat , & ! latent heat flux (W/m^2) + fsens , & ! sensible heat flux (W/m^2) + flat , & ! latent heat flux (W/m^2) fswabs , & ! shortwave flux absorbed in ice and ocean (W/m^2) - flw , & ! incoming longwave radiation (W/m^2) - flwout , & ! outgoing longwave radiation (W/m^2) - evap , & ! evaporative water flux (kg/m^2/s) - evaps , & ! evaporative water flux over snow (kg/m^2/s) + flw , & ! incoming longwave radiation (W/m^2) + flwout , & ! outgoing longwave radiation (W/m^2) + evap , & ! evaporative water flux (kg/m^2/s) + evaps , & ! evaporative water flux over snow(kg/m^2/s) evapi , & ! evaporative water flux over ice (kg/m^2/s) congel , & ! basal ice growth (m/step-->cm/day) snoice , & ! snow-ice formation (m/step-->cm/day) - Tref , & ! 2m atm reference temperature (K) - Qref , & ! 2m atm reference spec humidity (kg/kg) - Uref , & ! 10m atm reference wind speed (m/s) + Tref , & ! 2m atm reference temperature (K) + Qref , & ! 2m atm reference spec humidity (kg/kg) + Uref , & ! 10m atm reference wind speed (m/s) Cdn_atm , & ! atm drag coefficient Cdn_ocn , & ! ocn drag coefficient - hfreebd , & ! freeboard (m) + hfreebd , & ! freeboard (m) hdraft , & ! draft of ice + snow column (Stoessel1993) hridge , & ! ridge height distrdg , & ! distance between ridges @@ -2248,14 +2248,14 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & strocnxT , & ! ice-ocean stress, x-direction strocnyT , & ! ice-ocean stress, y-direction fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) - frzmlt , & ! freezing/melting potential (W/m^2) + frzmlt , & ! freezing/melting potential (W/m^2) rside , & ! fraction of ice that melts laterally - fside , & ! lateral heat flux (W/m^2) - sst , & ! sea surface temperature (C) - Tf , & ! freezing temperature (C) - Tbot , & ! ice bottom surface temperature (deg C) - Tsnice , & ! snow ice interface temperature (deg C) - sss , & ! sea surface salinity (ppt) + fside , & ! lateral heat flux (W/m^2) + sst , & ! sea surface temperature (C) + Tf , & ! freezing temperature (C) + Tbot , & ! ice bottom surface temperature (deg C) + Tsnice , & ! snow ice interface temperature (deg C) + sss , & ! sea surface salinity (ppt) meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) meltb , & ! basal ice melt (m/step-->cm/day) @@ -2268,86 +2268,88 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2) fswthru_idf , & ! nir dif shortwave penetrating to ocean (W/m^2) dsnow , & ! change in snow depth (m/step-->cm/day) - meltsliq , & ! mass of snow melt (kg/m^2) - fsloss ! rate of snow loss to leads (kg/m^2/s) - - real (kind=dbl_kind), dimension(:), optional, intent(inout) :: & - Qa_iso , & ! isotope specific humidity (kg/kg) - Qref_iso , & ! isotope 2m atm reference spec humidity (kg/kg) - fiso_atm , & ! isotope deposition rate (kg/m^2 s) - fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) - fiso_evap , & ! isotope evaporation (kg/m^2/s) - meltsliqn ! mass of snow melt (kg/m^2) - - real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: & - rsnwn , & ! snow grain radius (10^-6 m) - smicen , & ! tracer for mass of ice in snow (kg/m^3) - smliqn ! tracer for mass of liquid in snow (kg/m^3) - - real (kind=dbl_kind), optional, intent(in) :: & - HDO_ocn , & ! ocean concentration of HDO (kg/kg) - H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) - H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) + meltsliq , & ! mass of snow melt (kg/m^2) + fsloss ! rate of snow loss to leads (kg/m^2/s) + + real (kind=dbl_kind), dimension(:), intent(inout), optional :: & + Qa_iso , & ! isotope specific humidity (kg/kg) + Qref_iso , & ! isotope 2m atm ref spec humidity (kg/kg) + fiso_atm , & ! isotope deposition rate (kg/m^2 s) + fiso_ocn , & ! isotope flux to ocean (kg/m^2/s) + fiso_evap , & ! isotope evaporation (kg/m^2/s) + meltsliqn ! mass of snow melt (kg/m^2) + + real (kind=dbl_kind), dimension(:,:), intent(inout), optional :: & + rsnwn , & ! snow grain radius (10^-6 m) + smicen , & ! tracer for mass of ice in snow (kg/m^3) + smliqn ! tracer for mass of liq in snow (kg/m^3) + + real (kind=dbl_kind), intent(in), optional :: & + HDO_ocn , & ! ocean concentration of HDO (kg/kg) + H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg) + H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) zlvs ! atm level height for scalars (if different than zlvl) (m) real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice - vicen_init , & ! volume per unit area of ice (m) - vsnon_init , & ! volume per unit area of snow (m) + vicen_init , & ! volume per unit area of ice (m) + vsnon_init , & ! volume per unit area of snow (m) aicen , & ! concentration of ice - vicen , & ! volume per unit area of ice (m) - vsnon , & ! volume per unit area of snow (m) + vicen , & ! volume per unit area of ice (m) + vsnon , & ! volume per unit area of snow (m) Tsfc , & ! ice/snow surface temperature, Tsfcn alvl , & ! level ice area fraction vlvl , & ! level ice volume fraction apnd , & ! melt pond area fraction - hpnd , & ! melt pond depth (m) - ipnd , & ! melt pond refrozen lid thickness (m) + hpnd , & ! melt pond depth (m) + ipnd , & ! melt pond refrozen lid thickness (m) iage , & ! volume-weighted ice age FY , & ! area-weighted first-year ice area fsurfn , & ! net flux to top surface, excluding fcondtop - fcondtopn , & ! downward cond flux at top surface (W m-2) + fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) - flatn , & ! latent heat flux (W m-2) - fsensn , & ! sensible heat flux (W m-2) + flatn , & ! latent heat flux (W m-2) + fsensn , & ! sensible heat flux (W m-2) fsurfn_f , & ! net flux to top surface, excluding fcondtop - fcondtopn_f , & ! downward cond flux at top surface (W m-2) - flatn_f , & ! latent heat flux (W m-2) - fsensn_f , & ! sensible heat flux (W m-2) - fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) - fswthrun , & ! SW through ice to ocean (W/m^2) + fcondtopn_f , & ! downward cond flux at top surface (W m-2) + flatn_f , & ! latent heat flux (W m-2) + fsensn_f , & ! sensible heat flux (W m-2) + fswsfcn , & ! SW absorbed at ice/snow surface (W m-2) fswintn , & ! SW absorbed in ice interior, below surface (W m-2) - faero_atm , & ! aerosol deposition rate (kg/m^2 s) - faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) + faero_atm , & ! aerosol deposition rate (kg/m^2 s) + faero_ocn , & ! aerosol flux to ocean (kg/m^2/s) dhsn , & ! depth difference for snow on sea ice and pond ice ffracn , & ! fraction of fsurfn used to melt ipond - meltsn , & ! snow melt (m) - melttn , & ! top ice melt (m) - meltbn , & ! bottom ice melt (m) - congeln , & ! congelation ice growth (m) - snoicen , & ! snow-ice growth (m) + meltsn , & ! snow melt (m) + melttn , & ! top ice melt (m) + meltbn , & ! bottom ice melt (m) + congeln , & ! congelation ice growth (m) + snoicen , & ! snow-ice growth (m) dsnown ! change in snow thickness (m/step-->cm/day) - real (kind=dbl_kind), optional, dimension(:), intent(inout) :: & - fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) - fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) - fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) - fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + real (kind=dbl_kind), dimension(:), intent(in) :: & + fswthrun ! SW through ice to ocean (W/m^2) + + real (kind=dbl_kind), dimension(:), intent(in), optional :: & + fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) + fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) + fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) + fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) real (kind=dbl_kind), dimension(:,:), intent(inout) :: & - zqsn , & ! snow layer enthalpy (J m-3) - zqin , & ! ice layer enthalpy (J m-3) + zqsn , & ! snow layer enthalpy (J m-3) + zqin , & ! ice layer enthalpy (J m-3) zSin , & ! internal ice layer salinities Sswabsn , & ! SW radiation absorbed in snow layers (W m-2) - Iswabsn ! SW radiation absorbed in ice layers (W m-2) + Iswabsn ! SW radiation absorbed in ice layers (W m-2) real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: & - aerosno , & ! snow aerosol tracer (kg/m^2) - aeroice ! ice aerosol tracer (kg/m^2) + aerosno , & ! snow aerosol tracer (kg/m^2) + aeroice ! ice aerosol tracer (kg/m^2) - real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: & - isosno , & ! snow isotope tracer (kg/m^2) - isoice ! ice isotope tracer (kg/m^2) + real (kind=dbl_kind), dimension(:,:), intent(inout), optional :: & + isosno , & ! snow isotope tracer (kg/m^2) + isoice ! ice isotope tracer (kg/m^2) !autodocument_end ! local variables @@ -2365,7 +2367,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fswabsn , & ! shortwave absorbed by ice (W/m^2) flwoutn , & ! upward LW at surface (W/m^2) evapn , & ! flux of vapor, atmos to ice (kg m-2 s-1) - evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1) + evapsn , & ! flux of vapor, atmos to ice over snow (kg m-2 s-1) evapin , & ! flux of vapor, atmos to ice over ice (kg m-2 s-1) freshn , & ! flux of water, ice to ocean (kg/m^2/s) fsaltn , & ! flux of salt, ice to ocean (kg/m^2/s) @@ -2381,154 +2383,79 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & rfrac ! water fraction retained for melt ponds real (kind=dbl_kind), dimension(nslyr,ncat) :: & - massicen , & ! mass of ice in snow (kg/m^2) - massliqn ! mass of liquid in snow (kg/m^2) + massicen , & ! mass of ice in snow (kg/m^2) + massliqn ! mass of liquid in snow (kg/m^2) real (kind=dbl_kind), dimension(n_iso) :: & - Qrefn_iso , & ! isotope air sp hum reference level (kg/kg) - fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s) - fiso_evapn ! isotope evaporation (kg/m^2/s) + Qrefn_iso , & ! isotope air sp hum reference level (kg/kg) + fiso_ocnn , & ! isotope flux to ocean (kg/m^2/s) + fiso_evapn ! isotope evaporation (kg/m^2/s) - real (kind=dbl_kind), allocatable, dimension(:,:) :: & - l_isosno , & ! local snow isotope tracer (kg/m^2) - l_isoice ! local ice isotope tracer (kg/m^2) + real (kind=dbl_kind), dimension(nslyr) :: & + rsnw , & ! snow grain radius (10^-6 m) + smice , & ! tracer for mass of ice in snow (kg/m^3) + smliq ! tracer for mass of liquid in snow (kg/m^3) real (kind=dbl_kind), allocatable, dimension(:) :: & - l_Qa_iso , & ! local isotope specific humidity (kg/kg) - l_Qref_iso , & ! local isotope 2m atm reference spec humidity (kg/kg) - l_fiso_atm , & ! local isotope deposition rate (kg/m^2 s) - l_fiso_ocn , & ! local isotope flux to ocean (kg/m^2/s) - l_fiso_evap , & ! local isotope evaporation (kg/m^2/s) - l_meltsliqn ! mass of snow melt (kg/m^2) - - real (kind=dbl_kind), allocatable, dimension(:,:) :: & - l_rsnw , & ! snow grain radius (10^-6 m) - l_smice , & ! tracer for mass of ice in snow (kg/m^3) - l_smliq ! tracer for mass of liquid in snow (kg/m^3) - - real (kind=dbl_kind) :: & - l_fsloss , & ! rate of snow loss to leads (kg/m^2/s) - l_meltsliq , & ! mass of snow melt (kg/m^2) - l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg) - l_H2_16O_ocn, & ! local ocean concentration of H2_16O (kg/kg) - l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg) - - real (kind=dbl_kind) :: & - l_fswthru_vdr , & ! vis dir SW through ice to ocean (W/m^2) - l_fswthru_vdf , & ! vis dif SW through ice to ocean (W/m^2) - l_fswthru_idr , & ! nir dir SW through ice to ocean (W/m^2) - l_fswthru_idf ! nir dif SW through ice to ocean (W/m^2) - - real (kind=dbl_kind), dimension(:), allocatable :: & - l_fswthrun_vdr , & ! vis dir SW through ice to ocean (W/m^2) - l_fswthrun_vdf , & ! vis dif SW through ice to ocean (W/m^2) - l_fswthrun_idr , & ! nir dir SW through ice to ocean (W/m^2) - l_fswthrun_idf ! nir dif SW through ice to ocean (W/m^2) + l_meltsliqn ! mass of snow melt local (kg/m^2) + + real (kind=dbl_kind) :: & + l_fswthrun_vdr, & ! vis dir SW local n ice to ocean (W/m^2) + l_fswthrun_vdf, & ! vis dif SW local n ice to ocean (W/m^2) + l_fswthrun_idr, & ! nir dir SW local n ice to ocean (W/m^2) + l_fswthrun_idf, & ! nir dif SW local n ice to ocean (W/m^2) + l_meltsliq ! mass of snow melt local (kg/m^2) real (kind=dbl_kind) :: & - pond ! water retained in ponds (m) + pond ! water retained in ponds (m) + + logical (kind=log_kind), save :: & + first_call = .true. ! first call flag character(len=*),parameter :: subname='(icepack_step_therm1)' !----------------------------------------------------------------- - ! allocate local optional arguments + ! check optional arguments !----------------------------------------------------------------- - if (present(isosno) ) then - allocate(l_isosno(size(isosno,dim=1),size(isosno,dim=2))) - l_isosno = isosno - else - allocate(l_isosno(1,1)) - l_isosno = c0 - endif - - if (present(isoice) ) then - allocate(l_isoice(size(isoice,dim=1),size(isoice,dim=2))) - l_isoice = isoice - else - allocate(l_isoice(1,1)) - l_isoice = c0 - endif - - if (present(Qa_iso) ) then - allocate(l_Qa_iso(size(Qa_iso))) - l_Qa_iso = Qa_iso - else - allocate(l_Qa_iso(1)) - l_Qa_iso = c0 - endif - - if (present(Qref_iso) ) then - allocate(l_Qref_iso(size(Qref_iso))) - l_Qref_iso = Qref_iso - else - allocate(l_Qref_iso(1)) - l_Qref_iso = c0 - endif - - if (present(fiso_atm) ) then - allocate(l_fiso_atm(size(fiso_atm))) - l_fiso_atm = fiso_atm - else - allocate(l_fiso_atm(1)) - l_fiso_atm = c0 - endif - - if (present(fiso_ocn) ) then - allocate(l_fiso_ocn(size(fiso_ocn))) - l_fiso_ocn = fiso_ocn - else - allocate(l_fiso_ocn(1)) - l_fiso_ocn = c0 - endif - - if (present(fiso_evap) ) then - allocate(l_fiso_evap(size(fiso_evap))) - l_fiso_evap = fiso_evap - else - allocate(l_fiso_evap(1)) - l_fiso_evap = c0 + if (icepack_chkoptargflag(first_call)) then + if (tr_iso) then + if (.not.(present(isosno) .and. present(isoice) .and. & + present(fiso_atm) .and. present(fiso_ocn) .and. & + present(fiso_evap).and. & + present(Qa_iso) .and. present(Qref_iso) .and. & + present(HDO_ocn) .and. present(H2_16O_ocn) .and. & + present(H2_18O_ocn))) then + call icepack_warnings_add(subname//' error in iso arguments, tr_iso=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif + if (snwgrain) then + if (.not.(present(rsnwn) .and. & + present(smicen) .and. present(smliqn))) then + call icepack_warnings_add(subname//' error in snw arguments, snwgrain=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + endif + if ((present(fswthru_vdr) .and. .not.present(fswthrun_vdr)) .or. & + (present(fswthru_vdf) .and. .not.present(fswthrun_vdf)) .or. & + (present(fswthru_idr) .and. .not.present(fswthrun_idr)) .or. & + (present(fswthru_idf) .and. .not.present(fswthrun_idf))) then + call icepack_warnings_add(subname//' error in fswthru [iv]d[rf] arguments') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif endif - l_fsloss = c0 - if (present(fsloss) ) l_fsloss = fsloss - - l_HDO_ocn = c0 - if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn - - l_H2_16O_ocn = c0 - if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn - - l_H2_18O_ocn = c0 - if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn - - l_fswthru_vdr = c0 - if (present(fswthru_vdr)) l_fswthru_vdr = fswthru_vdr - - l_fswthru_vdf = c0 - if (present(fswthru_vdf)) l_fswthru_vdf = fswthru_vdf - - l_fswthru_idr = c0 - if (present(fswthru_idr)) l_fswthru_idr = fswthru_idr - - l_fswthru_idf = c0 - if (present(fswthru_idf)) l_fswthru_idf = fswthru_idf - - allocate(l_fswthrun_vdr(ncat)) - l_fswthrun_vdr = c0 - if (present(fswthrun_vdr)) l_fswthrun_vdr = fswthrun_vdr - - allocate(l_fswthrun_vdf(ncat)) - l_fswthrun_vdf = c0 - if (present(fswthrun_vdf)) l_fswthrun_vdf = fswthrun_vdf - - allocate(l_fswthrun_idr(ncat)) - l_fswthrun_idr = c0 - if (present(fswthrun_idr)) l_fswthrun_idr = fswthrun_idr + !----------------------------------------------------------------- + ! allocate local optional arguments + !----------------------------------------------------------------- - allocate(l_fswthrun_idf(ncat)) - l_fswthrun_idf = c0 - if (present(fswthrun_idf)) l_fswthrun_idf = fswthrun_idf + rsnw (:) = c0 + smice(:) = c0 + smliq(:) = c0 allocate(l_meltsliqn(ncat)) l_meltsliqn = c0 @@ -2536,23 +2463,12 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & l_meltsliq = c0 if (present(meltsliq )) l_meltsliq = meltsliq - allocate(l_rsnw(nslyr,ncat)) - l_rsnw = rsnw_fall - if (present(rsnwn)) l_rsnw = rsnwn - - allocate(l_smice(nslyr,ncat)) - l_smice = c0 - if (present(smicen)) l_smice = smicen - - allocate(l_smliq(nslyr,ncat)) - l_smliq = c0 - if (present(smliqn)) l_smliq = smliqn - !----------------------------------------------------------------- ! Initialize rate of snow loss to leads !----------------------------------------------------------------- - l_fsloss = fsnow * (c1 - aice) + if (present(fsloss)) & + fsloss = fsnow * (c1 - aice) !----------------------------------------------------------------- ! snow redistribution using snwlvlfac: precip factor @@ -2566,8 +2482,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & enddo worka = worka * (snwlvlfac/(c1+snwlvlfac)) / aice endif - l_fsloss = l_fsloss + fsnow* worka - fsnow = fsnow*(c1-worka) + if (present(fsloss)) & + fsloss = fsloss + fsnow* worka + fsnow = fsnow*(c1-worka) endif ! snwredist !----------------------------------------------------------------- @@ -2580,8 +2497,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & rnslyr = c1 / real(nslyr, dbl_kind) do n = 1, ncat do k = 1, nslyr - massicen(k,n) = l_smice(k,n) * vsnon(n) * rnslyr ! kg/m^2 - massliqn(k,n) = l_smliq(k,n) * vsnon(n) * rnslyr + massicen(k,n) = smicen(k,n) * vsnon(n) * rnslyr ! kg/m^2 + massliqn(k,n) = smliqn(k,n) * vsnon(n) * rnslyr enddo enddo endif @@ -2638,7 +2555,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & congeln(n) = c0 snoicen(n) = c0 dsnown (n) = c0 - meltsliqn(n) = c0 + l_meltsliqn(n) = c0 Trefn = c0 Qrefn = c0 @@ -2675,7 +2592,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & lhcoef, shcoef, & Cdn_atm, & Cdn_atm_ratio_n, & - Qa_iso=l_Qa_iso, & + Qa_iso=Qa_iso, & Qref_iso=Qrefn_iso, & uvel=uvel, vvel=vvel, & Uref=Urefn, zlvs=zlvs) @@ -2731,35 +2648,41 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (icepack_warnings_aborted(subname)) return endif - call thermo_vertical(nilyr=nilyr, nslyr=nslyr, & - dt=dt, aicen=aicen (n), & - vicen=vicen (n), vsnon=vsnon (n), & - Tsf=Tsfc (n), zSin=zSin (:,n), & - zqin=zqin (:,n), zqsn=zqsn (:,n), & - apond=apnd (n), hpond=hpnd (n), & - flw=flw, potT=potT, & - Qa=Qa, rhoa=rhoa, & - fsnow=fsnow, fpond=fpond, & - fbot=fbot, Tbot=Tbot, & - Tsnice=Tsnice, sss=sss, & - rsnw=l_rsnw (:,n), & - lhcoef=lhcoef, shcoef=shcoef, & - fswsfc=fswsfcn (n), fswint=fswintn (n), & - Sswabs=Sswabsn(:,n), Iswabs=Iswabsn(:,n), & - fsurfn=fsurfn (n), fcondtopn=fcondtopn(n), & - fcondbotn=fcondbotn(n), & - fsensn=fsensn (n), flatn=flatn (n), & - flwoutn=flwoutn, evapn=evapn, & - evapsn=evapsn, evapin=evapin, & - freshn=freshn, fsaltn=fsaltn, & - fhocnn=fhocnn, frain=frain, & - meltt=melttn (n), melts=meltsn (n), & - meltb=meltbn (n), meltsliq=l_meltsliqn(n),& - smice=l_smice (:,n), massice=massicen(:,n), & - smliq=l_smliq (:,n), massliq=massliqn(:,n), & - congel=congeln (n), snoice=snoicen (n), & - mlt_onset=mlt_onset, frz_onset=frz_onset, & - yday=yday, dsnow=dsnown (n), & + if (snwgrain) then + rsnw (:) = rsnwn (:,n) + smice(:) = smicen(:,n) + smliq(:) = smliqn(:,n) + endif + + call thermo_vertical(nilyr=nilyr, nslyr=nslyr, & + dt=dt, aicen=aicen (n), & + vicen=vicen (n), vsnon=vsnon (n), & + Tsf=Tsfc (n), zSin=zSin (:,n), & + zqin=zqin (:,n), zqsn=zqsn (:,n), & + apond=apnd (n), hpond=hpnd (n), & + flw=flw, potT=potT, & + Qa=Qa, rhoa=rhoa, & + fsnow=fsnow, fpond=fpond, & + fbot=fbot, Tbot=Tbot, & + Tsnice=Tsnice, sss=sss, & + rsnw=rsnw, & + lhcoef=lhcoef, shcoef=shcoef, & + fswsfc=fswsfcn (n), fswint=fswintn (n), & + Sswabs=Sswabsn(:,n), Iswabs=Iswabsn (:,n), & + fsurfn=fsurfn (n), fcondtopn=fcondtopn (n), & + fcondbotn=fcondbotn (n), & + fsensn=fsensn (n), flatn=flatn (n), & + flwoutn=flwoutn, evapn=evapn, & + evapsn=evapsn, evapin=evapin, & + freshn=freshn, fsaltn=fsaltn, & + fhocnn=fhocnn, frain=frain, & + meltt=melttn (n), melts=meltsn (n), & + meltb=meltbn (n), meltsliq=l_meltsliqn(n), & + smice=smice, massice=massicen (:,n), & + smliq=smliq, massliq=massliqn (:,n), & + congel=congeln (n), snoice=snoicen (n), & + mlt_onset=mlt_onset, frz_onset=frz_onset, & + yday=yday, dsnow=dsnown (n), & prescribed_ice=prescribed_ice) if (icepack_warnings_aborted(subname)) then @@ -2768,6 +2691,12 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & return endif + if (snwgrain) then + rsnwn (:,n) = rsnw (:) + smicen(:,n) = smice(:) + smliqn(:,n) = smliq(:) + endif + !----------------------------------------------------------------- ! Total absorbed shortwave radiation !----------------------------------------------------------------- @@ -2795,22 +2724,25 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & if (tr_iso) then call update_isotope (dt = dt, & - nilyr = nilyr, nslyr = nslyr, & - meltt = melttn(n),melts = meltsn(n), & - meltb = meltbn(n),congel=congeln(n), & - snoice=snoicen(n),evap=evapn, & - fsnow=fsnow, Tsfc=Tsfc(n), & + nilyr = nilyr, nslyr = nslyr, & + meltt = melttn(n),melts = meltsn(n), & + meltb = meltbn(n),congel=congeln(n), & + snoice=snoicen(n),evap=evapn, & + fsnow=fsnow, Tsfc=Tsfc(n), & Qref_iso=Qrefn_iso(:), & - isosno=l_isosno(:,n),isoice=l_isoice(:,n), & - aice_old=aicen_init(n),vice_old=vicen_init(n), & + isosno=isosno(:,n),isoice=isoice(:,n), & + aice_old=aicen_init(n), & + vice_old=vicen_init(n), & vsno_old=vsnon_init(n), & - vicen=vicen(n),vsnon=vsnon(n), & - aicen=aicen(n), & - fiso_atm=l_fiso_atm(:), & - fiso_evapn=fiso_evapn(:), & - fiso_ocnn=fiso_ocnn(:), & - HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn, & - H2_18O_ocn=l_H2_18O_ocn) + vicen=vicen(n), & + vsnon=vsnon(n), & + aicen=aicen(n), & + fiso_atm=fiso_atm(:), & + fiso_evapn=fiso_evapn(:), & + fiso_ocnn=fiso_ocnn(:), & + HDO_ocn=HDO_ocn, & + H2_16O_ocn=H2_16O_ocn, & + H2_18O_ocn=H2_18O_ocn) if (icepack_warnings_aborted(subname)) return endif @@ -2897,7 +2829,17 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & ! Increment area-weighted fluxes. !----------------------------------------------------------------- - if (aicen_init(n) > puny) & + if (aicen_init(n) > puny) then + + l_fswthrun_vdr = c0 + l_fswthrun_vdf = c0 + l_fswthrun_idr = c0 + l_fswthrun_idf = c0 + if (present(fswthrun_vdr)) l_fswthrun_vdr = fswthrun_vdr(n) + if (present(fswthrun_vdf)) l_fswthrun_vdf = fswthrun_vdf(n) + if (present(fswthrun_idr)) l_fswthrun_idr = fswthrun_idr(n) + if (present(fswthrun_idf)) l_fswthrun_idf = fswthrun_idf(n) + call merge_fluxes (aicen=aicen_init(n), & flw=flw, & strairxn=strairxn, strairyn=strairyn,& @@ -2912,10 +2854,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & freshn=freshn, fsaltn=fsaltn, & fhocnn=fhocnn, & fswthrun=fswthrun(n), & - fswthrun_vdr=l_fswthrun_vdr(n), & - fswthrun_vdf=l_fswthrun_vdf(n), & - fswthrun_idr=l_fswthrun_idr(n), & - fswthrun_idf=l_fswthrun_idf(n), & + fswthrun_vdr=l_fswthrun_vdr, & + fswthrun_vdf=l_fswthrun_vdf, & + fswthrun_idr=l_fswthrun_idr, & + fswthrun_idf=l_fswthrun_idf, & strairxT=strairxT, strairyT=strairyT,& Cdn_atm_ratio=Cdn_atm_ratio, & fsurf=fsurf, fcondtop=fcondtop,& @@ -2928,10 +2870,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fresh=fresh, fsalt=fsalt, & fhocn=fhocn, & fswthru=fswthru, & - fswthru_vdr=l_fswthru_vdr, & - fswthru_vdf=l_fswthru_vdf, & - fswthru_idr=l_fswthru_idr, & - fswthru_idf=l_fswthru_idf, & + fswthru_vdr=fswthru_vdr, & + fswthru_vdf=fswthru_vdf, & + fswthru_idr=fswthru_idr, & + fswthru_idf=fswthru_idf, & melttn=melttn (n), meltsn=meltsn(n), & meltbn=meltbn (n), congeln=congeln(n),& meltt=meltt, melts=melts, & @@ -2941,14 +2883,16 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & meltsliq=l_meltsliq, & meltsliqn=l_meltsliqn(n), & Uref=Uref, Urefn=Urefn, & - Qref_iso=l_Qref_iso, & + Qref_iso=Qref_iso, & Qrefn_iso=Qrefn_iso, & - fiso_ocn=l_fiso_ocn, & + fiso_ocn=fiso_ocn, & fiso_ocnn=fiso_ocnn, & - fiso_evap=l_fiso_evap, & + fiso_evap=fiso_evap, & fiso_evapn=fiso_evapn) - if (icepack_warnings_aborted(subname)) return + if (icepack_warnings_aborted(subname)) return + + endif enddo ! ncat @@ -2960,59 +2904,26 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & do n = 1, ncat if (vsnon(n) > puny) then do k = 1, nslyr - l_smice(k,n) = massicen(k,n) / (vsnon(n) * rnslyr) - l_smliq(k,n) = massliqn(k,n) / (vsnon(n) * rnslyr) - worka = l_smice(k,n) + l_smliq(k,n) + smicen(k,n) = massicen(k,n) / (vsnon(n) * rnslyr) + smliqn(k,n) = massliqn(k,n) / (vsnon(n) * rnslyr) + worka = smicen(k,n) + smliqn(k,n) if (worka > puny) then - l_smice(k,n) = rhos * l_smice(k,n) / worka - l_smliq(k,n) = rhos * l_smliq(k,n) / worka + smicen(k,n) = rhos * smicen(k,n) / worka + smliqn(k,n) = rhos * smliqn(k,n) / worka endif enddo else ! reset to default values do k = 1, nslyr - l_smice(k,n) = rhos - l_smliq(k,n) = c0 + smicen(k,n) = rhos + smliqn(k,n) = c0 enddo endif enddo endif - if (present(isosno )) isosno = l_isosno - if (present(isoice )) isoice = l_isoice - if (present(Qa_iso )) Qa_iso = l_Qa_iso - if (present(Qref_iso )) Qref_iso = l_Qref_iso - if (present(fiso_atm )) fiso_atm = l_fiso_atm - if (present(fiso_ocn )) fiso_ocn = l_fiso_ocn - if (present(fiso_evap )) fiso_evap = l_fiso_evap - if (present(fswthrun_vdr)) fswthrun_vdr = l_fswthrun_vdr - if (present(fswthrun_vdf)) fswthrun_vdf = l_fswthrun_vdf - if (present(fswthrun_idr)) fswthrun_idr = l_fswthrun_idr - if (present(fswthrun_idf)) fswthrun_idf = l_fswthrun_idf - if (present(fswthru_vdr )) fswthru_vdr = l_fswthru_vdr - if (present(fswthru_vdf )) fswthru_vdf = l_fswthru_vdf - if (present(fswthru_idr )) fswthru_idr = l_fswthru_idr - if (present(fswthru_idf )) fswthru_idf = l_fswthru_idf - if (present(fsloss )) fsloss = l_fsloss if (present(meltsliqn )) meltsliqn = l_meltsliqn if (present(meltsliq )) meltsliq = l_meltsliq - if (present(rsnwn )) rsnwn = l_rsnw - if (present(smicen )) smicen = l_smice - if (present(smliqn )) smliqn = l_smliq - deallocate(l_isosno) - deallocate(l_isoice) - deallocate(l_Qa_iso) - deallocate(l_Qref_iso) - deallocate(l_fiso_atm) - deallocate(l_fiso_ocn) - deallocate(l_fiso_evap) - deallocate(l_fswthrun_vdr) - deallocate(l_fswthrun_vdf) - deallocate(l_fswthrun_idr) - deallocate(l_fswthrun_idf) deallocate(l_meltsliqn) - deallocate(l_rsnw) - deallocate(l_smice) - deallocate(l_smliq) !----------------------------------------------------------------- ! Calculate ponds from the topographic scheme @@ -3033,6 +2944,8 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & endif !call ice_timer_stop(timer_ponds) + first_call = .false. + end subroutine icepack_step_therm1 !=======================================================================