MODULE physiq_mod IMPLICIT NONE CONTAINS SUBROUTINE physiq( $ ngrid,nlayer,nq $ ,firstcall,lastcall $ ,pday,ptime,ptimestep $ ,pplev,pplay,pphi $ ,pu,pv,pt,pq $ ,flxw $ ,pdu,pdv,pdt,pdq,pdpsrf) use watercloud_mod, only: watercloud, zdqcloud, zdqscloud use calchim_mod, only: calchim, ichemistry, zdqchim, zdqschim use watersat_mod, only: watersat use co2condens_mod, only: co2condens use co2cloud_mod, only: co2cloud use callradite_mod, only: callradite use callsedim_mod, only: callsedim use rocketduststorm_mod, only: rocketduststorm, dustliftday use calcstormfract_mod, only: calcstormfract use topmons_mod, only: topmons,topmons_setup use nltecool_mod, only: nltecool use nlte_tcool_mod, only: nlte_tcool use blendrad_mod, only: blendrad use nlthermeq_mod, only: nlthermeq use thermosphere_mod, only: thermosphere use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice, & igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice, & igcm_hdo_vap, igcm_hdo_ice, & igcm_ccn_mass, igcm_ccn_number, & igcm_ccnco2_mass, igcm_ccnco2_number, & igcm_ccnco2_h2o_mass_ice, & igcm_ccnco2_h2o_mass_ccn, & igcm_ccnco2_h2o_number, & igcm_ccnco2_meteor_mass, & igcm_ccnco2_meteor_number, & rho_ice_co2,nuiceco2_sed,nuiceco2_ref, & igcm_dust_mass, igcm_dust_number, igcm_h2o2, & nuice_ref, rho_ice, rho_dust, ref_r0, & igcm_he, igcm_stormdust_mass, & igcm_stormdust_number, igcm_topdust_mass, & igcm_topdust_number, & qperemin use comsoil_h, only: inertiedat, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start & tsoil, nsoilmx,!number of subsurface layers & mlayer,layer, ! soil mid layer depths & nqsoil,qsoil ! adsorption use geometry_mod, only: longitude, latitude, cell_area, & cell_area_for_lonlat_outputs,longitude_deg use comgeomfi_h, only: sinlon, coslon, sinlat, coslat use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, & zthe, z0, albedo_h2o_cap,albedo_h2o_frost, & frost_albedo_threshold,frost_metam_threshold, & tsurf, emis, & capcal, fluxgrd, qsurf, & hmons,summit,base,watercap,watercaptag, & perenial_co2ice use comsaison_h, only: dist_sol, declin, zls, & mu0, fract, local_time use solarlong_mod, only: solarlong use nirdata_mod, only: NIR_leedat use nirco2abs_mod, only: nirco2abs use slope_mod, only: theta_sl, psi_sl use conc_mod, only: rnew, cpnew, mmean use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec use dimradmars_mod, only: aerosol, totcloudfrac, & dtrad, fluxrad_sky, fluxrad, albedo, & naerkind, iaer_dust_doubleq, & iaer_stormdust_doubleq, iaer_h2o_ice, & flux_1AU use dust_param_mod, only: doubleq, lifting, callddevil, & tauscaling, odpref, dustbin, & dustscaling_mode, dust_rad_adjust, & freedust, reff_driven_IRtoVIS_scenario use turb_mod, only: q2, wstar, ustar, sensibFlux, & zmax_th, hfmax_th, turb_resolved use planete_h, only: aphelie, periheli, year_day, peri_day, & obliquit use planete_h, only: iniorbit USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad USE calldrag_noro_mod, ONLY: calldrag_noro USE vdifc_mod, ONLY: vdifc use param_v4_h, only: nreact,n_avog, & fill_data_thermos, allocate_param_thermos use iono_h, only: allocate_param_iono use compute_dtau_mod, only: compute_dtau use nonoro_gwd_ran_mod, only: nonoro_gwd_ran use check_fields_mod, only: check_physics_fields #ifdef MESOSCALE use comsoil_h, only: mlayer,layer use surfdat_h, only: z0_default use comm_wrf #else USE planetwide_mod, ONLY: planetwide_maxval, planetwide_minval, & planetwide_sumval use phyredem, only: physdem0, physdem1 use phyetat0_mod, only: phyetat0, tab_cntrl_mod use wstats_mod, only: callstats, wstats, mkstats use eofdump_mod, only: eofdump USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt USE mod_phys_lmdz_omp_data, ONLY: is_omp_master USE time_phylmdz_mod, ONLY: day_end #endif #ifdef CPP_XIOS use xios_output_mod, only: initialize_xios_output, & update_xios_timestep, & send_xios_field use wxios, only: wxios_context_init, xios_context_finalize #endif USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured, & regular_lonlat use ioipsl_getin_p_mod, only: getin_p use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, & subslope_dist,iflat,sky_slope, & major_slope,compute_meshgridavg, & ini_comslope_h use write_output_mod, only: write_output IMPLICIT NONE c======================================================================= c c subject: c -------- c c Organisation of the physical parametrisations of the LMD c martian atmospheric general circulation model. c c The GCM can be run without or with tracer transport c depending on the value of Logical "tracer" in file "callphys.def" c Tracers may be water vapor, ice OR chemical species OR dust particles c c SEE comments in initracer.F about numbering of tracer species... c c It includes: c c 1. Initialization: c 1.1 First call initializations c 1.2 Initialization for every call to physiq c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff. c 2. Compute radiative transfer tendencies c (longwave and shortwave) for CO2 and aerosols. c 3. Gravity wave and subgrid scale topography drag : c 4. Vertical diffusion (turbulent mixing): c 5. Convective adjustment c 6. Condensation and sublimation of carbon dioxide. c 7. TRACERS : c 7a. water, water ice, co2 ice (clouds) c 7b. call for photochemistry when tracers are chemical species c 7c. other scheme for tracer (dust) transport (lifting, sedimentation) c 7d. updates (CO2 pressure variations, surface budget) c 8. Contribution to tendencies due to thermosphere c 9. Surface and sub-surface temperature calculations c 10. Write outputs : c - "startfi", "histfi" (if it's time) c - Saving statistics (if "callstats = .true.") c - Dumping eof (if "calleofdump = .true.") c - Output any needed variables in "diagfi" c 11. Diagnostic: mass conservation of tracers c c author: c ------- c Frederic Hourdin 15/10/93 c Francois Forget 1994 c Christophe Hourdin 02/1997 c Subroutine completly rewritten by F.Forget (01/2000) c Introduction of the photochemical module: S. Lebonnois (11/2002) c Introduction of the thermosphere module: M. Angelats i Coll (2002) c Water ice clouds: Franck Montmessin (update 06/2003) c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) c Nb: See callradite.F for more information. c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization c c 10/16 J. Audouard: modifications for CO2 clouds scheme c arguments: c ---------- c c input: c ------ c ecri period (in dynamical timestep) to write output c ngrid Size of the horizontal grid. c All internal loops are performed on that grid. c nlayer Number of vertical layers. c nq Number of advected fields c firstcall True at the first call c lastcall True at the last call c pday Number of days counted from the North. Spring c equinoxe. c ptime Universal time (00 when downwards) INTEGER length PARAMETER (length=100) c Variables for the total dust for diagnostics REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust INTEGER iaer c local variables only used for diagnostic (output in file "diagfi" or "stats") c ----------------------------------------------------------------------------- REAL ps(ngrid), zt(ngrid,nlayer) REAL zu(ngrid,nlayer),zv(ngrid,nlayer) REAL zq(ngrid,nlayer,nq) REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid) REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid) character*2 str2 ! character*5 str5 real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer) real rdust(ngrid,nlayer) ! dust geometric mean radius (m) real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m) real rtopdust(ngrid,nlayer) ! topdust geometric mean radius (m) integer igmin, lmin logical tdiag real co2col(ngrid) ! CO2 column ! pplev and pplay are dynamical inputs and must not be modified in the physics. ! instead, use zplay and zplev : REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) ! REAL zstress(ngrid),cd real tmean, zlocal(nlayer) real rho(ngrid,nlayer) ! density real vmr(ngrid,nlayer) ! volume mixing ratio real rhopart(ngrid,nlayer) ! number density of a given species real colden(ngrid,nq) ! vertical column of tracers real mass(nq) ! global mass of tracers (g) REAL mtot(ngrid) ! Total mass of water vapor (kg/m2) REAL mstormdtot(ngrid) ! Total mass of stormdust tracer (kg/m2) REAL mdusttot(ngrid) ! Total mass of dust tracer (kg/m2) REAL icetot(ngrid) ! Total mass of water ice (kg/m2) REAL mtotco2(ngrid) ! Total mass of co2, including ice at the surface (kg/m2) REAL vaptotco2(ngrid) ! Total mass of co2 vapor (kg/m2) REAL icetotco2(ngrid) ! Total mass of co2 ice (kg/m2) REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2) REAL NccnCO2tot(ngrid) ! Total number of ccnCO2 (nbr/m2) REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2) REAL rave(ngrid) ! Mean water ice effective radius (m) REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1 REAL tauTES(ngrid) ! column optical depth at 825 cm-1 REAL Qabsice ! Water ice absorption coefficient REAL taucloudtes(ngrid) ! Cloud opacity at infrared ! reference wavelength using ! Qabs instead of Qext ! (direct comparison with TES) REAL mtotD(ngrid) ! Total mass of HDO vapor (kg/m2) REAL icetotD(ngrid) ! Total mass of HDO ice (kg/m2) REAL DoH_vap(ngrid,nlayer) !D/H ratio REAL DoH_ice(ngrid,nlayer) !D/H ratio REAL DoH_surf(ngrid) !D/H ratio surface REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s) REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s) REAL ndust(ngrid,nlayer) ! true n dust (kg/kg) REAL qdust(ngrid,nlayer) ! true q dust (kg/kg) REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg) REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg) c definition tendancies of stormdust tracers REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s) REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s) REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg) REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg) REAL wspeed(ngrid,nlayer+1) ! vertical velocity stormdust tracer REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer REAL dsodust(ngrid,nlayer) ! density scaled opacity for background dust REAL dsords(ngrid,nlayer) ! density scaled opacity for stormdust REAL dsotop(ngrid,nlayer) ! density scaled opacity for topdust c Test 1d/3d scavenging real h2otot(ngrid) real hdotot(ngrid) REAL satu(ngrid,nlayer) ! satu ratio for output REAL zqsat(ngrid,nlayer) ! saturation REAL satuco2(ngrid,nlayer) ! co2 satu ratio for output REAL zqsatco2(ngrid,nlayer) ! saturation co2 ! Added for new NLTE scheme real co2vmr_gcm(ngrid,nlayer) real n2vmr_gcm(ngrid,nlayer) real ovmr_gcm(ngrid,nlayer) real covmr_gcm(ngrid,nlayer) integer ierr_nlte real*8 varerr C Non-oro GW drag & Calcul of Brunt-Vaisala freq. (BV2) REAL ztetalev(ngrid,nlayer) real zdtetalev(ngrid,nlayer), zdzlev(ngrid,nlayer) REAL bv2(ngrid,nlayer) ! BV2 at zlev c Non-oro GW tendencies REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer) REAL d_t_hin(ngrid,nlayer) c Diagnostics 2D of gw_nonoro REAL zustrhi(ngrid), zvstrhi(ngrid) c Variables for PBL REAL zz1(ngrid) REAL lmax_th_out(ngrid) REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer) REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq) INTEGER lmax_th(ngrid),dimout,n_out,n CHARACTER(50) zstring REAL dtke_th(ngrid,nlayer+1) REAL zcdv(ngrid), zcdh(ngrid) REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out REAL u_out1(ngrid) REAL T_out1(ngrid) REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters] REAL tstar(ngrid) ! friction velocity and friction potential temp REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid) real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer) real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer) real qdusttotal0(ngrid),qdusttotal1(ngrid) c sub-grid scale water ice clouds (A. Pottier 2013) logical clearsky ! flux for the part without clouds real zdtswclf(ngrid,nlayer) real zdtlwclf(ngrid,nlayer) real fluxsurf_lwclf(ngrid) real fluxsurf_dn_swclf(ngrid,2),fluxsurf_up_swclf(ngrid,2) real fluxtop_lwclf(ngrid) real fluxtop_dn_swclf(ngrid,2),fluxtop_up_swclf(ngrid,2) real taucloudtesclf(ngrid) real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m) C test de conservation de la masse de CO2 REAL co2totA REAL co2totB REAL co2conservation c entrainment by mountain top dust flows above sub-grid scale topography REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons c when no startfi file is asked for init real alpha,lay1 ! coefficients for building layers integer iloop ! flags to trigger extra sanity checks logical,save :: check_physics_inputs=.false. logical,save :: check_physics_outputs=.false. !$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs) c Sub-grid scale slopes real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K] real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1] real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1] real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2] real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie integer :: islope real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting logical :: write_restart ! Variable for ice table REAL :: rhowater_surf(ngrid,nslope) ! Water density at the surface [kg/m^3] REAL :: rhowater_surf_sat(ngrid,nslope) ! Water density at the surface at saturation [kg/m^3] REAL :: rhowater_soil(ngrid,nsoilmx,nslope) ! Water density in soil layers [kg/m^3] REAL,PARAMETER :: alpha_clap_h2o = 28.9074 ! Coeff for Clapeyron law [/] REAL,PARAMETER :: beta_clap_h2o = -6143.7 ! Coeff for Clapeyron law [K] REAL :: pvap_surf(ngrid) ! Water vapor partial pressure in first layer [Pa] REAL,PARAMETER :: m_co2 = 44.01E-3 ! CO2 molecular mass [kg/mol] REAL,PARAMETER :: m_noco2 = 33.37E-3 ! Non condensible mol mass [kg/mol] REAL :: ztmp1,ztmp2 ! intermediate variables to compute the mean molar mass of the layer c======================================================================= pdq(:,:,:) = 0. c 1. Initialisation: c ----------------- c 1.1 Initialisation only at first call c --------------------------------------- IF (firstcall) THEN call getin_p("check_physics_inputs",check_physics_inputs) call getin_p("check_physics_outputs",check_physics_outputs) c variables set to 0 c ~~~~~~~~~~~~~~~~~~ aerosol(:,:,:)=0 dtrad(:,:)=0 #ifndef MESOSCALE fluxrad(:,:)=0 wstar(:)=0. #endif #ifdef CPP_XIOS ! Initialize XIOS context write(*,*) "physiq: call wxios_context_init" CALL wxios_context_init #endif c read startfi c ~~~~~~~~~~~~ #ifndef MESOSCALE ! GCM. Read netcdf initial physical parameters. CALL phyetat0 ("startfi.nc",0,0, & nsoilmx,ngrid,nlayer,nq,nqsoil, & day_ini,time_phys, & tsurf,tsoil,albedo,emis, & q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & watercap,perenial_co2ice, & def_slope,def_slope_mean,subslope_dist) ! Sky view: DO islope=1,nslope sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2. END DO ! Determine the 'flatest' slopes iflat = 1 DO islope=2,nslope IF(abs(def_slope_mean(islope)).lt. & abs(def_slope_mean(iflat)))THEN iflat = islope ENDIF ENDDO write(*,*)'Flat slope for islope = ',iflat write(*,*)'corresponding criterium = ',def_slope_mean(iflat) #else ! MESOSCALE. Supposedly everything is already set in modules. ! So we just check. And we fill day_ini write(*,*)"check: --- in physiq.F" write(*,*)"check: rad,cpp,g,r,rcp,daysec" write(*,*)rad,cpp,g,r,rcp,daysec write(*,*)'check: tsurf ',tsurf(1,:),tsurf(ngrid,:) write(*,*)'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:) write(*,*)'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx) write(*,*)'check: midlayer,layer ', mlayer(:),layer(:) write(*,*)'check: tracernames ', noms write(*,*)'check: emis ',emis(1,:),emis(ngrid,:) write(*,*)'check: q2 ',q2(1,1),q2(ngrid,nlayer+1) write(*,*)'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:) write(*,*)'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:) !!! day_ini = pday !!! a couple initializations (dummy for mesoscale) done in phyetat0 !!! --- maybe this should be done in update_inputs_physiq_mod tauscaling(:)=1.0 !! probably important totcloudfrac(:)=1.0 DO islope = 1,nslope albedo(:,1,islope)=albedodat(:) albedo(:,2,islope)=albedo(:,1,islope) inertiesoil(:,:,islope) = inertiedat(:,:) watercap(:,:)=0.0 ENDDO #endif #ifndef MESOSCALE if (.not.startphy_file) then ! starting without startfi.nc and with callsoil ! is not yet possible as soildepth default is not defined if (callsoil) then ! default mlayer distribution, following a power law: ! mlayer(k)=lay1*alpha**(k-1/2) lay1=2.e-4 alpha=2 do iloop=0,nsoilmx-1 mlayer(iloop)=lay1*(alpha**(iloop-0.5)) enddo lay1=sqrt(mlayer(0)*mlayer(1)) alpha=mlayer(1)/mlayer(0) do iloop=1,nsoilmx layer(iloop)=lay1*(alpha**(iloop-1)) enddo endif ! additionnal "academic" initialization of physics do islope = 1,nslope tsurf(:,islope)=pt(:,1) enddo write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!" do isoil=1,nsoilmx tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:) enddo write(*,*) "Physiq: initializing inertiedat !!" inertiedat(:,:)=400. inertiesoil(:,:,:)=400. write(*,*) "Physiq: initializing day_ini to pdat !" day_ini=pday endif #endif if (pday.ne.day_ini) then write(*,*) "PHYSIQ: ERROR: bad synchronization between ", & "physics and dynamics" write(*,*) "dynamics day [pday]: ",pday write(*,*) "physics day [day_ini]: ",day_ini call abort_physic("physiq","dynamics day /= physics day",1) endif write (*,*) 'In physiq day_ini =', day_ini c initialize tracers c ~~~~~~~~~~~~~~~~~~ CALL initracer(ngrid,nq,qsurf) c Initialize albedo and orbital calculation c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL surfini(ngrid,nslope,qsurf) CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) c initialize soil c ~~~~~~~~~~~~~~~ IF (callsoil) THEN c Thermal inertia feedback: IF (tifeedback) THEN DO islope = 1,nslope CALL soil_tifeedback(ngrid,nsoilmx, s qsurf(:,:,islope), s inertiesoil_tifeedback(:,:,islope)) ENDDO CALL soil(ngrid,nsoilmx,firstcall, s inertiesoil_tifeedback, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ELSE CALL soil(ngrid,nsoilmx,firstcall,inertiesoil, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ENDIF ! of IF (tifeedback) ELSE write(*,*) & 'PHYSIQ WARNING! Thermal conduction in the soil turned off' DO ig=1,ngrid capcal(ig,:)=1.e5 fluxgrd(ig,:)=0. ENDDO ENDIF icount=1 #ifndef MESOSCALE c Initialize thermospheric parameters c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (callthermos) then call fill_data_thermos call allocate_param_thermos(nlayer) call allocate_param_iono(nlayer,nreact) call param_read_e107 endif #endif c Initialize R and Cp as constant if (.not.callthermos .and. .not.photochem) then do l=1,nlayer do ig=1,ngrid rnew(ig,l)=r cpnew(ig,l)=cpp mmean(ig,l)=mugaz enddo enddo endif if(callnlte.and.nltemodel.eq.2) call nlte_setup if(callnirco2.and.nircorr.eq.1) call NIR_leedat IF (water.AND.(ngrid.NE.1)) THEN write(*,*)"physiq: water_param Surface water frost albedo:", . albedo_h2o_frost write(*,*)"physiq: water_param Surface watercap albedo:", . albedo_h2o_cap ENDIF #ifndef MESOSCALE if(ngrid.ne.1) then ! no need to compute slopes when in 1D; it is an input if (callslope) call getslopes(ngrid,phisfi) endif !1D if (ecritstart.GT.0) then call physdem0("restartfi.nc",longitude,latitude, & nsoilmx,ngrid,nlayer,nq, & ptimestep,pday,0.,cell_area, & albedodat,inertiedat,def_slope, & subslope_dist) else call physdem0("restartfi.nc",longitude,latitude, & nsoilmx,ngrid,nlayer,nq, & ptimestep,float(day_end),0.,cell_area, & albedodat,inertiedat,def_slope, & subslope_dist) endif c Initialize mountain mesh fraction for the entrainment by top flows param. c ~~~~~~~~~~~~~~~ if (topflows) call topmons_setup(ngrid) #endif #ifdef CPP_XIOS ! XIOS outputs write(*,*) "physiq firstcall: call initialize_xios_output" call initialize_xios_output(pday,ptime,ptimestep,daysec, & presnivs,pseudoalt,mlayer) #endif ENDIF ! (end of "if firstcall") if (check_physics_inputs) then ! Check the validity of input fields coming from the dynamics call check_physics_fields("begin physiq:",pt,pu,pv,pplev,pq) endif c --------------------------------------------------- c 1.2 Initializations done at every physical timestep: c --------------------------------------------------- c #ifdef CPP_XIOS ! update XIOS time/calendar call update_xios_timestep #endif c Initialize various variables c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ pdv(:,:)=0 pdu(:,:)=0 pdt(:,:)=0 pdq(:,:,:)=0 pdpsrf(:)=0 zflubid(:,:)=0 zdtsurf(:,:)=0 dqsurf(:,:,:)=0 dsodust(:,:)=0. dsords(:,:)=0. dsotop(:,:)=0. dwatercap(:,:)=0 call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf, & albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg) ! Dust scenario conversion coefficient from IRabs to VISext IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015 ! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true. #ifdef DUSTSTORM pq_tmp(:,:,:)=0 #endif igout=ngrid/2+1 zday=pday+ptime ! compute time, in sols (and fraction thereof) ! Compute local time at each grid point DO ig=1,ngrid local_time(ig)=modulo(1.+(zday-INT(zday)) & +(longitude_deg(ig)/15)/24,1.) ENDDO c Compute Solar Longitude (Ls) : c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (season) then call solarlong(zday,zls) else call solarlong(float(day_ini),zls) end if c Initialize pressure levels c ~~~~~~~~~~~~~~~~~~ zplev(:,:) = pplev(:,:) zplay(:,:) = pplay(:,:) ps(:) = pplev(:,1) c Compute geopotential at interlayers c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c ponderation des altitudes au niveau des couches en dp/p DO l=1,nlayer DO ig=1,ngrid zzlay(ig,l)=pphi(ig,l)/g ENDDO ENDDO DO ig=1,ngrid zzlev(ig,1)=0. zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km... ENDDO DO l=2,nlayer DO ig=1,ngrid z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l)) z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l)) zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) ENDDO ENDDO ! Potential temperature calculation not the same in physiq and dynamic c Compute potential temperature c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO l=1,nlayer DO ig=1,ngrid zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp zh(ig,l)=pt(ig,l)/zpopsk(ig,l) ENDDO ENDDO #ifndef MESOSCALE c----------------------------------------------------------------------- c 1.2.5 Compute mean mass, cp, and R c -------------------------------- if(photochem.or.callthermos) then call concentrations(ngrid,nlayer,nq, & zplay,pt,pdt,pq,pdq,ptimestep) endif #endif ! Compute vertical velocity (m/s) from vertical mass flux ! w = F / (rho*area) and rho = P/(r*T) ! but first linearly interpolate mass flux to mid-layers do l=1,nlayer-1 pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1)) enddo pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0 do l=1,nlayer pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) / & (pplay(1:ngrid,l)*cell_area(1:ngrid)) ! NB: here we use r and not rnew since this diagnostic comes ! from the dynamics enddo ! test for co2 conservation with co2 microphysics if (igcm_co2_ice.ne.0) then ! calculates the amount of co2 at the beginning of physics co2totA = 0. do ig=1,ngrid do l=1,nlayer co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g* & (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice) & +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep) end do do islope = 1,nslope co2totA = co2totA + qsurf(ig,igcm_co2,islope)* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) enddo end do else co2totA = 0. do ig=1,ngrid do l=1,nlayer co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g* & (pq(ig,l,igcm_co2) & +pdq(ig,l,igcm_co2)*ptimestep) end do do islope = 1,nslope co2totA = co2totA + qsurf(ig,igcm_co2,islope)* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) enddo end do endif ! of if (igcm_co2_ice.ne.0) c----------------------------------------------------------------------- c 2. Compute radiative tendencies : c------------------------------------ IF (callrad) THEN c Local Solar zenith angle c ~~~~~~~~~~~~~~~~~~~~~~~~ CALL orbite(zls,dist_sol,declin) IF (diurnal) THEN ztim1=SIN(declin) ztim2=COS(declin)*COS(2.*pi*(zday-.5)) ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) CALL solang(ngrid,sinlon,coslon,sinlat,coslat, & ztim1,ztim2,ztim3, mu0,fract) ELSE CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad) ENDIF ! of IF (diurnal) IF( MOD(icount-1,iradia).EQ.0) THEN c NLTE cooling from CO2 emission c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(callnlte) then if(nltemodel.eq.0.or.nltemodel.eq.1) then CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte) else if(nltemodel.eq.2) then co2vmr_gcm(1:ngrid,1:nlayer)= & pq(1:ngrid,1:nlayer,igcm_co2)* & mmean(1:ngrid,1:nlayer)/mmol(igcm_co2) n2vmr_gcm(1:ngrid,1:nlayer)= & pq(1:ngrid,1:nlayer,igcm_n2)* & mmean(1:ngrid,1:nlayer)/mmol(igcm_n2) covmr_gcm(1:ngrid,1:nlayer)= & pq(1:ngrid,1:nlayer,igcm_co)* & mmean(1:ngrid,1:nlayer)/mmol(igcm_co) ovmr_gcm(1:ngrid,1:nlayer)= & pq(1:ngrid,1:nlayer,igcm_o)* & mmean(1:ngrid,1:nlayer)/mmol(igcm_o) CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6, $ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, $ ovmr_gcm, zdtnlte,ierr_nlte,varerr ) if(ierr_nlte.gt.0) then write(*,*) $ 'WARNING: nlte_tcool output with error message', $ 'ierr_nlte=',ierr_nlte,'varerr=',varerr write(*,*)'I will continue anyway' endif zdtnlte(1:ngrid,1:nlayer)= & zdtnlte(1:ngrid,1:nlayer)/86400. endif ELSE zdtnlte(:,:)=0. ENDIF !end callnlte c Find number of layers for LTE radiation calculations IF(MOD(iphysiq*(icount-1),day_step).EQ.0) & CALL nlthermeq(ngrid,nlayer,zplev,zplay) c rocketstorm : compute dust storm mesh fraction IF (rdstorm) THEN CALL calcstormfract(ngrid,nlayer,nq,pq, & totstormfract) ENDIF c Note: Dustopacity.F has been transferred to callradite.F #ifdef DUSTSTORM !! specific case: save the quantity of dust before adding perturbation if (firstcall) then pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass) pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number) endif #endif c Call main radiative transfer scheme c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c Transfer through CO2 (except NIR CO2 absorption) c and aerosols (dust and water ice) ! callradite for background dust (out of the rdstorm fraction) clearatm=.true. !! callradite for background dust (out of the topflows fraction) nohmons=.true. c Radiative transfer c ------------------ ! callradite for the part with clouds clearsky=.false. ! part with clouds for both cases CLFvarying true/false CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq, & albedo_meshavg,emis_meshavg, & mu0,zplev,zplay,pt,tsurf_meshavg,fract,dist_sol,igout, & zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat), & fluxsurf_up_sw, & fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw, & tau_pref_scenario,tau_pref_gcm, & tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2, & qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract, & clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac) DO islope=1,nslope fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat) fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat) ENDDO ! case of sub-grid water ice clouds: callradite for the clear case IF (CLFvarying) THEN ! ---> PROBLEMS WITH ALLOCATED ARRAYS ! (temporary solution in callcorrk: do not deallocate ! if ! CLFvarying ...) ?? AP ?? clearsky=.true. CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq, & albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt, & tsurf_meshavg,fract, & dist_sol,igout,zdtlwclf,zdtswclf, & fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf, & fluxtop_lwclf,fluxtop_dn_swclf,fluxtop_up_swclf, & tau_pref_scenario,tau_pref_gcm,tau,aerosol, & dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, & taucloudtesclf,rdust, & rice,nuice,riceco2, nuiceco2, & qsurf_meshavg(:,igcm_co2), & rstormdust,rtopdust,totstormfract, & clearatm,dsords,dsotop, & nohmons,clearsky,totcloudfrac) clearsky = .false. ! just in case. ! Sum the fluxes and heating rates from cloudy/clear ! cases DO ig=1,ngrid tf_clf=totcloudfrac(ig) ntf_clf=1.-tf_clf DO islope=1,nslope fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig) & + tf_clf*fluxsurf_lw(ig,islope) fluxsurf_dn_sw(ig,1:2,islope) = & ntf_clf*fluxsurf_dn_swclf(ig,1:2) & + tf_clf*fluxsurf_dn_sw(ig,1:2,islope) ENDDO fluxsurf_up_sw(ig,1:2) = & ntf_clf*fluxsurf_up_swclf(ig,1:2) & + tf_clf*fluxsurf_up_sw(ig,1:2) fluxtop_lw(ig) = ntf_clf*fluxtop_lwclf(ig) & + tf_clf*fluxtop_lw(ig) fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2) & + tf_clf*fluxtop_dn_sw(ig,1:2) fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2) & + tf_clf*fluxtop_up_sw(ig,1:2) taucloudtes(ig) = ntf_clf*taucloudtesclf(ig) & + tf_clf*taucloudtes(ig) zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer) & + tf_clf*zdtlw(ig,1:nlayer) zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer) & + tf_clf*zdtsw(ig,1:nlayer) ENDDO ENDIF ! (CLFvarying) !============================================================================ #ifdef DUSTSTORM !! specific case: compute the added quantity of dust for perturbation if (firstcall) then pdq(1:ngrid,1:nlayer,igcm_dust_mass)= & pdq(1:ngrid,1:nlayer,igcm_dust_mass) & - pq_tmp(1:ngrid,1:nlayer,1) & + pq(1:ngrid,1:nlayer,igcm_dust_mass) pdq(1:ngrid,1:nlayer,igcm_dust_number)= & pdq(1:ngrid,1:nlayer,igcm_dust_number) & - pq_tmp(1:ngrid,1:nlayer,2) & + pq(1:ngrid,1:nlayer,igcm_dust_number) endif #endif c Outputs for basic check (middle of domain) c ------------------------------------------ write(*,'("Ls =",f11.6," check lat =",f10.6, & " lon =",f11.6)') & zls*180./pi,latitude(igout)*180/pi, & longitude(igout)*180/pi write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6, & " tau(",f4.0," Pa) =",f9.6)') & odpref,tau_pref_gcm(igout), & odpref,tau(igout,1)*odpref/zplev(igout,1) c --------------------------------------------------------- c Call slope parameterization for direct and scattered flux c --------------------------------------------------------- IF(callslope) THEN ! assume that in this case, nslope = 1 if(nslope.ne.1) then call abort_physic( & "physiq","callslope=true but nslope.ne.1",1) endif print *, 'Slope scheme is on and computing...' DO ig=1,ngrid sl_the = theta_sl(ig) IF (sl_the .ne. 0.) THEN ztim1=fluxsurf_dn_sw(ig,1,iflat) & +fluxsurf_dn_sw(ig,2,iflat) DO l=1,2 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15. sl_ra = pi*(1.0-sl_lct/12.) sl_lat = 180.*latitude(ig)/pi sl_tau = tau(ig,1) !il faudrait iaerdust(iaer) sl_alb = albedo(ig,l,iflat) sl_psi = psi_sl(ig) sl_fl0 = fluxsurf_dn_sw(ig,l,iflat) sl_di0 = 0. if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol sl_di0 = sl_di0/ztim1 sl_di0 = fluxsurf_dn_sw(ig,l,iflat)*sl_di0 endif ! you never know (roundup concern...) if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 !!!!!!!!!!!!!!!!!!!!!!!!!! CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, & sl_tau, sl_alb, sl_the, sl_psi, & sl_di0, sl_fl0, sl_flu ) !!!!!!!!!!!!!!!!!!!!!!!!!! fluxsurf_dn_sw(ig,l,1) = sl_flu ENDDO !!! compute correction on IR flux as well sky= (1.+cos(pi*theta_sl(ig)/180.))/2. fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky ENDIF ENDDO ELSE ! not calling subslope, nslope might be > 1 DO islope = 1,nslope sl_the=abs(def_slope_mean(islope)) IF (sl_the .gt. 1e-6) THEN IF(def_slope_mean(islope).ge.0.) THEN psi_sl(:) = 0. !Northward slope ELSE psi_sl(:) = 180. !Southward slope ENDIF DO ig=1,ngrid ztim1=fluxsurf_dn_sw(ig,1,islope) s +fluxsurf_dn_sw(ig,2,islope) DO l=1,2 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15. sl_ra = pi*(1.0-sl_lct/12.) sl_lat = 180.*latitude(ig)/pi sl_tau = tau(ig,1) !il faudrait iaerdust(iaer) sl_alb = albedo(ig,l,islope) sl_psi = psi_sl(ig) sl_fl0 = fluxsurf_dn_sw(ig,l,islope) sl_di0 = 0. if (mu0(ig) .gt. 0.) then sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol sl_di0 = sl_di0/ztim1 sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0 endif ! you never know (roundup concern...) if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 !!!!!!!!!!!!!!!!!!!!!!!!!! CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, & sl_tau, sl_alb, sl_the, sl_psi, & sl_di0, sl_fl0, sl_flu ) !!!!!!!!!!!!!!!!!!!!!!!!!! fluxsurf_dn_sw(ig,l,islope) = sl_flu ENDDO !!! compute correction on IR flux as well fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope) & *sky_slope(islope) ENDDO ENDIF ! sub grid is not flat ENDDO ! islope = 1,nslope ENDIF ! callslope c CO2 near infrared absorption c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ zdtnirco2(:,:)=0 if (callnirco2) then call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq, . mu0,fract,declin, zdtnirco2) endif c Radiative flux from the sky absorbed by the surface (W.m-2) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid DO islope = 1,nslope fluxrad_sky(ig,islope) = $ emis(ig,islope)*fluxsurf_lw(ig,islope) $ +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope)) $ +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope)) ENDDO ENDDO c Net atmospheric radiative heating rate (K.s-1) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(callnlte) THEN CALL blendrad(ngrid, nlayer, zplay, & zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad) ELSE DO l=1,nlayer DO ig=1,ngrid dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) & +zdtnirco2(ig,l) ENDDO ENDDO ENDIF ENDIF ! of if(mod(icount-1,iradia).eq.0) c Transformation of the radiative tendencies: c ------------------------------------------- c Net radiative surface flux (W.m-2) c ~~~~~~~~~~~~~~~~~~~~~~~~~~ c DO ig=1,ngrid DO islope = 1,nslope zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope) zplanck(ig)=emis(ig,islope)* $ stephan*zplanck(ig)*zplanck(ig) fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig) IF(callslope) THEN sky= (1.+cos(pi*theta_sl(ig)/180.))/2. fluxrad(ig,nslope)=fluxrad(ig,nslope)+ $ (1.-sky)*zplanck(ig) ELSE fluxrad(ig,islope)=fluxrad(ig,islope) + $ (1.-sky_slope(iflat))*emis(ig,iflat)* $ stephan*tsurf(ig,iflat)**4 ENDIF ENDDO ENDDO DO l=1,nlayer DO ig=1,ngrid pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) ENDDO ENDDO ENDIF ! of IF (callrad) c 3.1 Rocket dust storm c ------------------------------------------- IF (rdstorm) THEN clearatm=.false. pdqrds(:,:,:)=0. qdusttotal0(:)=0. qdusttotal1(:)=0. do ig=1,ngrid do l=1,nlayer zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential ! temperature tendency ! for diagnostics ! qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+ ! & pdq(ig,l,igcm_dust_mass)*ptimestep ! qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+ ! & pdq(ig,l,igcm_stormdust_mass)*ptimestep ! qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+ ! & qstormrds0(ig,l))*(zplev(ig,l)- ! & zplev(ig,l+1))/g enddo enddo ! call write_output('qdustrds0','qdust before rds', ! & 'kg/kg ',qdustrds0(:,:)) ! call write_output('qstormrds0','qstorm before rds', ! & 'kg/kg ',qstormrds0(:,:)) CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, & pq,pdq,pt,pdt,zplev,zplay,zzlev, & zzlay,zdtsw,zdtlw, c for radiative transfer & clearatm,icount,zday,zls, & tsurf_meshavg,qsurf_meshavg(:,igcm_co2), & igout,totstormfract,tauscaling, & dust_rad_adjust,IRtoVIScoef, & albedo_meshavg,emis_meshavg, c input sub-grid scale cloud & clearsky,totcloudfrac, c input sub-grid scale topography & nohmons, c output & pdqrds,wspeed,dsodust,dsords,dsotop, & tau_pref_scenario,tau_pref_gcm) c update the tendencies of both dust after vertical transport DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,igcm_stormdust_mass)= & pdq(ig,l,igcm_stormdust_mass)+ & pdqrds(ig,l,igcm_stormdust_mass) pdq(ig,l,igcm_stormdust_number)= & pdq(ig,l,igcm_stormdust_number)+ & pdqrds(ig,l,igcm_stormdust_number) pdq(ig,l,igcm_dust_mass)= & pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass) pdq(ig,l,igcm_dust_number)= & pdq(ig,l,igcm_dust_number)+ & pdqrds(ig,l,igcm_dust_number) ENDDO ENDDO do l=1,nlayer do ig=1,ngrid qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+ & pdq(ig,l,igcm_dust_mass)*ptimestep qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+ & pdq(ig,l,igcm_stormdust_mass)*ptimestep qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+ & qstormrds1(ig,l))*(zplev(ig,l)- & zplev(ig,l+1))/g enddo enddo c for diagnostics ! call write_output('qdustrds1','qdust after rds', ! & 'kg/kg ',qdustrds1(:,:)) ! call write_output('qstormrds1','qstorm after rds', ! & 'kg/kg ',qstormrds1(:,:)) ! ! call write_output('qdusttotal0','q sum before rds', ! & 'kg/m2 ',qdusttotal0(:)) ! call write_output('qdusttotal1','q sum after rds', ! & 'kg/m2 ',qdusttotal1(:)) ENDIF ! end of if(rdstorm) c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography c ------------------------------------------- IF (topflows) THEN clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains nohmons=.false. pdqtop(:,:,:)=0. CALL topmons(ngrid,nlayer,nq,ptime,ptimestep, & pq,pdq,pt,pdt,zplev,zplay,zzlev, & zzlay,zdtsw,zdtlw, & icount,zday,zls,tsurf(:,iflat), & qsurf_meshavg(:,igcm_co2), & igout,aerosol,tauscaling,dust_rad_adjust, & IRtoVIScoef,albedo_meshavg,emis_meshavg, & totstormfract,clearatm, & clearsky,totcloudfrac, & nohmons, & pdqtop,wtop,dsodust,dsords,dsotop, & tau_pref_scenario,tau_pref_gcm) c update the tendencies of both dust after vertical transport DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,igcm_topdust_mass)= & pdq(ig,l,igcm_topdust_mass)+ & pdqtop(ig,l,igcm_topdust_mass) pdq(ig,l,igcm_topdust_number)= & pdq(ig,l,igcm_topdust_number)+ & pdqtop(ig,l,igcm_topdust_number) pdq(ig,l,igcm_dust_mass)= & pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass) pdq(ig,l,igcm_dust_number)= & pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number) ENDDO ENDDO ENDIF ! end of if (topflows) c 3.3 Dust injection from the surface c ------------------------------------------- if (dustinjection.gt.0) then CALL compute_dtau(ngrid,nlayer, & zday,pplev,tau_pref_gcm, & ptimestep,local_time,IRtoVIScoef, & dustliftday) endif ! end of if (dustinjection.gt.0) c----------------------------------------------------------------------- c 4. Gravity wave and subgrid scale topography drag : c ------------------------------------------------- IF(calllott)THEN CALL calldrag_noro(ngrid,nlayer,ptimestep, & zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw) DO l=1,nlayer DO ig=1,ngrid pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) ENDDO ENDDO ENDIF c----------------------------------------------------------------------- c 5. Vertical diffusion (turbulent mixing): c ----------------------------------------- IF (calldifv) THEN DO ig=1,ngrid DO islope = 1,nslope zflubid(ig,islope)=fluxrad(ig,islope) & +fluxgrd(ig,islope) ENDDO ENDDO zdum1(:,:)=0 zdum2(:,:)=0 do l=1,nlayer do ig=1,ngrid zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) enddo enddo c ---------------------- c Treatment of a special case : using new surface layer (Richardson based) c without using the thermals in gcm and mesoscale can yield problems in c weakly unstable situations when winds are near to 0. For those cases, we add c a unit subgrid gustiness. Remember that thermals should be used we using the c Richardson based surface layer model. IF ( .not.calltherm . .and. callrichsl . .and. .not.turb_resolved) THEN DO ig=1, ngrid IF (zh(ig,1) .lt. tsurf_meshavg(ig)) THEN wstar(ig)=1. hfmax_th(ig)=0.2 ELSE wstar(ig)=0. hfmax_th(ig)=0. ENDIF ENDDO ENDIF c ---------------------- IF (tke_heat_flux .ne. 0.) THEN zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)* & (-alog(zplay(:,1)/zplev(:,1))) pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1) ENDIF c Calling vdif (Martian version WITH CO2 condensation) dwatercap_dif(:,:) = 0. zcdh(:) = 0. zcdv(:) = 0. CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk, $ ptimestep,capcal,lwrite, $ zplay,zplev,zzlay,zzlev,z0, $ pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,qsoil, $ zdum1,zdum2,zdh,pdq,zflubid, $ zdudif,zdvdif,zdhdif,zdtsdif,q2, & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th, & zcondicea_co2microp,sensibFlux, & dustliftday,local_time,watercap,dwatercap_dif) DO ig=1,ngrid zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:) dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:) ENDDO call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif, & albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp) IF (.not.turb_resolved) THEN DO l=1,nlayer DO ig=1,ngrid pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l) pdu(ig,l)=pdu(ig,l)+zdudif(ig,l) pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l) zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only ENDDO ENDDO DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) ENDDO ENDDO ENDDO DO iq=1, nq DO ig=1,ngrid dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:) ENDDO ENDDO ELSE write (*,*) '******************************************' write (*,*) '** LES mode: the difv part is only used to' write (*,*) '** - provide HFX and UST to the dynamics' write (*,*) '** - update TSURF' write (*,*) '******************************************' !! Specific treatment for lifting in turbulent-resolving mode (AC) IF (lifting .and. doubleq) THEN !! lifted dust is injected in the first layer. !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF. !! => lifted dust is not incremented before the sedimentation step. zdqdif(1:ngrid,1,1:nq)=0. zdqdif(1:ngrid,1,igcm_dust_number) = . -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number) zdqdif(1:ngrid,1,igcm_dust_mass) = . -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass) zdqdif(1:ngrid,2:nlayer,1:nq) = 0. DO iq=1, nq IF ((iq .ne. igcm_dust_mass) & .and. (iq .ne. igcm_dust_number)) THEN zdqsdif(:,iq,:)=0. ENDIF ENDDO ELSE zdqdif(1:ngrid,1:nlayer,1:nq) = 0. zdqsdif(1:ngrid,1:nq,1:nslope) = 0. ENDIF ENDIF ELSE DO ig=1,ngrid DO islope=1,nslope zdtsurf(ig,islope)=zdtsurf(ig,islope)+ s (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope) ENDDO ENDDO IF (turb_resolved) THEN write(*,*) 'Turbulent-resolving mode !' write(*,*) 'Please set calldifv to T in callphys.def' call abort_physic("physiq","turbulent-resolving mode",1) ENDIF ENDIF ! of IF (calldifv) c----------------------------------------------------------------------- c 6. Thermals : c ----------------------------- if(calltherm .and. .not.turb_resolved) then call calltherm_interface(ngrid,nlayer,nq,igcm_co2, $ zzlev,zzlay, $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, $ zplay,zplev,pphi,zpopsk, $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th, $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux) DO l=1,nlayer DO ig=1,ngrid pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l) pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l) pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l) q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep ENDDO ENDDO DO ig=1,ngrid q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep ENDDO DO iq=1,nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq) ENDDO ENDDO ENDDO lmax_th_out(:)=real(lmax_th(:)) else !of if calltherm lmax_th(:)=0 wstar(:)=0. hfmax_th(:)=0. lmax_th_out(:)=0. end if c----------------------------------------------------------------------- c 7. Dry convective adjustment: c ----------------------------- IF(calladj) THEN DO l=1,nlayer DO ig=1,ngrid zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ENDDO ENDDO zduadj(:,:)=0 zdvadj(:,:)=0 zdhadj(:,:)=0 CALL convadj(ngrid,nlayer,nq,ptimestep, $ zplay,zplev,zpopsk,lmax_th, $ pu,pv,zh,pq, $ pdu,pdv,zdh,pdq, $ zduadj,zdvadj,zdhadj, $ zdqadj) DO l=1,nlayer DO ig=1,ngrid pdu(ig,l)=pdu(ig,l)+zduadj(ig,l) pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l) pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l) zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only ENDDO ENDDO DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) ENDDO ENDDO ENDDO ENDIF ! of IF(calladj) c----------------------------------------------------- c 8. Non orographic Gravity waves : c ------------------------------------------------- IF (calllott_nonoro) THEN CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, & cpnew,rnew, & zplay, & zmax_th, ! max altitude reached by thermals (m) & pt, pu, pv, & pdt, pdu, pdv, & zustrhi,zvstrhi, & d_t_hin, d_u_hin, d_v_hin) ! Update tendencies pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer) & +d_t_hin(1:ngrid,1:nlayer) pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer) & +d_u_hin(1:ngrid,1:nlayer) pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer) & +d_v_hin(1:ngrid,1:nlayer) ENDIF ! of IF (calllott_nonoro) c----------------------------------------------------------------------- c 9. Specific parameterizations for tracers c: ----------------------------------------- c 9a. Water and ice c --------------- c --------------------------------------- c Water ice condensation in the atmosphere c ---------------------------------------- IF (water) THEN call watercloud(ngrid,nlayer,ptimestep, & zplev,zplay,pdpsrf,zzlay, pt,pdt, & pq,pdq,zdqcloud,zdtcloud, & nq,tau,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud,totcloudfrac) c Temperature variation due to latent heat release if (activice) then pdt(1:ngrid,1:nlayer) = & pdt(1:ngrid,1:nlayer) + & zdtcloud(1:ngrid,1:nlayer) endif ! increment water vapour and ice atmospheric tracers tendencies pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = & pdq(1:ngrid,1:nlayer,igcm_h2o_vap) + & zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap) pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = & pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + & zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice) if (hdo) then ! increment HDO vapour and ice atmospheric tracers tendencies pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = & pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + & zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap) pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = & pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + & zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice) endif !hdo ! increment dust and ccn masses and numbers ! We need to check that we have Nccn & Ndust > 0 ! This is due to single precision rounding problems if (microphys) then pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = & pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + & zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass) pdq(1:ngrid,1:nlayer,igcm_ccn_number) = & pdq(1:ngrid,1:nlayer,igcm_ccn_number) + & zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number) where (pq(:,:,igcm_ccn_mass) + & ptimestep*pdq(:,:,igcm_ccn_mass) < 0.) pdq(:,:,igcm_ccn_mass) = & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccn_number) = & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_ccn_number) + & ptimestep*pdq(:,:,igcm_ccn_number) < 0.) pdq(:,:,igcm_ccn_mass) = & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccn_number) = & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 end where endif if (scavenging) then pdq(1:ngrid,1:nlayer,igcm_dust_mass) = & pdq(1:ngrid,1:nlayer,igcm_dust_mass) + & zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass) pdq(1:ngrid,1:nlayer,igcm_dust_number) = & pdq(1:ngrid,1:nlayer,igcm_dust_number) + & zdqcloud(1:ngrid,1:nlayer,igcm_dust_number) where (pq(:,:,igcm_dust_mass) + & ptimestep*pdq(:,:,igcm_dust_mass) < 0.) pdq(:,:,igcm_dust_mass) = & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_dust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_dust_number) + & ptimestep*pdq(:,:,igcm_dust_number) < 0.) pdq(:,:,igcm_dust_mass) = & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_dust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where endif ! of if scavenging END IF ! of IF (water) c 9a bis. CO2 clouds (CL & JA) c --------------------------------------- c CO2 ice cloud condensation in the atmosphere c ---------------------------------------- c flag needed in callphys.def: c co2clouds=.true. is mandatory (default is .false.) c co2useh2o=.true. if you want to allow co2 condensation c on water ice particles c meteo_flux=.true. if you want to add a meteoritic c supply of CCN c CLFvaryingCO2=.true. if you want to have a sub-grid c temperature distribution c spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti c nuiceco2_sed=0.2 variance of the size distribution for the c sedimentation c nuiceco2_ref=0.2 variance of the size distribution for the c nucleation c imicroco2=50 micro-timestep is 1/50 of physical timestep zdqssed_co2(:) = 0. zdqssed_ccn(:,:) = 0. IF (co2clouds) THEN call co2cloud(ngrid,nlayer,ptimestep, & zplev,zplay,pdpsrf,zzlay,pt,pdt, & pq,pdq,zdqcloudco2,zdtcloudco2, & nq,tau,tauscaling,rdust,rice,riceco2,nuice, & rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2, & zdqssed_ccn,pdu,pu,zcondicea_co2microp) DO iq=1, nq DO ig=1,ngrid DO islope = 1,nslope dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+ & zdqssed_ccn(ig,iq)*cos(pi*def_slope_mean(islope)/180.) ENDDO !(islope) ENDDO ! (ig) ENDDO ! (iq)q) c Temperature variation due to latent heat release pdt(1:ngrid,1:nlayer) = & pdt(1:ngrid,1:nlayer) + & zdtcloudco2(1:ngrid,1:nlayer) ! increment dust and ccn masses and numbers ! We need to check that we have Nccn & Ndust > 0 ! This is due to single precision rounding problems ! increment dust tracers tendancies pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass) & + zdqcloudco2(:,:,igcm_dust_mass) pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number) & + zdqcloudco2(:,:,igcm_dust_number) pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2) & + zdqcloudco2(:,:,igcm_co2) pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice) & + zdqcloudco2(:,:,igcm_co2_ice) pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass) & + zdqcloudco2(:,:,igcm_ccnco2_mass) pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number) & + zdqcloudco2(:,:,igcm_ccnco2_number) if (meteo_flux) then pdq(:,:,igcm_ccnco2_meteor_mass) = & pdq(:,:,igcm_ccnco2_meteor_mass) + & zdqcloudco2(:,:,igcm_ccnco2_meteor_mass) pdq(:,:,igcm_ccnco2_meteor_number) = & pdq(:,:,igcm_ccnco2_meteor_number) & + zdqcloudco2(:,:,igcm_ccnco2_meteor_number) end if !Update water ice clouds values as well if (co2useh2o) then pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = & pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + & zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice) pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = & pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass) pdq(1:ngrid,1:nlayer,igcm_ccn_number) = & pdq(1:ngrid,1:nlayer,igcm_ccn_number) + & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number) pdq(:,:,igcm_ccnco2_h2o_mass_ice) = & pdq(:,:,igcm_ccnco2_h2o_mass_ice) + & zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) pdq(:,:,igcm_ccnco2_h2o_mass_ccn) = & pdq(:,:,igcm_ccnco2_h2o_mass_ccn) + & zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) pdq(:,:,igcm_ccnco2_h2o_number) = & pdq(:,:,igcm_ccnco2_h2o_number) + & zdqcloudco2(:,:,igcm_ccnco2_h2o_number) c Negative values? where (pq(:,:,igcm_ccn_mass) + & ptimestep*pdq(:,:,igcm_ccn_mass) < 0.) pdq(:,:,igcm_ccn_mass) = & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccn_number) = & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 end where c Negative values? where (pq(:,:,igcm_ccn_number) + & ptimestep*pdq(:,:,igcm_ccn_number) < 0.) pdq(:,:,igcm_ccn_mass) = & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccn_number) = & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_ccnco2_h2o_mass_ice) + & pq(:,:,igcm_ccnco2_h2o_mass_ccn) + & (pdq(:,:,igcm_ccnco2_h2o_mass_ice) + & pdq(:,:,igcm_ccnco2_h2o_mass_ccn) & )*ptimestep < 0.) pdq(:,:,igcm_ccnco2_h2o_mass_ice) = & - pq(:,:,igcm_ccnco2_h2o_mass_ice) & /ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_h2o_mass_ccn) = & - pq(:,:,igcm_ccnco2_h2o_mass_ccn) & /ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_h2o_number) = & - pq(:,:,igcm_ccnco2_h2o_number) & /ptimestep + 1.e-30 end where where (pq(:,:,igcm_ccnco2_h2o_number) + & (pdq(:,:,igcm_ccnco2_h2o_number) & )*ptimestep < 0.) pdq(:,:,igcm_ccnco2_h2o_mass_ice) = & - pq(:,:,igcm_ccnco2_h2o_mass_ice) & /ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_h2o_mass_ccn) = & - pq(:,:,igcm_ccnco2_h2o_mass_ccn) & /ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_h2o_number) = & - pq(:,:,igcm_ccnco2_h2o_number) & /ptimestep + 1.e-30 end where endif ! of if (co2useh2o) c Negative values? where (pq(:,:,igcm_ccnco2_mass) + & ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.) pdq(:,:,igcm_ccnco2_mass) = & - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_number) = & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_ccnco2_number) + & ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.) pdq(:,:,igcm_ccnco2_mass) = & - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_number) = & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 end where c Negative values? where (pq(:,:,igcm_dust_mass) + & ptimestep*pdq(:,:,igcm_dust_mass) < 0.) pdq(:,:,igcm_dust_mass) = & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_dust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_dust_number) + & ptimestep*pdq(:,:,igcm_dust_number) < 0.) pdq(:,:,igcm_dust_mass) = & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_dust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where if (meteo_flux) then where (pq(:,:,igcm_ccnco2_meteor_mass) + & ptimestep*pdq(:,:,igcm_ccnco2_meteor_mass) < 0.) pdq(:,:,igcm_ccnco2_meteor_mass) = & - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_meteor_number) = & - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_ccnco2_meteor_number) + & ptimestep*pdq(:,:,igcm_ccnco2_meteor_number) < 0.) pdq(:,:,igcm_ccnco2_meteor_mass) = & - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_ccnco2_meteor_number) = & - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30 end where end if END IF ! of IF (co2clouds) c 9b. Aerosol particles c ------------------- c ---------- c Dust devil : c ---------- IF(callddevil) then call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2, & zdqdev,zdqsdev) if (dustbin.ge.1) then do iq=1,nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq) ENDDO ENDDO enddo do iq=1,nq DO ig=1,ngrid DO islope = 1,nslope dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) + & zdqsdev(ig,iq)*cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO enddo endif ! of if (dustbin.ge.1) END IF ! of IF (callddevil) c ------------- c Sedimentation : acts also on water ice c ------------- IF (sedimentation) THEN zdqsed(1:ngrid,1:nlayer,1:nq)=0 zdqssed(1:ngrid,1:nq)=0 c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep c Zdqssed isn't call callsedim(ngrid,nlayer,ptimestep, & zplev,zzlev,zzlay,pt,pdt, & rdust,rstormdust,rtopdust, & rice,rsedcloud,rhocloud, & pq,pdq,zdqsed,zdqssed,nq, & tau,tauscaling) c Flux at the surface of co2 ice computed in co2cloud microtimestep IF (rdstorm) THEN c Storm dust cannot sediment to the surface DO ig=1,ngrid zdqsed(ig,1,igcm_stormdust_mass)= & zdqsed(ig,1,igcm_stormdust_mass)+ & zdqssed(ig,igcm_stormdust_mass) / & ((pplev(ig,1)-pplev(ig,2))/g) zdqsed(ig,1,igcm_stormdust_number)= & zdqsed(ig,1,igcm_stormdust_number)+ & zdqssed(ig,igcm_stormdust_number) / & ((pplev(ig,1)-pplev(ig,2))/g) zdqssed(ig,igcm_stormdust_mass)=0. zdqssed(ig,igcm_stormdust_number)=0. ENDDO ENDIF !rdstorm DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq) ENDDO ENDDO ENDDO DO iq=1, nq DO ig=1,ngrid DO islope = 1,nslope dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) + & zdqssed(ig,iq)*cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO ENDDO END IF ! of IF (sedimentation) c Add lifted dust to tendancies after sedimentation in the LES (AC) IF (turb_resolved) THEN DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) ENDDO ENDDO ENDDO DO iq=1, nq DO ig=1,ngrid dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:) ENDDO ENDDO ENDIF c c 9c. Chemical species c ------------------ #ifndef MESOSCALE c -------------- c photochemistry : c -------------- IF (photochem) then if (modulo(icount-1,ichemistry).eq.0) then ! compute chemistry every ichemistry physics step ! dust and ice surface area call surfacearea(ngrid, nlayer, naerkind, $ ptimestep, zplay, zzlay, $ pt, pq, pdq, nq, $ rdust, rice, tau, tauscaling, $ surfdust, surfice) ! call photochemistry DO ig = 1,ngrid qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig)) ENDDO call calchim(ngrid,nlayer,nq, & ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0, $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, $ zdqcloud,zdqscloud,tau(:,1), $ qsurf_tmp(:,igcm_co2), $ pu,pdu,pv,pdv,surfdust,surfice) endif ! of if (modulo(icount-1,ichemistry).eq.0) ! increment values of tracers: DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry ! tracers is zero anyways DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) ENDDO ENDDO ENDDO ! of DO iq=1,nq ! add condensation tendency for H2O2 if (igcm_h2o2.ne.0) then DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) & +zdqcloud(ig,l,igcm_h2o2) ENDDO ENDDO endif ! increment surface values of tracers: DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry ! tracers is zero anyways DO ig=1,ngrid DO islope = 1,nslope dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) + & zdqschim(ig,iq)*cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO ENDDO ! of DO iq=1,nq ! add condensation tendency for H2O2 if (igcm_h2o2.ne.0) then DO ig=1,ngrid DO islope = 1,nslope dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+ & zdqscloud(ig,igcm_h2o2)*cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO endif END IF ! of IF (photochem) #endif #ifndef MESOSCALE c----------------------------------------------------------------------- c 10. THERMOSPHERE CALCULATION c----------------------------------------------------------------------- if (callthermos) then call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol, $ mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay, & pt,pq,pu,pv,pdt,pdq, $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff, $ PhiEscH,PhiEscH2,PhiEscD) DO l=1,nlayer DO ig=1,ngrid dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l) pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) DO iq=1, nq pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) ENDDO ENDDO ENDDO endif ! of if (callthermos) #endif c----------------------------------------------------------------------- c 11. Carbon dioxide condensation-sublimation: c (should be the last atmospherical physical process to be computed) c ------------------------------------------- IF (tituscap) THEN !!! get the actual co2 seasonal cap from Titus observations CALL geticecover(ngrid, 180.*zls/pi, . 180.*longitude/pi, 180.*latitude/pi, . qsurf_tmp(:,igcm_co2) ) qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000. ENDIF IF (callcond) THEN zdtc(:,:) = 0. zdtsurfc(:,:) = 0. zduc(:,:) = 0. zdvc(:,:) = 0. zdqc(:,:,:) = 0. zdqsc(:,:,:) = 0. CALL co2condens(ngrid,nlayer,nq,nslope,ptimestep, $ capcal,zplay,zplev,tsurf,pt, $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, $ qsurf(:,igcm_co2,:),perenial_co2ice, $ albedo,emis,rdust, $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, $ fluxsurf_dn_sw,zls, $ zdqssed_co2,zcondicea_co2microp, & zdqsc) DO iq=1, nq DO ig=1,ngrid dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:) ENDDO ! (ig) ENDDO ! (iq) DO l=1,nlayer DO ig=1,ngrid pdt(ig,l)=pdt(ig,l)+zdtc(ig,l) pdv(ig,l)=pdv(ig,l)+zdvc(ig,l) pdu(ig,l)=pdu(ig,l)+zduc(ig,l) ENDDO ENDDO DO ig=1,ngrid zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:) ENDDO DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) ENDDO ENDDO ENDDO #ifndef MESOSCALE ! update surface pressure DO ig=1,ngrid ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep ENDDO ! update pressure levels DO l=1,nlayer DO ig=1,ngrid zplay(ig,l) = aps(l) + bps(l)*ps(ig) zplev(ig,l) = ap(l) + bp(l)*ps(ig) ENDDO ENDDO zplev(:,nlayer+1) = 0. ! update layers altitude DO l=2,nlayer DO ig=1,ngrid z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l)) z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l)) zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) ENDDO ENDDO #endif ENDIF ! of IF (callcond) c----------------------------------------------------------------------- c Updating tracer budget on surface c----------------------------------------------------------------------- DO iq=1, nq DO ig=1,ngrid DO islope = 1,nslope qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+ & ptimestep*dqsurf(ig,iq,islope) ENDDO ENDDO ! (ig) ENDDO ! (iq) c----------------------------------------------------------------------- c 12. Surface and sub-surface soil temperature c----------------------------------------------------------------------- c c c 12.1 Increment Surface temperature: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid DO islope = 1,nslope tsurf(ig,islope)=tsurf(ig,islope)+ & ptimestep*zdtsurf(ig,islope) ENDDO ENDDO c Prescribe a cold trap at south pole (except at high obliquity !!) c Temperature at the surface is set there to be the temperature c corresponding to equilibrium temperature between phases of CO2 IF (water) THEN !#ifndef MESOSCALE ! if (caps.and.(obliquit.lt.27.)) then => now done in co2condens ! NB: Updated surface pressure, at grid point 'ngrid', is ! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep ! tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095* ! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) ! tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*ps(ngrid))) ! endif !#endif c ------------------------------------------------------------- c Change of surface albedo in case of ground frost c everywhere except on the north permanent cap and in regions c covered by dry ice. c ALWAYS PLACE these lines after co2condens !!! c ------------------------------------------------------------- do ig=1,ngrid do islope = 1,nslope if ((qsurf(ig,igcm_co2,islope).eq.0).and. & (qsurf(ig,igcm_h2o_ice,islope) & .gt.frost_albedo_threshold)) then if ((watercaptag(ig)).and.(cst_cap_albedo)) then albedo(ig,1,islope) = albedo_h2o_cap albedo(ig,2,islope) = albedo_h2o_cap else albedo(ig,1,islope) = albedo_h2o_frost albedo(ig,2,islope) = albedo_h2o_frost endif !((watercaptag(ig)).and.(cst_cap_albedo)) then c write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice) c write(*,*) "physiq.F frost :" c & ,latitude(ig)*180./pi, longitude(ig)*180./pi endif enddo ! islope enddo ! of do ig=1,ngrid ENDIF ! of IF (water) c c 12.2 Compute soil temperatures and subsurface heat flux: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (callsoil) THEN c Thermal inertia feedback IF (tifeedback) THEN DO islope = 1,nslope CALL soil_tifeedback(ngrid,nsoilmx, s qsurf(:,:,islope), s inertiesoil_tifeedback(:,:,islope)) ENDDO CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ELSE CALL soil(ngrid,nsoilmx,.false.,inertiesoil, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ENDIF ENDIF c To avoid negative values IF (rdstorm) THEN where (pq(:,:,igcm_stormdust_mass) + & ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.) pdq(:,:,igcm_stormdust_mass) = & - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_stormdust_number) = & - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_stormdust_number) + & ptimestep*pdq(:,:,igcm_stormdust_number) < 0.) pdq(:,:,igcm_stormdust_mass) = & - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_stormdust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_dust_mass) + & ptimestep*pdq(:,:,igcm_dust_mass) < 0.) pdq(:,:,igcm_dust_mass) = & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_dust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where where (pq(:,:,igcm_dust_number) + & ptimestep*pdq(:,:,igcm_dust_number) < 0.) pdq(:,:,igcm_dust_mass) = & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 pdq(:,:,igcm_dust_number) = & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 end where ENDIF !(rdstorm) c----------------------------------------------------------------------- c J. Naar : Surface and sub-surface water ice c----------------------------------------------------------------------- c c c Increment Watercap (surface h2o reservoirs): c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid DO islope = 1,nslope watercap(ig,islope)=watercap(ig,islope)+ s ptimestep*dwatercap(ig,islope) ENDDO ENDDO IF (refill_watercap) THEN DO ig=1,ngrid DO islope = 1,nslope if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope) & .gt.frost_metam_threshold)) then watercap(ig,islope)=watercap(ig,islope) & +qsurf(ig,igcm_h2o_ice,islope) & - frost_metam_threshold qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold endif ! (watercaptag(ig).and. ENDDO ENDDO ENDIF ! (refill_watercap) THEN c----------------------------------------------------------------------- c 13. Write output files c ---------------------- call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf, & albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg) c ------------------------------- c Dynamical fields incrementation c ------------------------------- c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics) ! temperature, zonal and meridional wind DO l=1,nlayer DO ig=1,ngrid zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep ENDDO ENDDO ! tracers DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep ENDDO ENDDO ENDDO ! Density DO l=1,nlayer DO ig=1,ngrid rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l)) ENDDO ENDDO ! Potential Temperature DO ig=1,ngrid DO l=1,nlayer zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp ENDDO ENDDO c Compute surface stress : (NB: z0 is a common in surfdat.h) c DO ig=1,ngrid c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2 c zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2) c ENDDO c Sum of fluxes in solar spectral bands (for output only) fluxtop_dn_sw_tot(1:ngrid)=fluxtop_dn_sw(1:ngrid,1) + & fluxtop_dn_sw(1:ngrid,2) fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) + & fluxtop_up_sw(1:ngrid,2) fluxsurf_dn_sw_tot(1:ngrid,1:nslope)= & fluxsurf_dn_sw(1:ngrid,1,1:nslope) + & fluxsurf_dn_sw(1:ngrid,2,1:nslope) fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) + & fluxsurf_up_sw(1:ngrid,2) c ******* TEST ****************************************************** ztim1 = 999 DO l=1,nlayer DO ig=1,ngrid if (pt(ig,l).lt.ztim1) then ztim1 = pt(ig,l) igmin = ig lmin = l end if ENDDO ENDDO if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then write(*,*) 'PHYSIQ: stability WARNING :' write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin), & 'ig l =', igmin, lmin end if c ******************************************************************* c --------------------- c Outputs to the screen c --------------------- IF (lwrite) THEN write(*,*)'Global diagnostics for the physics' write(*,*)'Variables and their increments x and dx/dt * dt' WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps' WRITE(*,'(2f10.5,2f15.5)') s tsurf(igout,:),zdtsurf(igout,:)*ptimestep, s zplev(igout,1),pdpsrf(igout)*ptimestep WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT' WRITE(*,'(i4,6f10.5)') (l, s pu(igout,l),pdu(igout,l)*ptimestep, s pv(igout,l),pdv(igout,l)*ptimestep, s pt(igout,l),pdt(igout,l)*ptimestep, s l=1,nlayer) ENDIF ! of IF (lwrite) c ---------------------------------------------------------- c ---------------------------------------------------------- c INTERPOLATIONS IN THE SURFACE-LAYER c ---------------------------------------------------------- c ---------------------------------------------------------- n_out=0 ! number of elements in the z_out array. ! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set ! to 5 IF (n_out .ne. 0) THEN IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out)) IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out)) IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out)) z_out(:)=[3.,2.,1.,0.5,0.1] u_out(:,:)=0. T_out(:,:)=0. call pbl_parameters(ngrid,nlayer,ps,zplay,z0, & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf(:,iflat), & zh,z_out,n_out,T_out,u_out,ustar,tstar,L_mo,vhf,vvv) ! pourquoi ustar recalcule ici? fait dans vdifc. #ifndef MESOSCALE DO n=1,n_out write(zstring, '(F8.6)') z_out(n) call write_output('T_out_'//trim(zstring), & 'potential temperature at z_out','K',T_out(:,n)) call write_output('u_out_'//trim(zstring), & 'horizontal velocity norm at z_out','m/s',u_out(:,n)) ENDDO call write_output('u_star', & 'friction velocity','m/s',ustar) call write_output('teta_star', & 'friction potential temperature','K',tstar) call write_output('vvv', & 'Vertical velocity variance at zout','m',vvv) call write_output('vhf', & 'Vertical heat flux at zout','m',vhf) #else T_out1(:)=T_out(:,1) u_out1(:)=u_out(:,1) #endif ENDIF c ---------------------------------------------------------- c ---------------------------------------------------------- c END OF SURFACE LAYER INTERPOLATIONS c ---------------------------------------------------------- c ---------------------------------------------------------- #ifndef MESOSCALE c ------------------------------------------------------------------- c Writing NetCDF file "RESTARTFI" at the end of the run c ------------------------------------------------------------------- c Note: 'restartfi' is stored just before dynamics are stored c in 'restart'. Between now and the writting of 'restart', c there will have been the itau=itau+1 instruction and c a reset of 'time' (lastacll = .true. when itau+1= itaufin) c thus we store for time=time+dtvr ! default: not writing a restart file at this time step write_restart=.false. IF (ecritstart.GT.0) THEN ! For when we store multiple time steps in the restart file IF (MODULO(icount*iphysiq,ecritstart).EQ.0) THEN write_restart=.true. ENDIF ENDIF IF (lastcall) THEN ! Always write a restart at the end of the simulation write_restart=.true. ENDIF IF (write_restart) THEN IF (grid_type==unstructured) THEN !IF DYNAMICO ! When running Dynamico, no need to add a dynamics time step to ztime_fin IF (ptime.LE. 1.E-10) THEN ! Residual ptime occurs with Dynamico ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec) . - day_ini - time_phys ELSE ztime_fin = pday + ptime !+ ptimestep/(float(iphysiq)*daysec) . - day_ini - time_phys ENDIF if (ecritstart==0) then ztime_fin = ztime_fin-(day_end-day_ini) endif ELSE ! IF LMDZ if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change ztime_fin = pday - day_ini + ptime & + ptimestep/(float(iphysiq)*daysec) else !IF ONE RESTART final time in top of day_end ztime_fin = pday - day_ini-(day_end-day_ini) & + ptime + ptimestep/(float(iphysiq)*daysec) endif ENDIF ! of IF (grid_type==unstructured) write(*,'(A,I7,A,F12.5)') . 'PHYSIQ: writing in restartfi ; icount=', . icount,' date=',ztime_fin call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, . ptimestep,ztime_fin, . tsurf,tsoil,inertiesoil,albedo, . emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, . wstar,watercap,perenial_co2ice) ENDIF ! of IF (write_restart) #endif c IF (ngrid.NE.1) then c ------------------------------------------------------------------- c Calculation of diagnostic variables written in both stats and c diagfi files c ------------------------------------------------------------------- if (water) then do ig=1,ngrid if(mu0(ig).le.0.01) then fluxsurf_dir_dn_sw(ig) = 0. else fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_sol*mu0(ig)* & exp(-(tau(ig,iaer_dust_doubleq)+ & tau(ig,iaer_h2o_ice))/mu0(ig)) endif enddo endif ! Density-scaled opacities do ig=1,ngrid dsodust(ig,:) = & dsodust(ig,:)*tauscaling(ig) dsords(ig,:) = & dsords(ig,:)*tauscaling(ig) dsotop(ig,:) = & dsotop(ig,:)*tauscaling(ig) enddo if(doubleq) then do ig=1,ngrid IF (sedimentation) THEN dqdustsurf(ig) = & zdqssed(ig,igcm_dust_mass)*tauscaling(ig) dndustsurf(ig) = & zdqssed(ig,igcm_dust_number)*tauscaling(ig) ENDIF ndust(ig,:) = & zq(ig,:,igcm_dust_number)*tauscaling(ig) qdust(ig,:) = & zq(ig,:,igcm_dust_mass)*tauscaling(ig) enddo if (scavenging) then do ig=1,ngrid IF (sedimentation) THEN dqdustsurf(ig) = dqdustsurf(ig) + & zdqssed(ig,igcm_ccn_mass)*tauscaling(ig) dndustsurf(ig) = dndustsurf(ig) + & zdqssed(ig,igcm_ccn_number)*tauscaling(ig) ENDIF nccn(ig,:) = & zq(ig,:,igcm_ccn_number)*tauscaling(ig) qccn(ig,:) = & zq(ig,:,igcm_ccn_mass)*tauscaling(ig) enddo endif endif ! of (doubleq) if (rdstorm) then ! diagnostics of stormdust tendancies for 1D and 3D mstormdtot(:)=0 mdusttot(:)=0 qdusttotal(:,:)=0 do ig=1,ngrid rdsdqdustsurf(ig) = & zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig) rdsdndustsurf(ig) = & zdqssed(ig,igcm_stormdust_number)*tauscaling(ig) rdsndust(ig,:) = & pq(ig,:,igcm_stormdust_number)*tauscaling(ig) rdsqdust(ig,:) = & pq(ig,:,igcm_stormdust_mass)*tauscaling(ig) do l=1,nlayer mstormdtot(ig) = mstormdtot(ig) + & zq(ig,l,igcm_stormdust_mass) * & (zplev(ig,l) - zplev(ig,l+1)) / g mdusttot(ig) = mdusttot(ig) + & zq(ig,l,igcm_dust_mass) * & (zplev(ig,l) - zplev(ig,l+1)) / g qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust enddo enddo endif !(rdstorm) if (water) then mtot(:)=0 icetot(:)=0 rave(:)=0 tauTES(:)=0 IF (hdo) then mtotD(:)=0 icetotD(:)=0 ENDIF !hdo do ig=1,ngrid do l=1,nlayer mtot(ig) = mtot(ig) + & zq(ig,l,igcm_h2o_vap) * & (zplev(ig,l) - zplev(ig,l+1)) / g icetot(ig) = icetot(ig) + & zq(ig,l,igcm_h2o_ice) * & (zplev(ig,l) - zplev(ig,l+1)) / g IF (hdo) then mtotD(ig) = mtotD(ig) + & zq(ig,l,igcm_hdo_vap) * & (zplev(ig,l) - zplev(ig,l+1)) / g icetotD(ig) = icetotD(ig) + & zq(ig,l,igcm_hdo_ice) * & (zplev(ig,l) - zplev(ig,l+1)) / g ENDIF !hdo c Computing abs optical depth at 825 cm-1 in each c layer to simulate NEW TES retrieval Qabsice = min( & max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2 & ) opTES(ig,l)= 0.75 * Qabsice * & zq(ig,l,igcm_h2o_ice) * & (zplev(ig,l) - zplev(ig,l+1)) / g & / (rho_ice * rice(ig,l) * (1.+nuice_ref)) tauTES(ig)=tauTES(ig)+ opTES(ig,l) enddo c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0. enddo call watersat(ngrid*nlayer,zt,zplay,zqsat) satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:) if (scavenging) then Nccntot(:)= 0 Mccntot(:)= 0 rave(:)=0 do ig=1,ngrid do l=1,nlayer Nccntot(ig) = Nccntot(ig) + & zq(ig,l,igcm_ccn_number)*tauscaling(ig) & *(zplev(ig,l) - zplev(ig,l+1)) / g Mccntot(ig) = Mccntot(ig) + & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) & *(zplev(ig,l) - zplev(ig,l+1)) / g cccc Column integrated effective ice radius cccc is weighted by total ice surface area (BETTER than total ice mass) rave(ig) = rave(ig) + & tauscaling(ig) * & zq(ig,l,igcm_ccn_number) * & (zplev(ig,l) - zplev(ig,l+1)) / g * & rice(ig,l) * rice(ig,l)* (1.+nuice_ref) enddo rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75 & /max(pi*rave(ig),1.e-30) ! surface weight if (icetot(ig)*1e3.lt.0.01) rave(ig)=0. enddo else ! of if (scavenging) rave(:)=0 do ig=1,ngrid do l=1,nlayer rave(ig) = rave(ig) + & zq(ig,l,igcm_h2o_ice) * & (zplev(ig,l) - zplev(ig,l+1)) / g * & rice(ig,l) * (1.+nuice_ref) enddo rave(ig) = max(rave(ig) / & max(icetot(ig),1.e-30),1.e-30) ! mass weight enddo endif ! of if (scavenging) !Alternative A. Pottier weighting rave2(:) = 0. totrave2(:) = 0. do ig=1,ngrid do l=1,nlayer rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l) totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice) end do rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30) end do endif ! of if (water) if (co2clouds) then mtotco2(1:ngrid) = 0. icetotco2(1:ngrid) = 0. vaptotco2(1:ngrid) = 0. do ig=1,ngrid do l=1,nlayer vaptotco2(ig) = vaptotco2(ig) + & zq(ig,l,igcm_co2) * & (zplev(ig,l) - zplev(ig,l+1)) / g icetotco2(ig) = icetot(ig) + & zq(ig,l,igcm_co2_ice) * & (zplev(ig,l) - zplev(ig,l+1)) / g end do mtotco2(ig) = icetotco2(ig) + vaptotco2(ig) end do end if #ifndef MESOSCALE c ----------------------------------------------------------------- c WSTATS: Saving statistics c ----------------------------------------------------------------- c ("stats" stores and accumulates key variables in file "stats.nc" c which can later be used to make the statistic files of the run: c if flag "callstats" from callphys.def is .true.) call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) call wstats(ngrid,"tsurf","Surface temperature","K",2 & ,tsurf(:,iflat)) call wstats(ngrid,"co2ice","CO2 ice cover", & "kg.m-2",2,qsurf(:,igcm_co2,iflat)) call wstats(ngrid,"watercap","H2O ice cover", & "kg.m-2",2,watercap(:,iflat)) call wstats(ngrid,"tau_pref_scenario", & "prescribed visible dod at 610 Pa","NU", & 2,tau_pref_scenario) call wstats(ngrid,"tau_pref_gcm", & "visible dod at 610 Pa in the GCM","NU", & 2,tau_pref_gcm) call wstats(ngrid,"fluxsurf_lw", & "Thermal IR radiative flux to surface","W.m-2",2, & fluxsurf_lw(:,iflat)) call wstats(ngrid,"fluxsurf_dn_sw", & "Incoming Solar radiative flux to surface","W.m-2",2, & fluxsurf_dn_sw_tot(:,iflat)) call wstats(ngrid,"fluxsurf_up_sw", & "Reflected Solar radiative flux from surface","W.m-2",2, & fluxsurf_up_sw_tot) call wstats(ngrid,"fluxtop_lw", & "Thermal IR radiative flux to space","W.m-2",2, & fluxtop_lw) call wstats(ngrid,"fluxtop_dn_sw", & "Incoming Solar radiative flux from space","W.m-2",2, & fluxtop_dn_sw_tot) call wstats(ngrid,"fluxtop_up_sw", & "Outgoing Solar radiative flux to space","W.m-2",2, & fluxtop_up_sw_tot) call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) call wstats(ngrid,"v","Meridional (North-South) wind", & "m.s-1",3,zv) call wstats(ngrid,"w","Vertical (down-up) wind", & "m.s-1",3,pw) call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho) call wstats(ngrid,"pressure","Pressure","Pa",3,zplay) call wstats(ngrid,"q2", & "Boundary layer eddy kinetic energy", & "m2.s-2",3,q2) call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, & emis(:,iflat)) call wstats(ngrid,"fluxsurf_dir_dn_sw", & "Direct incoming SW flux at surface", & "W.m-2",2,fluxsurf_dir_dn_sw) if (calltherm) then call wstats(ngrid,"zmax_th","Height of thermals", & "m",2,zmax_th) call wstats(ngrid,"hfmax_th","Max thermals heat flux", & "K.m/s",2,hfmax_th) call wstats(ngrid,"wstar", & "Max vertical velocity in thermals", & "m/s",2,wstar) endif if (water) then vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap) & *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) call wstats(ngrid,"vmr_h2ovap", & "H2O vapor volume mixing ratio","mol/mol", & 3,vmr) vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice) & *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice) call wstats(ngrid,"vmr_h2oice", & "H2O ice volume mixing ratio","mol/mol", & 3,vmr) ! also store vmr_ice*rice for better diagnostics of rice vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)* & rice(1:ngrid,1:nlayer) call wstats(ngrid,"vmr_h2oice_rice", & "H2O ice mixing ratio times ice particule size", & "(mol/mol)*m", & 3,vmr) vmr=zqsat(1:ngrid,1:nlayer) & *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) call wstats(ngrid,"vmr_h2osat", & "saturation volume mixing ratio","mol/mol", & 3,vmr) call wstats(ngrid,"h2o_ice_s", & "surface h2o_ice","kg/m2", & 2,qsurf(1,igcm_h2o_ice,iflat)) call wstats(ngrid,'albedo', & 'albedo', & '',2,albedo(1,1,iflat)) call wstats(ngrid,"mtot", & "total mass of water vapor","kg/m2", & 2,mtot) call wstats(ngrid,"icetot", & "total mass of water ice","kg/m2", & 2,icetot) call wstats(ngrid,"reffice", & "Mean reff","m", & 2,rave) call wstats(ngrid,"Nccntot", & "condensation nuclei","Nbr/m2", & 2,Nccntot) call wstats(ngrid,"Mccntot", & "condensation nuclei mass","kg/m2", & 2,Mccntot) call wstats(ngrid,"rice", & "Ice particle size","m", & 3,rice) if (.not.activice) then call wstats(ngrid,"tauTESap", & "tau abs 825 cm-1","", & 2,tauTES) else call wstats(ngrid,'tauTES', & 'tau abs 825 cm-1', & '',2,taucloudtes) endif endif ! of if (water) if (co2clouds) then call wstats(ngrid,"mtotco2", & "total mass atm of co2","kg/m2", & 2,mtotco2) call wstats(ngrid,"icetotco2", & "total mass atm of co2 ice","kg/m2", & 2,icetotco2) call wstats(ngrid,"vaptotco2", & "total mass atm of co2 vapor","kg/m2", & 2,icetotco2) end if if (dustbin.ne.0) then call wstats(ngrid,'tau','taudust','SI',2,tau(1,1)) if (doubleq) then call wstats(ngrid,'dqsdust', & 'deposited surface dust mass', & 'kg.m-2.s-1',2,dqdustsurf) call wstats(ngrid,'dqndust', & 'deposited surface dust number', & 'number.m-2.s-1',2,dndustsurf) call wstats(ngrid,'reffdust','reffdust', & 'm',3,rdust*ref_r0) call wstats(ngrid,'dustq','Dust mass mr', & 'kg/kg',3,qdust) call wstats(ngrid,'dustN','Dust number', & 'part/kg',3,ndust) if (rdstorm) then call wstats(ngrid,'reffstormdust','reffdust', & 'm',3,rstormdust*ref_r0) call wstats(ngrid,'rdsdustq','Dust mass mr', & 'kg/kg',3,rdsqdust) call wstats(ngrid,'rdsdustN','Dust number', & 'part/kg',3,rdsndust) end if else do iq=1,dustbin write(str2(1:2),'(i2.2)') iq call wstats(ngrid,'q'//str2,'mix. ratio', & 'kg/kg',3,zq(1,1,iq)) call wstats(ngrid,'qsurf'//str2,'qsurf', & 'kg.m-2',2,qsurf(1,iq,iflat)) end do endif ! (doubleq) if (scavenging) then call wstats(ngrid,'ccnq','CCN mass mr', & 'kg/kg',3,qccn) call wstats(ngrid,'ccnN','CCN number', & 'part/kg',3,nccn) endif ! (scavenging) endif ! (dustbin.ne.0) if (photochem) then do iq=1,nq if (noms(iq) .ne. "dust_mass" .and. $ noms(iq) .ne. "dust_number" .and. $ noms(iq) .ne. "ccn_mass" .and. $ noms(iq) .ne. "ccn_number" .and. $ noms(iq) .ne. "ccnco2_mass" .and. $ noms(iq) .ne. "ccnco2_number" .and. $ noms(iq) .ne. "stormdust_mass" .and. $ noms(iq) .ne. "stormdust_number" .and. $ noms(iq) .ne. "topdust_mass" .and. $ noms(iq) .ne. "topdust_number") then ! volume mixing ratio vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq) & *mmean(1:ngrid,1:nlayer)/mmol(iq) call wstats(ngrid,"vmr_"//trim(noms(iq)), $ "Volume mixing ratio","mol/mol",3,vmr) if ((noms(iq).eq."o") $ .or. (noms(iq).eq."co2") $ .or. (noms(iq).eq."o3") $ .or. (noms(iq).eq."ar") $ .or. (noms(iq).eq."o2") $ .or. (noms(iq).eq."h2o_vap") ) then call write_output("vmr_"//trim(noms(iq)), $ "Volume mixing ratio","mol/mol",vmr(:,:)) end if ! number density (molecule.cm-3) rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq) & *rho(1:ngrid,1:nlayer)*n_avog/ & (1000*mmol(iq)) call wstats(ngrid,"num_"//trim(noms(iq)), $ "Number density","cm-3",3,rhopart) call write_output("num_"//trim(noms(iq)), $ "Number density","cm-3",rhopart(:,:)) ! vertical column (molecule.cm-2) do ig = 1,ngrid colden(ig,iq) = 0. end do do l=1,nlayer do ig=1,ngrid colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) $ *(zplev(ig,l)-zplev(ig,l+1)) $ *6.022e22/(mmol(iq)*g) end do end do call wstats(ngrid,"c_"//trim(noms(iq)), $ "column","mol cm-2",2,colden(1,iq)) call write_output("c_"//trim(noms(iq)), $ "column","mol cm-2",colden(:,iq)) ! global mass (g) call planetwide_sumval(colden(:,iq)/6.022e23 $ *mmol(iq)*1.e4*cell_area(:),mass(iq)) call write_output("mass_"//trim(noms(iq)), $ "global mass","g",mass(iq)) end if ! of if (noms(iq) .ne. "dust_mass" ...) end do ! of do iq=1,nq end if ! of if (photochem) IF(lastcall.and.callstats) THEN write (*,*) "Writing stats..." call mkstats(ierr) ENDIF c (Store EOF for Mars Climate database software) IF (calleofdump) THEN CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps) ENDIF #endif !endif of ifndef MESOSCALE #ifdef MESOSCALE !! see comm_wrf. !! not needed when an array is already in a shared module. !! --> example : hfmax_th, zmax_th CALL allocate_comm_wrf(ngrid,nlayer) !state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s" comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer) !state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s" comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer) !state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2" comm_SWDOWNZ(1:ngrid) = fluxsurf_dn_sw_tot(1:ngrid) !state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" "" comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid) !state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m" comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer) !state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2" IF (igcm_dust_mass .ne. 0) THEN comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass) ELSE comm_QSURFDUST(1:ngrid) = 0. ENDIF !state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic" comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice !state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2" comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice !state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm" IF (igcm_h2o_ice .ne. 0) THEN comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6 . * zq(1:ngrid,1:nlayer,igcm_h2o_ice) . * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice) ELSE comm_VMR_ICE(1:ngrid,1:nlayer) = 0. ENDIF !state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" "" if (activice) then comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid) else comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid) endif !state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m" comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer) !! calculate sensible heat flux in W/m2 for outputs !! -- the one computed in vdifc is not the real one !! -- vdifc must have been called if (.not.callrichsl) then sensibFlux(1:ngrid) = zflubid(1:ngrid) . - capcal(1:ngrid)*zdtsdif(1:ngrid) else sensibFlux(1:ngrid) = & (pplay(1:ngrid,1)/(r*pt(1:ngrid,1)))*cpp & *sqrt(pu(1:ngrid,1)*pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1) & +(log(1.+0.7*wstar(1:ngrid) + 2.3*wstar(1:ngrid)**2))**2) & *zcdh(1:ngrid)*(tsurf(1:ngrid)-zh(1:ngrid,1)) endif #else #ifndef MESOINI c ========================================================== c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing c any variable for diagnostic (output with period c "ecritphy", set in "run.def") c ========================================================== c WRITEDIAGFI can ALSO be called from any other subroutines c for any variables !! call write_output("emis","Surface emissivity","", & emis(:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("emis_slope"//str2, & "Surface emissivity","",emis(:,islope)) ENDDO call write_output("zzlay","Midlayer altitude", & "m",zzlay(:,:)) call write_output("zzlev","Interlayer altitude", & "m",zzlev(:,1:nlayer)) call write_output("pphi","Geopotential","m2s-2", & pphi(:,:)) call write_output("phisfi","Surface geopotential", & "m2s-2",phisfi(:)) if (grid_type == regular_lonlat) then call write_output("area","Mesh area","m2", & cell_area_for_lonlat_outputs) else ! unstructured grid (e.g. dynamico) call write_output("area","Mesh area","m2",cell_area) endif call write_output("tsurf","Surface temperature","K", & tsurf(:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("tsurf_slope"//str2, & "Surface temperature","K", & tsurf(:,islope)) ENDDO call write_output("ps","surface pressure","Pa",ps(:)) call write_output("co2ice","co2 ice thickness" & ,"kg.m-2",qsurf(:,igcm_co2,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("co2ice_slope"//str2,"co2 ice thickness" & ,"kg.m-2",qsurf(:,igcm_co2,islope)) ENDDO call write_output("watercap","Perennial water ice thickness" & ,"kg.m-2",watercap(:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("watercap_slope"//str2, & "Perennial water ice thickness" & ,"kg.m-2",watercap(:,islope)) ENDDO call write_output("temp_layer1","temperature in layer 1", & "K",zt(:,1)) call write_output("temp7","temperature in layer 7", & "K",zt(:,7)) call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2", & fluxsurf_lw(:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("fluxsurf_lw_slope"//str2, & "fluxsurf_lw","W.m-2", & fluxsurf_lw(:,islope)) ENDDO call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw", & "W.m-2",fluxsurf_dn_sw_tot(:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("fluxsurf_dn_sw_slope"//str2, & "fluxsurf_dn_sw", & "W.m-2",fluxsurf_dn_sw_tot(:,islope)) ENDDO call write_output("fluxtop_dn_sw","fluxtop_dn_sw", & "W.m-2",fluxtop_dn_sw(:,1) + fluxtop_dn_sw(:,2)) call write_output("fluxtop_lw","fluxtop_lw","W.m-2", & fluxtop_lw(:)) call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2", & fluxtop_up_sw_tot(:)) call write_output("temp","temperature","K",zt(:,:)) call write_output("Sols","Time","sols",zday) call write_output("Ls","Solar longitude","deg", & zls*180./pi) call write_output("u","Zonal wind","m.s-1",zu(:,:)) call write_output("v","Meridional wind","m.s-1",zv(:,:)) call write_output("w","Vertical wind","m.s-1",pw(:,:)) call write_output("rho","density","kg.m-3",rho(:,:)) call write_output("pressure","Pressure","Pa",zplay(:,:)) call write_output("zplev","Interlayer pressure","Pa", & zplev(:,1:nlayer)) call write_output('sw_htrt','sw heat. rate', & 'K/s',zdtsw(:,:)) call write_output('lw_htrt','lw heat. rate', & 'K/s',zdtlw(:,:)) call write_output("local_time","Local time", & 'sol',local_time(:)) if (.not.activice) then CALL write_output('tauTESap', & 'tau abs 825 cm-1', & '',tauTES(:)) else CALL write_output('tauTES', & 'tau abs 825 cm-1', & '',taucloudtes(:)) endif #else !!! this is to ensure correct initialisation of mesoscale model call write_output("tsurf","Surface temperature","K", & tsurf(:,iflat)) call write_output("ps","surface pressure","Pa",ps(:)) call write_output("co2ice","co2 ice thickness","kg.m-2", & qsurf(:,igcm_co2,iflat)) call write_output("temp","temperature","K",zt(:,:)) call write_output("u","Zonal wind","m.s-1",zu(:,:)) call write_output("v","Meridional wind","m.s-1",zv(:,:)) call write_output("emis","Surface emissivity","", & emis(:,iflat)) call write_output("tsoil","Soil temperature", & "K",tsoil(:,:,iflat)) call write_output("inertiedat","Soil inertia", & "K",inertiedat(:,:)) #endif c ---------------------------------------------------------- c Outputs of the CO2 cycle c ---------------------------------------------------------- if (igcm_co2.ne.0) then call write_output("co2","co2 mass mixing ratio", & "kg.kg-1",zq(:,:,igcm_co2)) if (co2clouds) then call write_output('ccnqco2','CCNco2 mmr', & 'kg.kg-1',zq(:,:,igcm_ccnco2_mass)) call write_output('ccnNco2','CCNco2 number', & 'part.kg-1',zq(:,:,igcm_ccnco2_number)) call write_output('co2_ice','co2_ice mmr in atm', & 'kg.kg-1',zq(:,:,igcm_co2_ice)) call write_output("mtotco2","total mass atm of co2", & "kg.m-2",mtotco2(:)) call write_output("icetotco2","total mass atm of co2 ice", & "kg.m-2", icetotco2(:)) call write_output("vaptotco2","total mass atm of co2 "// & "vapor","kg.m-2", vaptotco2(:)) if (co2useh2o) then call write_output('ccnqco2_h2o_m_ice', & 'CCNco2_h2o_mass_ice mmr', & 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice)) call write_output('ccnqco2_h2o_m_ccn', & 'CCNco2_h2o_mass_ccn mmr', & 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn)) call write_output('ccnNco2_h2o','CCNco2_h2o number', & 'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number)) end if if (meteo_flux) then call write_output('ccnqco2_meteor','CCNco2_meteor mmr', & 'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass)) call write_output('ccnNco2_meteor','CCNco2_meteor number', & 'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number)) end if end if ! of if (co2clouds) end if ! of if (igcm_co2.ne.0) ! Output He tracer, if there is one if (igcm_he.ne.0) then call write_output("he","helium mass mixing ratio", & "kg/kg",zq(:,:,igcm_he)) vmr = zq(1:ngrid,1:nlayer,igcm_he) & * mmean(1:ngrid,1:nlayer)/mmol(igcm_he) call write_output('vmr_he','helium vol. mixing ratio', & 'mol/mol',vmr(:,:)) end if c ---------------------------------------------------------- c Outputs of the water cycle c ---------------------------------------------------------- if (water) then #ifdef MESOINI !!!! waterice = q01, voir readmeteo.F90 call write_output('q01',noms(igcm_h2o_ice), & 'kg/kg', & zq(:,:,igcm_h2o_ice)) !!!! watervapor = q02, voir readmeteo.F90 call write_output('q02',noms(igcm_h2o_vap), & 'kg/kg', & zq(:,:,igcm_h2o_vap)) !!!! surface waterice qsurf02 (voir readmeteo) call write_output('qsurf02','surface tracer', & 'kg.m-2', & qsurf(:,igcm_h2o_ice,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('qsurf02_slope'//str2, & 'surface tracer','kg.m-2', & qsurf(:,igcm_h2o_ice,islope)) ENDDO #endif call write_output('mtot', & 'total mass of water vapor', & 'kg/m2',mtot(:)) call write_output('icetot', & 'total mass of water ice', & 'kg/m2',icetot(:)) vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice) & * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice) call write_output('vmr_h2oice','h2o ice vmr', & 'mol/mol',vmr(:,:)) vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap) & * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) call write_output('vmr_h2ovap','h2o vap vmr', & 'mol/mol',vmr(:,:)) call write_output('reffice', & 'Mean reff', & 'm',rave(:)) call write_output('h2o_ice','h2o_ice','kg/kg', & zq(:,:,igcm_h2o_ice)) call write_output('h2o_vap','h2o_vap','kg/kg', & zq(:,:,igcm_h2o_vap)) if (hdo) then vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice) & *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice) CALL write_output('vmr_hdoice','hdo ice vmr', & 'mol/mol',vmr(:,:)) vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap) & *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap) CALL write_output('vmr_hdovap','hdo vap vmr', & 'mol/mol',vmr(:,:)) call write_output('hdo_ice','hdo_ice','kg/kg', & zq(:,:,igcm_hdo_ice)) call write_output('hdo_vap','hdo_vap','kg/kg', & zq(:,:,igcm_hdo_vap)) CALL write_output('mtotD', & 'total mass of HDO vapor', & 'kg/m2',mtotD(:)) CALL write_output('icetotD', & 'total mass of HDO ice', & 'kg/m2',icetotD(:)) C Calculation of the D/H ratio do l=1,nlayer do ig=1,ngrid if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/ & zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6) else DoH_vap(ig,l) = 0. endif enddo enddo do l=1,nlayer do ig=1,ngrid if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/ & zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6) else DoH_ice(ig,l) = 0. endif enddo enddo CALL write_output('DoH_vap', & 'D/H ratio in vapor', & ' ',DoH_vap(:,:)) CALL write_output('DoH_ice', & 'D/H ratio in ice', & '',DoH_ice(:,:)) endif !hdo !A. Pottier ! CALL write_output('rmoym', ! & 'alternative reffice', ! & 'm',rave2(:)) call write_output('h2o_saturation', & 'h2o vap saturation ratio','',satu(:,:)) if (scavenging) then CALL write_output("Nccntot", & "condensation nuclei","Nbr/m2", & Nccntot(:)) CALL write_output("Mccntot", & "mass condensation nuclei","kg/m2", & Mccntot(:)) endif call write_output('rice','Water ice particle size', & 'm',rice(:,:)) call write_output('h2o_ice_s', & 'surface h2o_ice', & 'kg.m-2',qsurf(:,igcm_h2o_ice,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('h2o_ice_s_slope'//str2, & 'surface h2o_ice', & 'kg.m-2',qsurf(:,igcm_h2o_ice,islope)) ENDDO if (hdo) then call write_output('hdo_ice_s', & 'surface hdo_ice', & 'kg.m-2',qsurf(:,igcm_hdo_ice,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('hdo_ice_s_slope'//str2, & 'surface hdo_ice', & 'kg.m-2',qsurf(:,igcm_hdo_ice,islope)) ENDDO do ig=1,ngrid if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/ & qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6 else DoH_surf(ig) = 0. endif enddo call write_output('DoH_surf', & 'surface D/H', & '',DoH_surf(:)) endif ! hdo CALL write_output('albedo', & 'albedo', & '',albedo(:,1,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope CALL write_output('albedo_slope'//str2, & 'albedo', & '',albedo(:,1,islope)) ENDDO if (tifeedback) then call write_output("soiltemp", & "Soil temperature","K", & tsoil(:,:,iflat)) call write_output('soilti', & 'Soil Thermal Inertia', & 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('soilti_slope'//str2, & 'Soil Thermal Inertia', & 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope)) ENDDO endif !A. Pottier if (CLFvarying) then !AP14 nebulosity call write_output('totcloudfrac', & 'Total cloud fraction', & ' ',totcloudfrac(:)) end if !clf varying end if !(water) c ---------------------------------------------------------- c Outputs of the dust cycle c ---------------------------------------------------------- call write_output('tau_pref_scenario', & 'Prescribed visible dust optical depth at 610Pa', & 'NU',tau_pref_scenario(:)) call write_output('tau_pref_gcm', & 'Visible dust optical depth at 610Pa in the GCM', & 'NU',tau_pref_gcm(:)) if (reff_driven_IRtoVIS_scenario) then call write_output('IRtoVIScoef', & 'Conversion coeff for dust tau from abs9.3um to ext0.67um', & '/',IRtoVIScoef(:)) endif if (dustbin.ne.0) then #ifndef MESOINI if (doubleq) then call write_output('dqsdust', & 'deposited surface dust mass', & 'kg.m-2.s-1',dqdustsurf(:)) call write_output('dqndust', & 'deposited surface dust number', & 'number.m-2.s-1',dndustsurf(:)) call write_output('reffdust','reffdust', & 'm',rdust(:,:)*ref_r0) call write_output('dustq','Dust mass mr', & 'kg/kg',qdust(:,:)) call write_output('dustN','Dust number', & 'part/kg',ndust(:,:)) select case (trim(dustiropacity)) case ("tes") call write_output('dsodust_TES', & 'density scaled extinction opacity of std dust at 9.3um(TES)', & 'm2.kg-1',dsodust(:,:)) call write_output('dso_TES', & 'density scaled extinction opacity of all dust at 9.3um(TES)', & 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:)) case ("mcs") call write_output('dsodust', & 'density scaled extinction opacity of std dust at 21.6um(MCS)', & 'm2.kg-1',dsodust(:,:)) call write_output('dso', & 'density scaled extinction opacity of all dust at 21.6um(MCS)', & 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:)) end select else ! (doubleq=.false.) do iq=1,dustbin write(str2(1:2),'(i2.2)') iq call write_output('q'//str2,'mix. ratio', & 'kg/kg',zq(:,:,iq)) call write_output('qsurf'//str2,'qsurf', & 'kg.m-2',qsurf(:,iq,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('qsurf_slope'//str2,'qsurf', & 'kg.m-2',qsurf(:,iq,islope)) ENDDO end do endif ! (doubleq) if (rdstorm) then ! writediagfi tendencies stormdust tracers call write_output('reffstormdust','reffstormdust', & 'm',rstormdust(:,:)*ref_r0) call write_output('mstormdtot', & 'total mass of stormdust only', & 'kg.m-2',mstormdtot(:)) call write_output('mdusttot', & 'total mass of dust only', & 'kg.m-2',mdusttot(:)) call write_output('rdsdqsdust', & 'deposited surface stormdust mass', & 'kg.m-2.s-1',rdsdqdustsurf(:)) call write_output('rdsdustq','storm Dust mass mr', & 'kg/kg',rdsqdust(:,:)) call write_output('rdsdustqmodel','storm Dust massmr', & 'kg/kg',pq(:,:,igcm_stormdust_mass)) call write_output('rdsdustN','storm Dust number', & 'part/kg',rdsndust(:,:)) call write_output("stormfract", & "fraction of the mesh, with stormdust","none", & totstormfract(:)) ! call write_output('qsurf', ! & 'stormdust injection', ! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat)) ! do islope=1,nslope ! write(str2(1:2),'(i2.2)') islope ! call write_output('qsurf_slope'//str2, ! & 'stormdust injection', ! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,islope)) ! ENDDO ! call write_output('pdqsurf', ! & 'tendancy stormdust mass at surface', ! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat)) ! do islope=1,nslope ! write(str2(1:2),'(i2.2)') islope ! call write_output('pdqsurf_slope'//str2, ! & 'tendancy stormdust mass at surface', ! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope)) ! ENDDO call write_output('wspeed_stormdust', & 'vertical velocity of stormdust', & 'm/s',wspeed(:,:)) call write_output('zdqsed_dust_mass' & ,'sedimentation tendency of background dust mmr' & ,'kg/kg.s-1', & zdqsed(:,:,igcm_dust_mass)) call write_output('zdqssed_dust_mass' & ,'sedimentation tendency of background dust on surface' & ,'kg.m-2.s-1', & zdqssed(:,igcm_dust_mass)) call write_output('zdqsed_stormdust_mass' & ,'sedimentation tendency of stormdust dust mmr' & ,'kg/kg.s-1', & zdqsed(:,:,igcm_stormdust_mass)) call write_output('zdqsed_dust_number' & ,'sedimentation tendency of background dust number' & ,'nbr/kg.s-1', & zdqsed(:,:,igcm_dust_number)) call write_output('rdust','rdust', & 'm',rdust(:,:)) call write_output('rstormdust','rstormdust', & 'm',rstormdust(:,:)) select case (trim(dustiropacity)) case ("tes") call write_output('dsords_TES', & 'density scaled extinction opacity of stormdust at 9.3um(TES)', & 'm2.kg-1',dsords(:,:)) case ("mcs") call write_output('dsords', & 'density scaled extinction opacity of stormdust at 21.6um(MCS)', & 'm2.kg-1',dsords(:,:)) end select endif ! (rdstorm) if (topflows) then call write_output('refftopdust', & 'Topdust dust effective radius', & 'm',rtopdust(:,:)*ref_r0) call write_output('topdustq','top Dust mass mr', & 'kg/kg',pq(:,:,igcm_topdust_mass)) call write_output('topdustN','top Dust number', & 'part/kg',pq(:,:,igcm_topdust_number)) select case (trim(dustiropacity)) case ("tes") call write_output('dsotop_TES', & 'density scaled extinction opacity of topdust at 9.3um(TES)', & 'm2.kg-1',dsotop(:,:)) case ("mcs") call write_output('dsotop', & 'density scaled extinction opacity of topdust at 21.6um(MCS)', & 'm2.kg-1',dsotop(:,:)) end select endif ! (topflows) if (dustscaling_mode==2) then call write_output("dust_rad_adjust", & "radiative adjustment coefficient for dust", & "",dust_rad_adjust(:)) endif ! if (scavenging) then ! these outputs should be in the scavenging routine ! call write_output('ccnq','CCN mass mr', ! & 'kg/kg',qccn(:,:)) ! call write_output('ccnN','CCN number', ! & 'part/kg',nccn(:,:)) ! call write_output('surfccnq','Surf nuclei mass mr', ! & 'kg.m-2',qsurf(:,igcm_ccn_mass,iflat)) ! do islope=1,nslope ! write(str2(1:2),'(i2.2)') islope ! call write_output('surfccnq_slope'//str2, ! & 'Surf nuclei mass mr', ! & 'kg.m-2',qsurf(:,igcm_ccn_mass,islope)) ! ENDDO ! call write_output('surfccnN','Surf nuclei number', ! & 'kg.m-2',qsurf(:,igcm_ccn_number,iflat)) ! do islope=1,nslope ! write(str2(1:2),'(i2.2)') islope ! call write_output('surfccnN_slope'//str2, ! & 'Surf nuclei number', ! & 'kg.m-2',qsurf(:,igcm_ccn_number,islope)) ! ENDDO ! endif ! (scavenging) #else ! !!! to initialize mesoscale we need scaled variables ! !!! because this must correspond to starting point for tracers ! call write_output('dustq','Dust mass mr', ! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass)) ! call write_output('dustN','Dust number', ! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number)) ! call write_output('ccn','Nuclei mass mr', ! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass)) ! call write_output('ccnN','Nuclei number', ! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number)) if (freedust) then call write_output('dustq','Dust mass mr', & 'kg/kg',qdust) call write_output('dustN','Dust number', & 'part/kg',ndust) call write_output('ccn','CCN mass mr', & 'kg/kg',qccn) call write_output('ccnN','CCN number', & 'part/kg',nccn) else call write_output('dustq','Dust mass mr', & 'kg/kg',pq(:,:,igcm_dust_mass)) call write_output('dustN','Dust number', & 'part/kg',pq(:,:,igcm_dust_number)) call write_output('ccn','Nuclei mass mr', & 'kg/kg',pq(:,:,igcm_ccn_mass)) call write_output('ccnN','Nuclei number', & 'part/kg',pq(:,:,igcm_ccn_number)) endif #endif end if ! (dustbin.ne.0) c ---------------------------------------------------------- c Thermospheric outputs c ---------------------------------------------------------- if(callthermos) then call write_output("quv","UV heating","K/s", $ zdteuv(:,:)) call write_output("cond","Thermal conduction","K/s", $ zdtconduc(:,:)) !H, H2 and D escape fluxes call write_output("PhiH","H escape flux","s-1", $ PhiEscH) call write_output("PhiH2","H2 escape flux","s-1", $ PhiEscH2) call write_output("PhiD","D escape flux","s-1", $ PhiEscD) endif !(callthermos) call write_output("q15um","15 um cooling","K/s", $ zdtnlte(:,:)) call write_output("qnir","NIR heating","K/s", $ zdtnirco2(:,:)) c ---------------------------------------------------------- c ---------------------------------------------------------- c PBL OUTPUS c ---------------------------------------------------------- c ---------------------------------------------------------- c ---------------------------------------------------------- c Outputs of thermals c ---------------------------------------------------------- if (calltherm) then call write_output('lmax_th', & 'index of vertical extension of thermals', & 'grid level',lmax_th_out(:)) call write_output('zmax_th', & 'vertical extension of thermals','m', & zmax_th(:)) call write_output('hfmax_th', & 'maximum heat flux in thermals','K.m/s', & hfmax_th(:)) call write_output('wstar', & 'maximum thermals vertical velocity','m/s', & wstar(:)) end if c ---------------------------------------------------------- c ---------------------------------------------------------- c END OF PBL OUTPUS c ---------------------------------------------------------- c ---------------------------------------------------------- c ---------------------------------------------------------- c Output in netcdf file "diagsoil.nc" for subterranean c variables (output every "ecritphy", as for writediagfi) c ---------------------------------------------------------- ! Write soil temperature call write_output("soiltemp","Soil temperature","K", & tsoil(:,:,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output("soiltemp_slope"//str2, & "Soil temperature","K", & tsoil(:,:,islope)) ENDDO !PREVIOUSLY IN 1D ONLY call write_output("dtrad","rad. heat. rate", & "K.s-1",dtrad(:,:)) if (rdstorm) then call write_output('aerosol_dust','opacity of env. dust','' & ,aerosol(:,:,iaer_dust_doubleq)) call write_output('aerosol_stormdust', & 'opacity of storm dust','' & ,aerosol(:,:,iaer_stormdust_doubleq)) call write_output('dqsdifdustq', &'tendency due to vertical diffusion of background dust on surface' & ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('dqsdifdustq_slope'//str2, &'tendency due to vertical diffusion of background dust on surface' & ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope)) ENDDO call write_output('dqsdifrdsq', & 'tendency due to vertical diffusion of stormdust on surface', & 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat)) do islope=1,nslope write(str2(1:2),'(i2.2)') islope call write_output('dqsdifrdsq_slope'//str2, & 'tendency due to vertical diffusion of stormdust on surface', & 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope)) ENDDO endif !(rdstorm) if(water) then if (.not.scavenging) then call write_output('zdqcloud_ice','cloud ice', & 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)) call write_output('zdqcloud_vap','cloud vap', & 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap)) call write_output('zdqcloud','cloud', & 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice) & +zdqcloud(:,:,igcm_h2o_vap)) IF (hdo) THEN call write_output('zdqcloud_iceD','cloud ice hdo', & 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice)) call write_output('zdqcloud_vapD','cloud vap hdo', & 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap)) ENDIF ! hdo endif !not.scavenging ! Output needed by the PEM DO ig = 1,ngrid ztmp1 =(1/m_co2 - 1/m_noco2) ztmp2=1/m_noco2 pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2) & * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)*1.e-3)*ps(ig) DO islope = 1,nslope ! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005 rhowater_surf_sat(ig,islope) = & exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o) & / tsurf(ig,islope) & * mmol(igcm_h2o_vap)/(mugaz*r) if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then ! we consider to be at saturation above 1.e-4 kg.m-2 rhowater_surf(ig,islope) = rhowater_surf_sat(ig,islope) else ! otherwise, use vapor partial pressure rhowater_surf(ig,islope) = pvap_surf(ig) & / tsurf(ig,islope) & * mmol(igcm_h2o_vap)/(mugaz*r) endif DO isoil = 1,nsoilmx rhowater_soil(ig,isoil,islope) = & exp(beta_clap_h2o/tsoil(ig,isoil,islope)+alpha_clap_h2o) & / tsoil(ig,isoil,islope) & * mmol(igcm_h2o_vap)/(mugaz*r) ENDDO ENDDO ENDDO CALL write_output("Waterdensity_soil", & "rhowater_soil",'kg.m-3', & rhowater_soil(:,:,iflat)) CALL write_output("Waterdensity_surface", & "rhowater_surface",'kg.m-3', & rhowater_surf(:,iflat)) DO islope = 1,nslope write(str2(1:2),'(i2.2)') islope CALL write_output("Waterdensity_soil_slope"//str2, & "rhowater_soil_slope"//str2,'kg.m-3', & rhowater_soil(:,:,islope)) CALL write_output("Waterdensity_surface"//str2, & "rhowater_surface"//str2,'kg.m-3', & rhowater_surf(:,islope)) ENDDO CALL write_output("h2o_layer1","h2o mass mr in the first layer", & 'kg/kg',zq(:,1,igcm_h2o_vap)) CALL write_output("co2_layer1","co2 mass mr in the first layer", & 'kg/kg',zq(:,1,igcm_co2)) ENDIF ! of IF (water) !PREVIOUSLY IN 1D ONLY c ========================================================== c END OF WRITEDIAGFI c ========================================================== #endif ! of ifdef MESOSCALE c ELSE ! if(ngrid.eq.1) c#ifndef MESOSCALE c write(*, c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)') c & zls*180./pi,odpref,tau_pref_scenario c#endif c END IF ! if(ngrid.ne.1) ! test for co2 conservation with co2 microphysics if (igcm_co2_ice.ne.0) then co2totB = 0. ! added by C.M. do ig=1,ngrid do l=1,nlayer co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g* & (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice) & +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep) enddo co2totB = co2totB + qsurf(ig,igcm_co2,iflat) enddo else co2totB = 0. ! added by C.M. do ig=1,ngrid do l=1,nlayer co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g* & (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep) enddo co2totB = co2totB + qsurf(ig,igcm_co2,iflat) enddo endif ! of if (igcm_co2_ice.ne.0) co2conservation = (co2totA-co2totB)/co2totA call write_output( 'co2conservation', & 'Total CO2 mass conservation in physic', & 'kg', co2conservation) ! XIOS outputs #ifdef CPP_XIOS ! Send fields to XIOS: (NB these fields must also be defined as ! in context_lmdz_physics.xml to be correctly used) CALL send_xios_field("controle",tab_cntrl_mod,1) CALL send_xios_field("ap",ap,1) CALL send_xios_field("bp",bp,1) CALL send_xios_field("aps",aps,1) CALL send_xios_field("bps",bps,1) if (lastcall.and.is_omp_master) then write(*,*) "physiq lastcall: call xios_context_finalize" call xios_context_finalize endif #endif if (check_physics_outputs) then ! Check the validity of updated fields at the end of the physics step call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq) endif icount=icount+1 END SUBROUTINE physiq END MODULE physiq_mod