! source file: /Users/oschlies/UVIC/master/source/embm/embmio.F subroutine embmout (is, ie, js, je) !======================================================================= ! output routine for energy-moisture balance model ! input: ! is, ie, js, je = starting and ending indicies for i and j ! based on code by: A. Fanning and M. Eby !======================================================================= implicit none include "param.h" include "calendar.h" include "csbc.h" include "atm.h" include "solve.h" include "coord.h" include "grdvar.h" include "ice.h" include "evp.h" include "mtlm.h" include "cembm.h" include "iounit.h" include "scalar.h" include "switch.h" include "tmngr.h" include "riv.h" include "cregin.h" real tmp_at(imt,jmt,nat) real sm(imt,jmt), st(imt,jmt), hs(imt,jmt), ro(imt,jmt) real rntatsl, rntatil character(120) :: fname integer i, ie, is, id_xt, id_yt, iou, j, je, js, n integer ndx, ntrec, it(10), ib(10), ic(10) real avgper, ca, time, wt, wtp1 real p_alb(is:ie,js:je), a_alb(is:ie,js:je), s_alb(is:ie,js:je) real sat(is:ie,js:je) real c100, c400, C2K, tmp c100 = 100. c400 = 400. C2K = 273.15 !----------------------------------------------------------------------- ! write atmospheric diagnostics !----------------------------------------------------------------------- if (tsits .and. ntatia .ne. 0) then call ta_embm_tsi (2) if (iotsi .eq. stdout .or. iotsi .lt. 0) then write (*,'(1x, a3, i7, 1x, a32, 3(a,1pe13.6))') & 'ts=',itt, stamp, ' iterations =', tai_maxit &, ' TAbar=', tai_sat, ' QAbar=', tai_shum endif if (iotsi .ne. stdout) then avgper = dtatm*ntatia/daylen time = relyr + year0 rntatil = 0. if (ntatil .gt. 0) rntatil = 1./float(ntatil) tai_hsno = tai_hsno + tai_LYING_SNOW*1000.*rntatil/rhosno call def_tsi call def_tsi_embm (fname) call embm_tsi_out (fname, avgper, time, stamp, tai_sat &, tai_shum, tai_precip, tai_evap, tai_ohice &, tai_oaice, tai_hsno, tai_lhice, tai_laice &, tai_co2ccn, tai_dc14ccn, tai_cfc11ccn &, tai_cfc12ccn, tai_maxit, tai_sst, tai_sss &, tai_nsat, tai_ssat, tai_nshum, tai_sshum &, tai_nprecip, tai_sprecip, tai_nevap &, tai_sevap, tai_nohice, tai_sohice &, tai_noaice, tai_soaice, tai_nhsno &, tai_shsno, tai_nlhice, tai_slhice &, tai_nlaice, tai_slaice, ntrec) endif call ta_embm_tsi (0) endif if (timavgts .and. ntatsa .ne. 0) then !----------------------------------------------------------------------- ! write atmospheric time averaged data !----------------------------------------------------------------------- ! calculate average values call ta_embm_snap (is, ie, js, je, 2) ! write time averaged data rntatsl = 0. if (ntatsl .gt. 0) rntatsl = 1./float(ntatsl) call unloadland (POINTS, TA_M, imt, jmt, land_map, sm) call unloadland (POINTS, TA_TSTAR_GB, imt, jmt, land_map, st) call unloadland (POINTS, TA_LYING_SNOW, imt, jmt, land_map, hs) call unloadland (POINTS, TA_SURF_ROFF, imt, jmt, land_map, ro) do j=js,je do i=is,ie sat(i,j) = ta_at(i,j,isat) - elev(i,j)*rlapse & - hicel(i,j,2)*rlapse if (ta_solins(i,j) .gt. 0) then ! calculate incoming swr reaching the surface tmp = ta_solins(i,j) - ta_arswr(i,j) - ta_absin(i,j) s_alb(i,j) = 1. - ta_dnswr(i,j)/tmp a_alb(i,j) = ta_arswr(i,j)/ta_solins(i,j) ! calculate outgoing swr leaving the atmosphere tmp = tmp - ta_dnswr(i,j) - ta_absout(i,j) p_alb(i,j) = (ta_arswr(i,j) + tmp)/ta_solins(i,j) else s_alb(i,j) = -1. a_alb(i,j) = -1. p_alb(i,j) = -1. endif if (land_map(i,j) .eq. 0) then sm(i,j) = ta_soilm(i,j) st(i,j) = ta_surf(i,j) else ! convert from kg m-2 to cm (converted back later) sm(i,j) = sm(i,j)*rntatsl*0.1 ! convert from K to C (converted back later) st(i,j) = st(i,j)*rntatsl - 273.15 ! convert from kg/m2/s to g/cm2/s (converted back later) ta_runoff(i,j) = ro(i,j)*rntatsl*0.1 ! convert from kg/m2 to cm (converted back later) ta_hsno(i,j,1) = hs(i,j)*rntatsl*0.1/rhosno endif enddo enddo avgper = dtatm*ntatsa/daylen time = relyr + year0 call def_tavg call def_tavg_embm (fname) call embm_snap_out (fname, is, ie, js, je, imt, jmt, nat, ncat &, xt, yt, xu, yu, dxt, dyt, dxu, dyu, avgper, time, stamp, mapat &, ta_at, sat, ta_precip, ta_evap, ta_disch, ta_outlwr, ta_uplwr &, ta_upsens, ta_dnswr, ta_upltnt, ta_solins, p_alb, a_alb, s_alb &, elev, ta_psno, ta_ws, ta_runoff &, ta_wx_q, ta_wy_q &, ta_wx_t, ta_wy_t &, sm, st &, ta_tice, ta_hice, ta_aice, ta_hsno &, ta_uice, ta_vice, ta_xint, ta_yint &, tlat, tlon, ulat, ulon &, ta_aicel(is,js), ta_hicel(is,js) &, tmsk, mskhr, nriv, ntrec) write (*,'(a,i5,a,a,a,a)') '=> Atm time means #' &, ntrec, ' written to ',trim(fname),' on ', stamp ! zero time average accumulators call ta_embm_snap (is, ie, js, je, 0) endif if (snapts) then !----------------------------------------------------------------------- ! write atmospheric snapshot !----------------------------------------------------------------------- ! write snapshot tmp_at(is:ie,js:je,1:nat) = at(is:ie,js:je,2,1:nat) call unloadland (POINTS, M, imt, jmt, land_map, sm) call unloadland (POINTS, TSTAR_GB, imt, jmt, land_map, st) call unloadland (POINTS, LYING_SNOW, imt, jmt, land_map, hs) call unloadland (POINTS, SURF_ROFF, imt, jmt, land_map, ro) do j=js,je do i=is,ie sat(i,j) = at(i,j,2,isat) - elev(i,j)*rlapse & - hicel(i,j,2)*rlapse if (solins(i,j) .gt. 0) then ! calculate incoming swr reaching the surface tmp = solins(i,j)*sbc(i,j,iaca)*pass s_alb(i,j) = 1. - dnswr(i,j)/tmp a_alb(i,j) = 1. - sbc(i,j,iaca) ! calculate outgoing swr leaving the atmosphere tmp = tmp - dnswr(i,j) p_alb(i,j) = a_alb(i,j) + tmp/solins(i,j) else s_alb(i,j) = -1. a_alb(i,j) = -1. p_alb(i,j) = -1. endif if (land_map(i,j) .eq. 0) then sm(i,j) = soilm(i,j,2) st(i,j) = surf(i,j) hs(i,j) = hsno(i,j,2) else ! convert from kg m-2 to cm (converted back later) sm(i,j) = sm(i,j)*0.1 ! convert from K to C (converted back later) st(i,j) = st(i,j) - 273.15 ! convert from kg/m2 to cm (converted back later) hs(i,j) = hs(i,j)*0.1/rhosno endif enddo enddo avgper = 0. time = relyr + year0 call def_snap call def_snap_embm (fname) call embm_snap_out (fname, is, ie, js, je, imt, jmt, nat, ncat &, xt, yt, xu, yu, dxt, dyt, dxu, dyu, avgper, time, stamp, mapat &, tmp_at, sat, precip, evap, disch, outlwr, uplwr, upsens, dnswr &, upltnt, solins, p_alb, a_alb, s_alb, elev, psno &, sbc(is,js,iws), sbc(is,js,iro) &, sbc(is,js,iwxq), sbc(is,js,iwyq) &, sbc(is,js,iwxt), sbc(is,js,iwyt) &, sm, st &, tice(is,js), hice(is,js,2), aice(is,js,2), hs(is,js) &, uice, vice, xint, yint &, tlat, tlon, ulat, ulon &, aicel(is,js,2), hicel(is,js,2) &, tmsk, mskhr, nriv, ntrec) write (*,'(a,i5,a,a,a,a)') '=> Atm snapshot #' &, ntrec, ' written to ',trim(fname),' on ', stamp endif !----------------------------------------------------------------------- ! write atmospheric restart !----------------------------------------------------------------------- if (eoyear) call eoyr (1, imt, 1, jmt) if (restrt) then if (restts) then call def_rest (0) call def_rest_embm (0, fname) call embm_rest_out (fname, is, ie, js, je) endif if (eorun) then call def_rest (1) call def_rest_embm (1, fname) call embm_rest_out (fname, is, ie, js, je) endif endif return end subroutine ta_embm_snap (is, ie, js, je, m) !======================================================================= ! atmospheric data time averaging ! input: ! is, ie, js, je = starting and ending indicies for i and j ! m = flag (0 = zero, 1 = accumulate, 2 = write) ! based on code by: M. Eby !======================================================================= implicit none include "param.h" include "cembm.h" include "atm.h" include "riv.h" include "ice.h" include "csbc.h" integer i, ie, is, j, je, js, k, m, n, ndx real rntatsa !----------------------------------------------------------------------- ! time averaged data !----------------------------------------------------------------------- if (m .eq. 0.) then ! zero ntatsa = 0 ta_at(:,:,:) = 0. ta_solins(:,:) = 0. ta_arswr(:,:) = 0. ta_dnswr(:,:) = 0. ta_absin(:,:) = 0. ta_absout(:,:) = 0. ta_uplwr(:,:) = 0. ta_upsens(:,:) = 0. ta_upltnt(:,:) = 0. ta_outlwr(:,:) = 0. ta_precip(:,:) = 0. ta_psno(:,:) = 0. ta_evap(:,:) = 0. ta_disch(:,:) = 0. ta_ws(:,:) = 0. ta_runoff(:,:) = 0. ta_wx_q(:,:) = 0. ta_wy_q(:,:) = 0. ta_wx_t(:,:) = 0. ta_wy_t(:,:) = 0. ta_soilm(:,:) = 0. ta_surf(:,:) = 0. ta_tice(:,:,:) = 0. ta_hice(:,:,:) = 0. ta_aice(:,:,:) = 0. ta_hsno(:,:,:) = 0. ta_uice(:,:) = 0. ta_vice(:,:) = 0. ta_pice(:,:) = 0. ta_xint(:,:) = 0. ta_yint(:,:) = 0. ta_aicel(:,:) = 0. ta_hicel(:,:) = 0. ta_psum(:) = 0. elseif (m .eq. 1) then ! accumulate ntatsa = ntatsa + 1 do n=1,nat ta_at(:,:,n) = ta_at(:,:,n) + at(:,:,2,n) enddo ta_solins(:,:) = ta_solins(:,:) + solins(:,:) ta_arswr(:,:) = ta_arswr(:,:) + solins(:,:) & *(1. - sbc(:,:,iaca)) ta_dnswr(:,:) = ta_dnswr(:,:) + dnswr(:,:) ta_absin(:,:) = ta_absin(:,:) + solins(:,:)*sbc(:,:,iaca) & *scatter ta_absout(:,:) = ta_absout(:,:) + (solins(:,:)*sbc(:,:,iaca) & *pass - dnswr(:,:))*scatter ta_uplwr(:,:) = ta_uplwr(:,:) + uplwr(:,:) ta_upsens(:,:) = ta_upsens(:,:) + upsens(:,:) ta_upltnt(:,:) = ta_upltnt(:,:) + upltnt(:,:) ta_outlwr(:,:) = ta_outlwr(:,:) + outlwr(:,:) ta_precip(:,:) = ta_precip(:,:) + precip(:,:) ta_psno(:,:) = ta_psno(:,:) + psno(:,:) ta_evap(:,:) = ta_evap(:,:) + evap(:,:) ta_disch(:,:) = ta_disch(:,:) + disch(:,:) ta_ws(:,:) = ta_ws(:,:) + sbc(:,:,iws) ta_runoff(:,:) = ta_runoff(:,:) + runoff(:,:) ta_wx_q(:,:) = ta_wx_q(:,:) + sbc(:,:,iwxq) ta_wy_q(:,:) = ta_wy_q(:,:) + sbc(:,:,iwyq) ta_wx_t(:,:) = ta_wx_t(:,:) + sbc(:,:,iwxt) ta_wy_t(:,:) = ta_wy_t(:,:) + sbc(:,:,iwyt) ta_soilm(:,:) = ta_soilm(:,:) + soilm(:,:,2) ta_surf(:,:) = ta_surf(:,:) + surf(:,:) ta_tice(:,:,1) = ta_tice(:,:,1) + tice(:,:) ta_hice(:,:,1) = ta_hice(:,:,1) + hice(:,:,2) ta_aice(:,:,1) = ta_aice(:,:,1) + aice(:,:,2) ta_hsno(:,:,1) = ta_hsno(:,:,1) + hsno(:,:,2) ta_uice(:,:) = ta_uice(:,:) + uice(:,:) ta_vice(:,:) = ta_vice(:,:) + vice(:,:) ta_pice(:,:) = ta_pice(:,:) + pice(:,:) ta_xint(:,:) = ta_xint(:,:) + xint(:,:) ta_yint(:,:) = ta_yint(:,:) + yint(:,:) ta_aicel(:,:) = ta_aicel(:,:) + aicel(:,:,2) ta_hicel(:,:) = ta_hicel(:,:) + hicel(:,:,2) ta_psum(:) = ta_psum(:) + psum(:) elseif (m .eq. 2 .and. ntatsa .ne. 0) then ! average rntatsa = 1./float(ntatsa) ta_at(:,:,:) = ta_at(:,:,:)*rntatsa ta_solins(:,:) = ta_solins(:,:)*rntatsa ta_arswr(:,:) = ta_arswr(:,:)*rntatsa ta_dnswr(:,:) = ta_dnswr(:,:)*rntatsa ta_absin(:,:) = ta_absin(:,:)*rntatsa ta_absout(:,:) = ta_absout(:,:)*rntatsa ta_uplwr(:,:) = ta_uplwr(:,:)*rntatsa ta_upsens(:,:) = ta_upsens(:,:)*rntatsa ta_upltnt(:,:) = ta_upltnt(:,:)*rntatsa ta_outlwr(:,:) = ta_outlwr(:,:)*rntatsa ta_precip(:,:) = ta_precip(:,:)*rntatsa ta_psno(:,:) = ta_psno(:,:)*rntatsa ta_evap(:,:) = ta_evap(:,:)*rntatsa ta_disch(:,:) = ta_disch(:,:)*rntatsa ta_ws(:,:) = ta_ws(:,:)*rntatsa ta_runoff(:,:) = ta_runoff(:,:)*rntatsa ta_wx_q(:,:) = ta_wx_q(:,:)*rntatsa ta_wy_q(:,:) = ta_wy_q(:,:)*rntatsa ta_wx_t(:,:) = ta_wx_t(:,:)*rntatsa ta_wy_t(:,:) = ta_wy_t(:,:)*rntatsa ta_soilm(:,:) = ta_soilm(:,:)*rntatsa ta_surf(:,:) = ta_surf(:,:)*rntatsa ta_hsno(:,:,:) = ta_hsno(:,:,:)*rntatsa ta_tice(:,:,:) = ta_tice(:,:,:)*rntatsa ta_hice(:,:,:) = ta_hice(:,:,:)*rntatsa ta_aice(:,:,:) = ta_aice(:,:,:)*rntatsa ta_uice(:,:) = ta_uice(:,:)*rntatsa ta_vice(:,:) = ta_vice(:,:)*rntatsa ta_pice(:,:) = ta_pice(:,:)*rntatsa ta_xint(:,:) = ta_xint(:,:)*rntatsa ta_yint(:,:) = ta_yint(:,:)*rntatsa ta_aicel(:,:) = ta_aicel(:,:)*rntatsa ta_hicel(:,:) = ta_hicel(:,:)*rntatsa ta_psum(:) = ta_psum(:)*rntatsa endif return end subroutine ta_embm_tsi (m) !======================================================================= ! atmospheric data time integral averaging ! input: ! m = flag (0 = zero, 1 = accumulate, 2 = write) ! based on code by: M. Eby !======================================================================= implicit none include "param.h" include "csbc.h" include "grdvar.h" include "cembm.h" include "atm.h" include "ice.h" include "solve.h" integer maxit, m, n real rntatia, sum, tmp, sat(imt,jmt), dmsk(imt,jmt) !----------------------------------------------------------------------- ! time averaged integrated data !----------------------------------------------------------------------- if (m .eq. 0.) then ! zero ntatia = 0 tai_maxit = 0. tai_sat = 0. tai_shum = 0. tai_co2ccn = 0. tai_dc14ccn = 0. tai_cfc11ccn = 0. tai_cfc12ccn = 0. tai_precip = 0. tai_evap = 0. tai_ohice = 0. tai_oaice = 0. tai_hsno = 0. tai_lhice = 0. tai_laice = 0. tai_sst = 0. tai_sss = 0. tai_nsat = 0. tai_ssat = 0. tai_nshum = 0. tai_sshum = 0. tai_nprecip = 0. tai_sprecip = 0. tai_nevap = 0. tai_sevap = 0. tai_nohice = 0. tai_sohice = 0. tai_noaice = 0. tai_soaice = 0. tai_nhsno = 0. tai_shsno = 0. tai_nlhice = 0. tai_slhice = 0. tai_nlaice = 0. tai_slaice = 0. elseif (m .eq. 1) then ! accumulate ntatia = ntatia + 1 sat(:,:) = at(:,:,2,isat) - elev(:,:)*rlapse & - hicel(:,:,2)*rlapse maxit = 0 ! set data mask dmsk(:,:) = 1. do n=1,nat if (itout(n) .gt. maxit) maxit = itout(n) enddo tai_maxit = tai_maxit + maxit call areaavg (sat, dmsk, tmp) tai_sat = tai_sat + tmp call areaavg (at(1,1,2,ishum), dmsk, tmp) tai_shum = tai_shum + tmp tai_co2ccn = tai_co2ccn + co2ccn tai_dc14ccn = tai_dc14ccn + dc14ccn call areaavg (precip, dmsk, tmp) tai_precip = tai_precip + tmp call areaavg (evap, dmsk, tmp) tai_evap = tai_evap + tmp call areatot (hsno(1,1,2), dmsk, tmp) tai_hsno = tai_hsno + tmp ! mask for ocean where (tmsk(:,:) .lt. 0.5) dmsk(:,:) = 0. call areatot (hice(1,1,2), dmsk, tmp) tai_ohice = tai_ohice + tmp call areatot (aice(1,1,2), dmsk, tmp) tai_oaice = tai_oaice + tmp ! mask for land dmsk(:,:) = 1. where (tmsk(:,:) .ge. 0.5) dmsk(:,:) = 0. call areatot (hicel(1,1,2), dmsk, tmp) tai_lhice = tai_lhice + tmp call areatot (aicel(1,1,2), dmsk, tmp) tai_laice = tai_laice + tmp ! mask for ocean dmsk(:,:) = 1. where (tmsk(:,:) .lt. 0.5) dmsk(:,:) = 0. call areaavg (sbc(1,1,isst), dmsk, tmp) tai_sst = tai_sst + tmp call areaavg (sbc(1,1,isss), dmsk, tmp) tai_sss = tai_sss + tmp ! northern hemisphere mask dmsk(:,:) = 1. where (tlat(:,:) .lt. 0.) dmsk(:,:) = 0. call areaavg (sat, dmsk, tmp) tai_nsat = tai_nsat + tmp call areaavg (at(1,1,2,ishum), dmsk, tmp) tai_nshum = tai_nshum + tmp call areaavg (precip, dmsk, tmp) tai_nprecip = tai_nprecip + tmp call areaavg (evap, dmsk, tmp) tai_nevap = tai_nevap + tmp call areatot (hsno(1,1,2), dmsk, tmp) tai_nhsno = tai_nhsno + tmp ! mask for northern hemisphere ocean where (tmsk(:,:) .lt. 0.5) dmsk(:,:) = 0. call areatot (hice(1,1,2), dmsk, tmp) tai_nohice = tai_nohice + tmp call areatot (aice(1,1,2), dmsk, tmp) tai_noaice = tai_noaice + tmp ! mask for northern hemisphere land dmsk(:,:) = 1. where (tlat(:,:) .lt. 0. .or. tmsk(:,:) .ge. 0.5) dmsk(:,:) = 0. call areatot (hicel(1,1,2), dmsk, tmp) tai_nlhice = tai_nlhice + tmp call areatot (aicel(1,1,2), dmsk, tmp) tai_nlaice = tai_nlaice + tmp ! southern hemisphere mask dmsk(:,:) = 1. where (tlat(:,:) .gt. 0.) dmsk(:,:) = 0. call areaavg (sat, dmsk, tmp) tai_ssat = tai_ssat + tmp call areaavg (at(1,1,2,ishum), dmsk, tmp) tai_sshum = tai_sshum + tmp call areaavg (precip, dmsk, tmp) tai_sprecip = tai_sprecip + tmp call areaavg (evap, dmsk, tmp) tai_sevap = tai_sevap + tmp call areatot (hsno(1,1,2), dmsk, tmp) tai_shsno = tai_shsno + tmp ! mask for southern hemisphere ocean where (tmsk(:,:) .lt. 0.5) dmsk(:,:) = 0. call areatot (hice(1,1,2), dmsk, tmp) tai_sohice = tai_sohice + tmp call areatot (aice(1,1,2), dmsk, tmp) tai_soaice = tai_soaice + tmp ! mask for southern hemisphere land dmsk(:,:) = 1. where (tlat(:,:) .gt. 0. .or. tmsk(:,:) .ge. 0.5) dmsk(:,:) = 0. call areatot (hicel(1,1,2), dmsk, tmp) tai_slhice = tai_slhice + tmp call areatot (aicel(1,1,2), dmsk, tmp) tai_slaice = tai_slaice + tmp elseif (m .eq. 2 .and. ntatia .ne. 0) then ! average rntatia = 0. if (ntatia .gt. 0.) rntatia = 1./float(ntatia) tai_maxit = tai_maxit*rntatia tai_sat = tai_sat*rntatia tai_shum = tai_shum*rntatia tai_co2ccn = tai_co2ccn*rntatia tai_dc14ccn = tai_dc14ccn*rntatia tai_cfc11ccn = tai_cfc11ccn*rntatia tai_cfc12ccn = tai_cfc12ccn*rntatia tai_precip = tai_precip*rntatia tai_evap = tai_evap*rntatia tai_ohice = tai_ohice*rntatia tai_oaice = tai_oaice*rntatia tai_hsno = tai_hsno*rntatia tai_lhice = tai_lhice*rntatia tai_laice = tai_laice*rntatia tai_sst = tai_sst*rntatia tai_sss = tai_sss*rntatia tai_nsat = tai_nsat*rntatia tai_ssat = tai_ssat*rntatia tai_nshum = tai_nshum*rntatia tai_sshum = tai_sshum*rntatia tai_nprecip = tai_nprecip*rntatia tai_sprecip = tai_sprecip*rntatia tai_nevap = tai_nevap*rntatia tai_sevap = tai_sevap*rntatia tai_nohice = tai_nohice*rntatia tai_sohice = tai_sohice*rntatia tai_noaice = tai_noaice*rntatia tai_soaice = tai_soaice*rntatia tai_nhsno = tai_nhsno*rntatia tai_shsno = tai_shsno*rntatia tai_nlhice = tai_nlhice*rntatia tai_slhice = tai_slhice*rntatia tai_nlaice = tai_nlaice*rntatia tai_slaice = tai_slaice*rntatia endif return end