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,alpha_hmons 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, ! soil thermal inertia & tsoil, nsoilmx,!number of subsurface layers & mlayer,layer ! soil mid layer depths use geometry_mod, only: longitude, latitude, cell_area, & 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, co2ice, emis, & capcal, fluxgrd, qsurf, & hmons,summit,base,watercap,watercaptag use comsaison_h, only: dist_sol, declin, zls, & mu0, fract, local_time 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 use dust_param_mod, only: doubleq, lifting, callddevil, & tauscaling, odpref, dustbin, & dustscaling_mode, dust_rad_adjust, & freedust 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 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 use ioipsl_getin_p_mod, only: getin_p 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_sw_tot(ngrid), fluxsurf_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_swclf(ngrid,2) real fluxtop_lwclf(ngrid) real fluxtop_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 c entrainment by slope winds above sb-grid scale topography REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons REAL hmax,hmin REAL hsummit(ngrid) 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======================================================================= 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, & day_ini,time_phys, & tsurf,tsoil,albedo,emis, & q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, & watercap) #else ! MESOSCALE. Supposedly everything is already set in modules. ! So we just check. And we fill day_ini print*,"check: --- in physiq.F" print*,"check: rad,cpp,g,r,rcp,daysec" print*,rad,cpp,g,r,rcp,daysec PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid) PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx) PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx) PRINT*,'check: midlayer,layer ', mlayer(:),layer(:) PRINT*,'check: tracernames ', noms PRINT*,'check: emis ',emis(1),emis(ngrid) PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1) PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq) PRINT*,'check: co2 ',co2ice(1),co2ice(ngrid) !!! 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 albedo(:,1)=albedodat(:) albedo(:,2)=albedo(:,1) watercap(:)=0.0 #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 write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!" tsurf(:)=pt(:,1) 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. 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 ~~~~~~~~~~~~~~~~~~ IF (tracer) THEN CALL initracer(ngrid,nq,qsurf) ENDIF ! end tracer c Initialize albedo and orbital calculation c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL surfini(ngrid,co2ice,qsurf) CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) c initialize soil c ~~~~~~~~~~~~~~~ IF (callsoil) THEN c Thermal inertia feedback: IF (tifeedback) THEN CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil) CALL soil(ngrid,nsoilmx,firstcall,inertiesoil, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ELSE CALL soil(ngrid,nsoilmx,firstcall,inertiedat, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ENDIF ! of IF (tifeedback) ELSE PRINT*, & '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 (tracer.AND.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) ! no need to create a restart file in 1d if (ecritstart.GT.0) then call physdem0("restartfi.nc",longitude,latitude, & nsoilmx,ngrid,nlayer,nq, & ptimestep,pday,0.,cell_area, & albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, & hmons,summit,base) else call physdem0("restartfi.nc",longitude,latitude, & nsoilmx,ngrid,nlayer,nq, & ptimestep,float(day_end),0.,cell_area, & albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, & hmons,summit,base) endif endif c Initialize mountain mesh fraction for the entrainment by slope wind param. c ~~~~~~~~~~~~~~~ if (slpwind) then !! alpha_hmons calculation if (ngrid.gt.1) then call planetwide_maxval(hmons,hmax ) call planetwide_minval(hmons,hmin ) do ig=1,ngrid alpha_hmons(ig)= 0.5*(hmons(ig)-hmin)/(hmax-hmin) enddo else hmin=0. hmax=23162.1 !set here the height of the sub-grid scaled topography do ig=1,ngrid alpha_hmons(ig)= (hmons(ig)-hmin)/(hmax-hmin) !0.1*(hmons(ig)-hmin)/(hmax-hmin) print*,"1D, hmons=",hmons(ig),"alpha=",alpha_hmons(ig) enddo endif ! (ngrid.gt.1) endif ! if (slpwind) #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 #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 co2totA = co2totA + co2ice(ig) 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 co2totA = co2totA + co2ice(ig) 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 clearatm=.true. !! callradite for background dust in the case of slope wind entrainment 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, & emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout, & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & fluxtop_sw,tau_pref_scenario,tau_pref_gcm, & tau,aerosol,dsodust,tauscaling,dust_rad_adjust, & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice, & rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop, & alpha_hmons,nohmons,clearsky,totcloudfrac) ! 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,emis,mu0,zplev,zplay,pt,tsurf,fract, & dist_sol,igout,zdtlwclf,zdtswclf,fluxsurf_lwclf, & fluxsurf_swclf,fluxtop_lwclf,fluxtop_swclf, & tau_pref_scenario,tau_pref_gcm, & tau,aerosol,dsodust,tauscaling,dust_rad_adjust, & taucloudtesclf,rdust, & rice,nuice,riceco2, nuiceco2,co2ice,rstormdust, & rtopdust,totstormfract, & clearatm,dsords,dsotop,alpha_hmons,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 fluxsurf_lw(ig) = ntf_clf*fluxsurf_lwclf(ig) & + tf_clf*fluxsurf_lw(ig) fluxsurf_sw(ig,1) = ntf_clf*fluxsurf_swclf(ig,1) & + tf_clf*fluxsurf_sw(ig,1) fluxsurf_sw(ig,2) = ntf_clf*fluxsurf_swclf(ig,2) & + tf_clf*fluxsurf_sw(ig,2) fluxtop_lw(ig) = ntf_clf*fluxtop_lwclf(ig) & + tf_clf*fluxtop_lw(ig) fluxtop_sw(ig,1) = ntf_clf*fluxtop_swclf(ig,1) & + tf_clf*fluxtop_sw(ig,1) fluxtop_sw(ig,2) = ntf_clf*fluxtop_swclf(ig,2) & + tf_clf*fluxtop_sw(ig,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 print *, 'Slope scheme is on and computing...' DO ig=1,ngrid sl_the = theta_sl(ig) IF (sl_the .ne. 0.) THEN ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2) 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) sl_psi = psi_sl(ig) sl_fl0 = fluxsurf_sw(ig,l) sl_di0 = 0. if (mu0(ig) .gt. 0.) then sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) sl_di0 = sl_di0*1370./dist_sol/dist_sol sl_di0 = sl_di0/ztim1 sl_di0 = fluxsurf_sw(ig,l)*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_sw(ig,l) = 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 ENDIF 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 fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig) $ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1)) $ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) 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 zplanck(ig)=tsurf(ig)*tsurf(ig) zplanck(ig)=emis(ig)* $ stephan*zplanck(ig)*zplanck(ig) fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig) IF(callslope) THEN sky= (1.+cos(pi*theta_sl(ig)/180.))/2. fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) ENDIF 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 writediagfi(ngrid,'qdustrds0','qdust before rds', & 'kg/kg ',3,qdustrds0) call writediagfi(ngrid,'qstormrds0','qstorm before rds', & 'kg/kg ',3,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,igout,totstormfract, & tauscaling,dust_rad_adjust, c input sub-grid scale cloud & clearsky,totcloudfrac, c input sub-grid scale topography & nohmons,alpha_hmons, 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 writediagfi(ngrid,'qdustrds1','qdust after rds', & 'kg/kg ',3,qdustrds1) call writediagfi(ngrid,'qstormrds1','qstorm after rds', & 'kg/kg ',3,qstormrds1) call writediagfi(ngrid,'qdusttotal0','q sum before rds', & 'kg/m2 ',2,qdusttotal0) call writediagfi(ngrid,'qdusttotal1','q sum after rds', & 'kg/m2 ',2,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 (slpwind) THEN if (ngrid.gt.1) then hsummit(:)=summit(:)-phisfi(:)/g else hsummit(:)=14000. endif 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,igout,aerosol, & tauscaling,dust_rad_adjust, & totstormfract,clearatm, & clearsky,totcloudfrac, & nohmons,hsummit, & 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 (slpwind) 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,dustliftday,local_time) 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 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) 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(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,nq,co2ice,zpopsk, $ ptimestep,capcal,lwrite, $ zplay,zplev,zzlay,zzlev,z0, $ pu,pv,zh,pq,tsurf,emis,qsurf, $ 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 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 if (tracer) 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 end if ! of if (tracer) 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(1:ngrid,igcm_dust_number) zdqdif(1:ngrid,1,igcm_dust_mass) = . -zdqsdif(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) = 0. ENDIF ENDIF ELSE DO ig=1,ngrid zdtsurf(ig)=zdtsurf(ig)+ s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) 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, $ tracer,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 if (tracer) then 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 endif 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 if(tracer) then 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 end if 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: ----------------------------------------- if (tracer) then 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, co2ice) DO iq=1, nq DO ig=1,ngrid dqsurf(ig,iq)=dqsurf(ig,iq)+zdqssed_ccn(ig,iq) ENDDO ! (ig) ENDDO ! (iq) 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 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq) 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 dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq) 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 call calchim(ngrid,nlayer,nq, & ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0, $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, $ zdqcloud,zdqscloud,tau(:,1),co2ice, $ 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 dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) ENDDO ENDDO ! of DO iq=1,nq ! add condensation tendency for H2O2 if (igcm_h2o2.ne.0) then DO ig=1,ngrid dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) & +zdqscloud(ig,igcm_h2o2) ENDDO endif END IF ! of IF (photochem) #endif endif ! of if (tracer) #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,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, co2ice ) co2ice = co2ice * 10000. ENDIF IF (callcond) THEN zdtc(:,:) = 0. zdtsurfc(:) = 0. zduc(:,:) = 0. zdvc(:,:) = 0. zdqc(:,:,:) = 0. CALL co2condens(ngrid,nlayer,nq,ptimestep, $ capcal,zplay,zplev,tsurf,pt, $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, $ co2ice,albedo,emis,rdust, $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, $ fluxsurf_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 IF (tracer) THEN 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 ENDIF ! of IF (tracer) #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----------------------------------------------------------------------- IF (tracer) THEN DO iq=1, nq DO ig=1,ngrid qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq) ENDDO ! (ig) ENDDO ! (iq) ENDIF c----------------------------------------------------------------------- c 12. Surface and sub-surface soil temperature c----------------------------------------------------------------------- c c c 12.1 Increment Surface temperature: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 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 (tracer.AND.water.AND.(ngrid.NE.1)) 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 if ((co2ice(ig).eq.0).and. & (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then if ((watercaptag(ig)).and.(cst_cap_albedo)) then albedo(ig,1) = albedo_h2o_cap albedo(ig,2) = albedo_h2o_cap else albedo(ig,1) = albedo_h2o_frost albedo(ig,2) = 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 ! of do ig=1,ngrid ENDIF ! of IF (tracer.AND.water.AND.(ngrid.NE.1)) c c 12.2 Compute soil temperatures and subsurface heat flux: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (callsoil) THEN c Thermal inertia feedback IF (tifeedback) THEN CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil) CALL soil(ngrid,nsoilmx,.false.,inertiesoil, s ptimestep,tsurf,tsoil,capcal,fluxgrd) ELSE CALL soil(ngrid,nsoilmx,.false.,inertiedat, 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 watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig) ENDDO IF (refill_watercap) THEN DO ig=1,ngrid if (watercaptag(ig).and. & (qsurf(ig,igcm_h2o_ice).gt.frost_metam_threshold)) then watercap(ig)=watercap(ig)+qsurf(ig,igcm_h2o_ice) & - frost_metam_threshold qsurf(ig,igcm_h2o_ice) = frost_metam_threshold endif ! (watercaptag(ig).and. ENDDO ENDIF ! (refill_watercap) THEN c----------------------------------------------------------------------- c 13. Write output files c ---------------------- 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) DO ig=1,ngrid fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) ENDDO 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 PRINT*,'Global diagnostics for the physics' PRINT*,'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,zh,z_out,n_out, & T_out,u_out,ustar,tstar,L_mo,vhf,vvv) ! pourquoi ustar recalcule ici? fait dans vdifc. #ifndef MESOSCALE IF (ngrid .eq. 1) THEN dimout=0 ELSE dimout=2 ENDIF DO n=1,n_out write(zstring, '(F8.6)') z_out(n) call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring), & 'potential temperature at z_out','K',dimout,T_out(:,n)) call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring), & 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n)) ENDDO call WRITEDIAGFI(ngrid,'u_star', & 'friction velocity','m/s',dimout,ustar) call WRITEDIAGFI(ngrid,'teta_star', & 'friction potential temperature','K',dimout,tstar) ! call WRITEDIAGFI(ngrid,'L', ! & 'Monin Obukhov length','m',dimout,L_mo) call WRITEDIAGFI(ngrid,'vvv', & 'Vertical velocity variance at zout','m',dimout,vvv) call WRITEDIAGFI(ngrid,'vhf', & 'Vertical heat flux at zout','m',dimout,vhf) #else T_out1(:)=T_out(:,1) u_out1(:)=u_out(:,1) #endif ENDIF c ---------------------------------------------------------- c ---------------------------------------------------------- c END OF SURFACE LAYER INTERPOLATIONS c ---------------------------------------------------------- c ---------------------------------------------------------- IF (ngrid.NE.1) THEN #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 IF( ((ecritstart.GT.0) .and. . (MOD(icount*iphysiq,ecritstart).EQ.0)) . .or. lastcall ) 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 write(*,'(A,I7,A,F12.5)') . 'PHYSIQ: Ecriture du fichier restartfi ; icount=', . icount,' date=',ztime_fin call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, . ptimestep,ztime_fin, . tsurf,tsoil,co2ice,albedo,emis, . q2,qsurf,tauscaling,totcloudfrac,wstar, . watercap) ENDIF #endif c ------------------------------------------------------------------- c Calculation of diagnostic variables written in both stats and c diagfi files c ------------------------------------------------------------------- if (tracer) then ! 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 dqdustsurf(ig) = & zdqssed(ig,igcm_dust_mass)*tauscaling(ig) dndustsurf(ig) = & zdqssed(ig,igcm_dust_number)*tauscaling(ig) 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 dqdustsurf(ig) = dqdustsurf(ig) + & zdqssed(ig,igcm_ccn_mass)*tauscaling(ig) dndustsurf(ig) = dndustsurf(ig) + & zdqssed(ig,igcm_ccn_number)*tauscaling(ig) 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 endif ! of if (tracer) #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) call wstats(ngrid,"co2ice","CO2 ice cover", & "kg.m-2",2,co2ice) call wstats(ngrid,"watercap","H2O ice cover", & "kg.m-2",2,watercap) 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) call wstats(ngrid,"fluxsurf_sw", & "Solar radiative flux to surface","W.m-2",2, & fluxsurf_sw_tot) call wstats(ngrid,"fluxtop_lw", & "Thermal IR radiative flux to space","W.m-2",2, & fluxtop_lw) call wstats(ngrid,"fluxtop_sw", & "Solar radiative flux to space","W.m-2",2, & fluxtop_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) c call wstats(ngrid,"ssurf","Surface stress","N.m-2", c & 2,zstress) c call wstats(ngrid,"sw_htrt","sw heat.rate", c & "W.m-2",3,zdtsw) c call wstats(ngrid,"lw_htrt","lw heat.rate", c & "W.m-2",3,zdtlw) 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 (tracer) then 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)) call wstats(ngrid,'albedo', & 'albedo', & '',2,albedo(1,1)) 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 c call wstats(ngrid,'qsurf','qsurf', c & 'kg.m-2',2,qsurf(1,igcm_dust_mass)) c call wstats(ngrid,'Nsurf','N particles', c & 'N.m-2',2,qsurf(1,igcm_dust_number)) c call wstats(ngrid,'dqsdev','ddevil lift', c & 'kg.m-2.s-1',2,zdqsdev(1,1)) c call wstats(ngrid,'dqssed','sedimentation', c & 'kg.m-2.s-1',2,zdqssed(1,1)) c call wstats(ngrid,'dqsdif','diffusion', c & 'kg.m-2.s-1',2,zdqsdif(1,1)) 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)) 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") 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 writediagfi(ngrid,"vmr_"//trim(noms(iq)), $ "Volume mixing ratio","mol/mol",3,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 writediagfi(ngrid,"num_"//trim(noms(iq)), $ "Number density","cm-3",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 writediagfi(ngrid,"c_"//trim(noms(iq)), $ "column","mol cm-2",2,colden(1,iq)) ! global mass (g) call planetwide_sumval(colden(:,iq)/6.022e23 $ *mmol(iq)*1.e4*cell_area(:),mass(iq)) call writediagfi(ngrid,"mass_"//trim(noms(iq)), $ "global mass","g",0,mass(iq)) end if ! of if (noms(iq) .ne. "dust_mass" ...) end do ! of do iq=1,nq end if ! of if (photochem) end if ! of if (tracer) 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_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 WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, & emis) c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, & tsurf) call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness" & ,"kg.m-2",2,co2ice) call WRITEDIAGFI(ngrid,"watercap","Water ice thickness" & ,"kg.m-2",2,watercap) call WRITEDIAGFI(ngrid,"temp_layer1","temperature in layer 1", & "K",2,zt(1,1)) call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7", & "K",2,zt(1,7)) call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, & fluxsurf_lw) call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, & fluxsurf_sw_tot) call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, & fluxtop_lw) call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, & fluxtop_sw_tot) call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) call WRITEDIAGFI(ngrid,"Sols","Time","sols",0,[zday]) call WRITEDIAGFI(ngrid,"Ls","Solar longitude","deg", & 0,[zls*180./pi]) call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) call WRITEDIAGFI(ngrid,"rho","density","kg.m-3",3,rho) c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay) c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, c & zstress) call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate', & 'w.m-2',3,zdtsw) call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate', & 'w.m-2',3,zdtlw) call writediagfi(ngrid,"local_time","Local time", & 'sol',2,local_time) if (.not.activice) then CALL WRITEDIAGFI(ngrid,'tauTESap', & 'tau abs 825 cm-1', & '',2,tauTES) else CALL WRITEDIAGFI(ngrid,'tauTES', & 'tau abs 825 cm-1', & '',2,taucloudtes) endif #else !!! this is to ensure correct initialisation of mesoscale model call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, & tsurf) call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, & co2ice) call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, & emis) call WRITEDIAGFI(ngrid,"tsoil","Soil temperature", & "K",3,tsoil) call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia", & "K",3,inertiedat) #endif c ---------------------------------------------------------- c Outputs of the CO2 cycle c ---------------------------------------------------------- if (tracer.and.(igcm_co2.ne.0)) then call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio", & "kg.kg-1",3,zq(:,:,igcm_co2)) if (co2clouds) then call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mmr', & 'kg.kg-1',3,zq(:,:,igcm_ccnco2_mass)) call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number', & 'part.kg-1',3,zq(:,:,igcm_ccnco2_number)) call WRITEDIAGFI(ngrid,'co2_ice','co2_ice mmr in atm', & 'kg.kg-1', 3, zq(:,:,igcm_co2_ice)) call WRITEDIAGFI(ngrid,"mtotco2","total mass atm of co2", & "kg.m-2",2, mtotco2) call WRITEDIAGFI(ngrid,"icetotco2","total mass atm of co2 ice", & "kg.m-2", 2, icetotco2) call WRITEDIAGFI(ngrid,"vaptotco2","total mass atm of co2 & vapor","kg.m-2", 2, vaptotco2) call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, & emis) if (co2useh2o) then call WRITEDIAGFI(ngrid,'ccnqco2_h2o_m_ice', & 'CCNco2_h2o_mass_ice mmr', & 'kg.kg-1',3,zq(:,:,igcm_ccnco2_h2o_mass_ice)) call WRITEDIAGFI(ngrid,'ccnqco2_h2o_m_ccn', & 'CCNco2_h2o_mass_ccn mmr', & 'kg.kg-1',3,zq(:,:,igcm_ccnco2_h2o_mass_ccn)) call WRITEDIAGFI(ngrid,'ccnNco2_h2o','CCNco2_h2o number', & 'part.kg-1',3,zq(:,:,igcm_ccnco2_h2o_number)) end if if (meteo_flux) then call WRITEDIAGFI(ngrid,'ccnqco2_meteor','CCNco2_meteor mmr', & 'kg.kg-1',3,zq(:,:,igcm_ccnco2_meteor_mass)) call WRITEDIAGFI(ngrid,'ccnNco2_meteor','CCNco2_meteor number', & 'part.kg-1',3,zq(:,:,igcm_ccnco2_meteor_number)) end if end if ! of if (co2clouds) end if ! of if (tracer.and.(igcm_co2.ne.0)) ! Output He tracer, if there is one if (tracer.and.(igcm_he.ne.0)) then call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio", & "kg/kg",3,zq(1,1,igcm_he)) vmr = zq(1:ngrid,1:nlayer,igcm_he) & * mmean(1:ngrid,1:nlayer)/mmol(igcm_he) call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio', & 'mol/mol',3,vmr) end if c ---------------------------------------------------------- c Outputs of the water cycle c ---------------------------------------------------------- if (tracer) then if (water) then #ifdef MESOINI !!!! waterice = q01, voir readmeteo.F90 call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice), & 'kg/kg',3, & zq(1:ngrid,1:nlayer,igcm_h2o_ice)) !!!! watervapor = q02, voir readmeteo.F90 call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap), & 'kg/kg',3, & zq(1:ngrid,1:nlayer,igcm_h2o_vap)) !!!! surface waterice qsurf02 (voir readmeteo) call WRITEDIAGFI(ngrid,'qsurf02','surface tracer', & 'kg.m-2',2, & qsurf(1:ngrid,igcm_h2o_ice)) #endif call WRITEDIAGFI(ngrid,'mtot', & 'total mass of water vapor', & 'kg/m2',2,mtot) call WRITEDIAGFI(ngrid,'icetot', & 'total mass of water ice', & 'kg/m2',2,icetot) vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice) & * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice) call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr', & 'mol/mol',3,vmr) vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap) & * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr', & 'mol/mol',3,vmr) call WRITEDIAGFI(ngrid,'reffice', & 'Mean reff', & 'm',2,rave) call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg', & 3,zq(:,:,igcm_h2o_ice)) call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg', & 3,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 WRITEDIAGFI(ngrid,'vmr_hdoice','hdo ice vmr', & 'mol/mol',3,vmr) vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap) & *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap) CALL WRITEDIAGFI(ngrid,'vmr_hdovap','hdo vap vmr', & 'mol/mol',3,vmr) call WRITEDIAGFI(ngrid,'hdo_ice','hdo_ice','kg/kg', & 3,zq(:,:,igcm_hdo_ice)) call WRITEDIAGFI(ngrid,'hdo_vap','hdo_vap','kg/kg', & 3,zq(:,:,igcm_hdo_vap)) CALL WRITEDIAGFI(ngrid,'mtotD', & 'total mass of HDO vapor', & 'kg/m2',2,mtotD) CALL WRITEDIAGFI(ngrid,'icetotD', & 'total mass of HDO ice', & 'kg/m2',2,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 WRITEDIAGFI(ngrid,'DoH_vap', & 'D/H ratio in vapor', & ' ',3,DoH_vap) CALL WRITEDIAGFI(ngrid,'DoH_ice', & 'D/H ratio in ice', & '',3,DoH_ice) endif !hdo !A. Pottier CALL WRITEDIAGFI(ngrid,'rmoym', & 'alternative reffice', & 'm',2,rave2) call WRITEDIAGFI(ngrid,'saturation', & 'h2o vap saturation ratio','dimless',3,satu) if (scavenging) then CALL WRITEDIAGFI(ngrid,"Nccntot", & "condensation nuclei","Nbr/m2", & 2,Nccntot) CALL WRITEDIAGFI(ngrid,"Mccntot", & "mass condensation nuclei","kg/m2", & 2,Mccntot) endif call WRITEDIAGFI(ngrid,'rice','Ice particle size', & 'm',3,rice) call WRITEDIAGFI(ngrid,'h2o_ice_s', & 'surface h2o_ice', & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) if (hdo) then call WRITEDIAGFI(ngrid,'hdo_ice_s', & 'surface hdo_ice', & 'kg.m-2',2,qsurf(1,igcm_hdo_ice)) do ig=1,ngrid if (qsurf(ig,igcm_h2o_ice).gt.qperemin) then DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/ & qsurf(ig,igcm_h2o_ice) )/155.76e-6 else DoH_surf(ig) = 0. endif enddo call WRITEDIAGFI(ngrid,'DoH_surf', & 'surface D/H', & '',2,DoH_surf) endif ! hdo CALL WRITEDIAGFI(ngrid,'albedo', & 'albedo', & '',2,albedo(1,1)) if (tifeedback) then call WRITEDIAGSOIL(ngrid,"soiltemp", & "Soil temperature","K", & 3,tsoil) call WRITEDIAGSOIL(ngrid,'soilti', & 'Soil Thermal Inertia', & 'J.s-1/2.m-2.K-1',3,inertiesoil) endif !A. Pottier if (CLFvarying) then !AP14 nebulosity call WRITEDIAGFI(ngrid,'totcloudfrac', & 'Total cloud fraction', & ' ',2,totcloudfrac) end if !clf varying end if !(water) if (water.and..not.photochem) then iq = nq c write(str2(1:2),'(i2.2)') iq c call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud', c & 'kg.m-2',2,zdqscloud(1,iq)) c call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim', c & 'kg/kg',3,zdqchim(1,1,iq)) c call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif', c & 'kg/kg',3,zdqdif(1,1,iq)) c call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj', c & 'kg/kg',3,zdqadj(1,1,iq)) c call WRITEDIAGFI(ngrid,'dqc'//str2,'var c', c & 'kg/kg',3,zdqc(1,1,iq)) end if !(water.and..not.photochem) end if !tracer c ---------------------------------------------------------- c Outputs of the dust cycle c ---------------------------------------------------------- call WRITEDIAGFI(ngrid,'tau_pref_scenario', & 'Prescribed visible dust optical depth at 610Pa', & 'NU',2,tau_pref_scenario) call WRITEDIAGFI(ngrid,'tau_pref_gcm', & 'Visible dust optical depth at 610Pa in the GCM', & 'NU',2,tau_pref_gcm) if (tracer.and.(dustbin.ne.0)) then call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1)) #ifndef MESOINI if (doubleq) then c call WRITEDIAGFI(ngrid,'qsurf','qsurf', c & 'kg.m-2',2,qsurf(1,igcm_dust_mass)) c call WRITEDIAGFI(ngrid,'Nsurf','N particles', c & 'N.m-2',2,qsurf(1,igcm_dust_number)) c call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift', c & 'kg.m-2.s-1',2,zdqsdev(1,1)) c call WRITEDIAGFI(ngrid,'dqssed','sedimentation', c & 'kg.m-2.s-1',2,zdqssed(1,1)) c call WRITEDIAGFI(ngrid,'dqsdif','diffusion', c & 'kg.m-2.s-1',2,zdqsdif(1,1)) c call WRITEDIAGFI(ngrid,'sedice','sedimented ice', c & 'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice)) c call WRITEDIAGFI(ngrid,'subice','sublimated ice', c & 'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice)) call WRITEDIAGFI(ngrid,'dqsdust', & 'deposited surface dust mass', & 'kg.m-2.s-1',2,dqdustsurf) call WRITEDIAGFI(ngrid,'dqndust', & 'deposited surface dust number', & 'number.m-2.s-1',2,dndustsurf) call WRITEDIAGFI(ngrid,'reffdust','reffdust', & 'm',3,rdust*ref_r0) call WRITEDIAGFI(ngrid,'dustq','Dust mass mr', & 'kg/kg',3,qdust) call WRITEDIAGFI(ngrid,'dustN','Dust number', & 'part/kg',3,ndust) select case (trim(dustiropacity)) case ("tes") call WRITEDIAGFI(ngrid,'dsodust', & 'density scaled extinction opacity of std dust at 9.3um(TES)', & 'm2.kg-1',3,dsodust) call WRITEDIAGFI(ngrid,'dso', & 'density scaled extinction opacity of all dust at 9.3um(TES)', & 'm2.kg-1',3,dsodust+dsords+dsotop) case ("mcs") call WRITEDIAGFI(ngrid,'dsodust', & 'density scaled extinction opacity of std dust at 21.6um(MCS)', & 'm2.kg-1',3,dsodust) call WRITEDIAGFI(ngrid,'dso', & 'density scaled extinction opacity of all dust at 21.6um(MCS)', & 'm2.kg-1',3,dsodust+dsords+dsotop) end select else ! (doubleq=.false.) do iq=1,dustbin write(str2(1:2),'(i2.2)') iq call WRITEDIAGFI(ngrid,'q'//str2,'mix. ratio', & 'kg/kg',3,zq(1,1,iq)) call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf', & 'kg.m-2',2,qsurf(1,iq)) end do endif ! (doubleq) if (rdstorm) then ! writediagfi tendencies stormdust tracers call WRITEDIAGFI(ngrid,'reffstormdust','reffstormdust', & 'm',3,rstormdust*ref_r0) call WRITEDIAGFI(ngrid,'mstormdtot', & 'total mass of stormdust only', & 'kg.m-2',2,mstormdtot) call WRITEDIAGFI(ngrid,'mdusttot', & 'total mass of dust only', & 'kg.m-2',2,mdusttot) call WRITEDIAGFI(ngrid,'rdsdqsdust', & 'deposited surface stormdust mass', & 'kg.m-2.s-1',2,rdsdqdustsurf) call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr', & 'kg/kg',3,rdsqdust) call WRITEDIAGFI(ngrid,'rdsdustqmodel','storm Dust massmr', & 'kg/kg',3,pq(:,:,igcm_stormdust_mass)) call WRITEDIAGFI(ngrid,'rdsdustN','storm Dust number', & 'part/kg',3,rdsndust) call WRITEDIAGFI(ngrid,"stormfract", & "fraction of the mesh, with stormdust","none", & 2,totstormfract) call WRITEDIAGFI(ngrid,'qsurf', & 'stormdust injection', & 'kg.m-2',2,qsurf(:,igcm_stormdust_mass)) call WRITEDIAGFI(ngrid,'pdqsurf', & 'tendancy stormdust mass at surface', & 'kg.m-2',2,dqsurf(:,igcm_stormdust_mass)) call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust', & 'm/s',3,wspeed(:,1:nlayer)) call WRITEDIAGFI(ngrid,'zdqsed_dustq' & ,'sedimentation q','kg.m-2.s-1',3, & zdqsed(:,:,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'zdqssed_dustq' & ,'sedimentation q','kg.m-2.s-1',2, & zdqssed(:,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'zdqsed_rdsq' & ,'sedimentation q','kg.m-2.s-1',3, & zdqsed(:,:,igcm_stormdust_mass)) call WRITEDIAGFI(ngrid,'rdust','rdust', & 'm',3,rdust) call WRITEDIAGFI(ngrid,'rstormdust','rstormdust', & 'm',3,rstormdust) call WRITEDIAGFI(ngrid,'totaldustq','total dust mass', & 'kg/kg',3,qdusttotal) select case (trim(dustiropacity)) case ("tes") call WRITEDIAGFI(ngrid,'dsords', & 'density scaled extinction opacity of stormdust at 9.3um(TES)', & 'm2.kg-1',3,dsords) case ("mcs") call WRITEDIAGFI(ngrid,'dsords', & 'density scaled extinction opacity of stormdust at 21.6um(MCS)', & 'm2.kg-1',3,dsords) end select endif ! (rdstorm) if (slpwind) then call WRITEDIAGFI(ngrid,'refftopdust','refftopdust', & 'm',3,rtopdust*ref_r0) call WRITEDIAGFI(ngrid,'topdustq','top Dust mass mr', & 'kg/kg',3,pq(:,:,igcm_topdust_mass)) call WRITEDIAGFI(ngrid,'topdustN','top Dust number', & 'part/kg',3,pq(:,:,igcm_topdust_number)) select case (trim(dustiropacity)) case ("tes") call WRITEDIAGFI(ngrid,'dsotop', & 'density scaled extinction opacity of topdust at 9.3um(TES)', & 'm2.kg-1',3,dsotop) case ("mcs") call WRITEDIAGFI(ngrid,'dsotop', & 'density scaled extinction opacity of topdust at 21.6um(MCS)', & 'm2.kg-1',3,dsotop) end select endif ! (slpwind) if (dustscaling_mode==2) then call writediagfi(ngrid,"dust_rad_adjust", & "radiative adjustment coefficient for dust", & "",2,dust_rad_adjust) endif if (scavenging) then call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr', & 'kg/kg',3,qccn) call WRITEDIAGFI(ngrid,'ccnN','CCN number', & 'part/kg',3,nccn) call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr', & 'kg.m-2',2,qsurf(1,igcm_ccn_mass)) call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number', & 'kg.m-2',2,qsurf(1,igcm_ccn_number)) endif ! (scavenging) c if (submicron) then c call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr', c & 'kg/kg',3,pq(1,1,igcm_dust_submicron)) c endif ! (submicron) #else ! !!! to initialize mesoscale we need scaled variables ! !!! because this must correspond to starting point for tracers ! call WRITEDIAGFI(ngrid,'dustq','Dust mass mr', ! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass)) ! call WRITEDIAGFI(ngrid,'dustN','Dust number', ! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number)) ! call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr', ! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass)) ! call WRITEDIAGFI(ngrid,'ccnN','Nuclei number', ! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number)) if (freedust) then call WRITEDIAGFI(ngrid,'dustq','Dust mass mr', & 'kg/kg',3,qdust) call WRITEDIAGFI(ngrid,'dustN','Dust number', & 'part/kg',3,ndust) call WRITEDIAGFI(ngrid,'ccn','CCN mass mr', & 'kg/kg',3,qccn) call WRITEDIAGFI(ngrid,'ccnN','CCN number', & 'part/kg',3,nccn) else call WRITEDIAGFI(ngrid,'dustq','Dust mass mr', & 'kg/kg',3,pq(1,1,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'dustN','Dust number', & 'part/kg',3,pq(1,1,igcm_dust_number)) call WRITEDIAGFI(ngrid,'ccn','Nuclei mass mr', & 'kg/kg',3,pq(1,1,igcm_ccn_mass)) call WRITEDIAGFI(ngrid,'ccnN','Nuclei number', & 'part/kg',3,pq(1,1,igcm_ccn_number)) endif #endif end if ! (tracer.and.(dustbin.ne.0)) c ---------------------------------------------------------- c GW non-oro outputs c ---------------------------------------------------------- if(calllott_nonoro) then call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", $ 3,d_u_hin/ptimestep) call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", $ 3,d_v_hin/ptimestep) endif !(calllott_nonoro) c ---------------------------------------------------------- c Thermospheric outputs c ---------------------------------------------------------- if(callthermos) then call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s", $ 3,zdtnlte) call WRITEDIAGFI(ngrid,"quv","UV heating","K/s", $ 3,zdteuv) call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s", $ 3,zdtconduc) call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s", $ 3,zdtnirco2) !H, H2 and D escape fluxes call WRITEDIAGFI(ngrid,"PhiH","H escape flux","s-1", $ 0,[PhiEscH]) call WRITEDIAGFI(ngrid,"PhiH2","H2 escape flux","s-1", $ 0,[PhiEscH2]) call WRITEDIAGFI(ngrid,"PhiD","D escape flux","s-1", $ 0,[PhiEscD]) ! call wstats(ngrid,"PhiH","H escape flux","s-1", ! $ 0,[PhiEscH]) ! call wstats(ngrid,"PhiH2","H2 escape flux","s-1", ! $ 0,[PhiEscH2]) ! call wstats(ngrid,"PhiD","D escape flux","s-1", ! $ 0,[PhiEscD]) ! call wstats(ngrid,"q15um","15 um cooling","K/s", ! $ 3,zdtnlte) ! call wstats(ngrid,"quv","UV heating","K/s", ! $ 3,zdteuv) ! call wstats(ngrid,"cond","Thermal conduction","K/s", ! $ 3,zdtconduc) ! call wstats(ngrid,"qnir","NIR heating","K/s", ! $ 3,zdtnirco2) endif !(callthermos) call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s", $ 3,zdtnlte) call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s", $ 3,zdtnirco2) c ---------------------------------------------------------- c ---------------------------------------------------------- c PBL OUTPUS c ---------------------------------------------------------- c ---------------------------------------------------------- c ---------------------------------------------------------- c Outputs of thermals c ---------------------------------------------------------- if (calltherm) then ! call WRITEDIAGFI(ngrid,'dtke', ! & 'tendance tke thermiques','m**2/s**2', ! & 3,dtke_th) ! call WRITEDIAGFI(ngrid,'d_u_ajs', ! & 'tendance u thermiques','m/s', ! & 3,pdu_th*ptimestep) ! call WRITEDIAGFI(ngrid,'d_v_ajs', ! & 'tendance v thermiques','m/s', ! & 3,pdv_th*ptimestep) ! if (tracer) then ! if (nq .eq. 2) then ! call WRITEDIAGFI(ngrid,'deltaq_th', ! & 'delta q thermiques','kg/kg', ! & 3,ptimestep*pdq_th(:,:,2)) ! end if ! end if call WRITEDIAGFI(ngrid,'zmax_th', & 'hauteur du thermique','m', & 2,zmax_th) call WRITEDIAGFI(ngrid,'hfmax_th', & 'maximum TH heat flux','K.m/s', & 2,hfmax_th) call WRITEDIAGFI(ngrid,'wstar', & 'maximum TH vertical velocity','m/s', & 2,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 writediagsoil(ngrid,"soiltemp","Soil temperature","K", ! & 3,tsoil) ! Write surface temperature ! call writediagsoil(ngrid,"tsurf","Surface temperature","K", ! & 2,tsurf) c ========================================================== c END OF WRITEDIAGFI c ========================================================== #endif ! of ifdef MESOSCALE ELSE ! if(ngrid.eq.1) #ifndef MESOSCALE write(*, & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)') & zls*180./pi,odpref,tau_pref_scenario c ---------------------------------------------------------------------- c Output in grads file "g1d" (ONLY when using testphys1d) c (output at every X physical timestep) c ---------------------------------------------------------------------- c c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2') c CALL writeg1d(ngrid,1,tsurf,'tsurf','K') c CALL writeg1d(ngrid,1,ps,'ps','Pa') c CALL writeg1d(ngrid,nlayer,zt,'T','K') c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1') c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1') c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') ! THERMALS STUFF 1D if(calltherm) then call WRITEDIAGFI(ngrid,'lmax_th', & 'hauteur du thermique','point', & 0,lmax_th_out) call WRITEDIAGFI(ngrid,'zmax_th', & 'hauteur du thermique','m', & 0,zmax_th) call WRITEDIAGFI(ngrid,'hfmax_th', & 'maximum TH heat flux','K.m/s', & 0,hfmax_th) call WRITEDIAGFI(ngrid,'wstar', & 'maximum TH vertical velocity','m/s', & 0,wstar) end if ! of if (calltherm) call WRITEDIAGFI(ngrid,'w','vertical velocity' & ,'m/s',1,pw) call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, & tsurf) call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu) call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv) call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay) call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev) call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho) call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate", & "K.s-1",1,dtrad) call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate', & 'w.m-2',1,zdtsw) call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate', & 'w.m-2',1,zdtlw) call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness" & ,"kg.m-2",0,co2ice) if (igcm_co2.ne.0) then call co2sat(ngrid*nlayer,zt,zqsatco2) do ig=1,ngrid do l=1,nlayer satuco2(ig,l) = zq(ig,l,igcm_co2)* & (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l) c write(*,*) "In PHYSIQMOD, pt,zt,time ",pt(ig,l) c & ,zt(ig,l),ptime enddo enddo endif call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps) call WRITEDIAGFI(ngrid,'temp','Temperature ', & 'K JA',1,zt) c call WRITEDIAGFI(ngrid,'temp2','Temperature ', c & 'K JA2',1,pt) if(tracer) then c CALL writeg1d(ngrid,1,tau,'tau','SI') do iq=1,nq c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') call WRITEDIAGFI(ngrid,trim(noms(iq)), & trim(noms(iq)),'kg/kg',1,zq(1,1,iq)) end do if (doubleq) then call WRITEDIAGFI(ngrid,'rdust','rdust', & 'm',1,rdust) endif ! doubleq 1D if (rdstorm) then call writediagfi(1,'aerosol_dust','opacity of env. dust','' & ,1,aerosol(:,:,iaer_dust_doubleq)) call writediagfi(1,'aerosol_stormdust', & 'opacity of storm dust','' & ,1,aerosol(:,:,iaer_stormdust_doubleq)) call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion', & 'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion', & 'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass)) call WRITEDIAGFI(ngrid,'mstormdtot', & 'total mass of stormdust only', & 'kg.m-2',0,mstormdtot) call WRITEDIAGFI(ngrid,'mdusttot', & 'total mass of dust only', & 'kg.m-2',0,mdusttot) call WRITEDIAGFI(ngrid,'tau_pref_scenario', & 'Prescribed dust ref opt depth at 610 Pa', & 'NU',0,tau_pref_scenario) call WRITEDIAGFI(ngrid,'tau_pref_gcm', & 'Dust ref opt depth at 610 Pa in the GCM', & 'NU',0,tau_pref_gcm) call WRITEDIAGFI(ngrid,'rdsdqsdust', & 'deposited surface stormdust mass', & 'kg.m-2.s-1',0,rdsdqdustsurf) call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr', & 'kg/kg',1,rdsqdust) call WRITEDIAGFI(ngrid,"stormfract", & "fraction of the mesh,with stormdust", & "none",0,totstormfract) call WRITEDIAGFI(ngrid,'rdsqsurf', & 'stormdust at surface', & 'kg.m-2',0,qsurf(:,igcm_stormdust_mass)) call WRITEDIAGFI(ngrid,'qsurf', & 'dust mass at surface', & 'kg.m-2',0,qsurf(:,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust', & 'm/s',1,wspeed) call WRITEDIAGFI(ngrid,'totaldustq','total dust mass', & 'kg/kg',1,qdusttotal) call WRITEDIAGFI(ngrid,'dsords', & 'density scaled opacity of stormdust', & 'm2.kg-1',1,dsords) call WRITEDIAGFI(ngrid,'zdqsed_dustq' & ,'sedimentation q','kg.m-2.s-1',1, & zdqsed(1,:,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'zdqsed_rdsq' & ,'sedimentation q','kg.m-2.s-1',1, & zdqsed(1,:,igcm_stormdust_mass)) endif !(rdstorm 1D) if (water.AND.tifeedback) then call WRITEDIAGFI(ngrid,"soiltemp", & "Soil temperature","K", & 1,tsoil) call WRITEDIAGFI(ngrid,'soilti', & 'Soil Thermal Inertia', & 'J.s-1/2.m-2.K-1',1,inertiesoil) endif end if cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF (water) THEN if (.not.activice) then tauTES=0 do l=1,nlayer Qabsice = min( & max(0.4e6*rice(1,l)*(1.+nuice_ref)-0.05 ,0.),1.2 & ) opTES(1,l)= 0.75 * Qabsice * & zq(1,l,igcm_h2o_ice) * & (zplev(1,l) - zplev(1,l+1)) / g & / (rho_ice * rice(1,l) * (1.+nuice_ref)) tauTES=tauTES+ opTES(1,l) enddo CALL WRITEDIAGFI(ngrid,'tauTESap', & 'tau abs 825 cm-1', & '',0,tauTES) else CALL WRITEDIAGFI(ngrid,'tauTES', & 'tau abs 825 cm-1', & '',0,taucloudtes) endif mtot = 0 icetot = 0 h2otot = qsurf(1,igcm_h2o_ice) if (hdo) THEN mtotD = 0 icetotD = 0 hdotot = qsurf(1,igcm_hdo_ice) ENDIF !hdo do l=1,nlayer mtot = mtot + zq(1,l,igcm_h2o_vap) & * (zplev(1,l) - zplev(1,l+1)) / g icetot = icetot + zq(1,l,igcm_h2o_ice) & * (zplev(1,l) - zplev(1,l+1)) / g if (hdo) THEN mtotD = mtotD + zq(1,l,igcm_hdo_vap) & * (zplev(1,l) - zplev(1,l+1)) / g icetotD = icetotD + zq(1,l,igcm_hdo_ice) & * (zplev(1,l) - zplev(1,l+1)) / g ENDIF !hdo end do h2otot = h2otot+mtot+icetot IF (hdo) then hdotot = hdotot+mtotD+icetotD ENDIF ! hdo CALL WRITEDIAGFI(ngrid,'h2otot', & 'h2otot', & 'kg/m2',0,h2otot) CALL WRITEDIAGFI(ngrid,'mtot', & 'mtot', & 'kg/m2',0,mtot) CALL WRITEDIAGFI(ngrid,'icetot', & 'icetot', & 'kg/m2',0,icetot) IF (hdo) THEN CALL WRITEDIAGFI(ngrid,'mtotD', & 'mtotD', & 'kg/m2',0,mtotD) CALL WRITEDIAGFI(ngrid,'icetotD', & 'icetotD', & 'kg/m2',0,icetotD) CALL WRITEDIAGFI(ngrid,'hdotot', & 'hdotot', & 'kg/m2',0,hdotot) C Calculation of the D/H ratio do l=1,nlayer if (zq(1,l,igcm_h2o_vap).gt.qperemin) then DoH_vap(1,l) = 0.5*( zq(1,l,igcm_hdo_vap)/ & zq(1,l,igcm_h2o_vap) )/155.76e-6 else DoH_vap(1,l) = 0. endif enddo do l=1,nlayer if (zq(1,l,igcm_h2o_ice).gt.qperemin) then DoH_ice(1,l) = 0.5*( zq(1,l,igcm_hdo_ice)/ & zq(1,l,igcm_h2o_ice) )/155.76e-6 else DoH_ice(1,l) = 0. endif enddo CALL WRITEDIAGFI(ngrid,'DoH_vap', & 'D/H ratio in vapor', & ' ',1,DoH_vap) CALL WRITEDIAGFI(ngrid,'DoH_ice', & 'D/H ratio in ice', & '',1,DoH_ice) ENDIF !Hdo if (scavenging) then rave = 0 do l=1,nlayer cccc Column integrated effective ice radius cccc is weighted by total ice surface area (BETTER) rave = rave + tauscaling(1) * & zq(1,l,igcm_ccn_number) * & (zplev(1,l) - zplev(1,l+1)) / g * & rice(1,l) * rice(1,l)* (1.+nuice_ref) enddo rave=icetot*0.75/max(rave*pi*rho_ice,1.e-30) ! surface weight Nccntot= 0 call watersat(ngrid*nlayer,zt,zplay,zqsat) do l=1,nlayer Nccntot = Nccntot + & zq(1,l,igcm_ccn_number)*tauscaling(1) & *(zplev(1,l) - zplev(1,l+1)) / g satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l) satu(1,l) = (max(satu(1,l),float(1))-1) ! & * zq(1,l,igcm_h2o_vap) * ! & (zplev(1,l) - zplev(1,l+1)) / g enddo call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1, & satu) CALL WRITEDIAGFI(ngrid,'Nccntot', & 'Nccntot', & 'nbr/m2',0,Nccntot) call WRITEDIAGFI(ngrid,'zdqsed_dustq' & ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass)) call WRITEDIAGFI(ngrid,'zdqsed_dustN' &,'sedimentation N','Nbr.m-2.s-1',1, & zdqsed(1,:,igcm_dust_number)) else ! of if (scavenging) cccc Column integrated effective ice radius cccc is weighted by total ice mass (LESS GOOD) rave = 0 do l=1,nlayer rave = rave + zq(1,l,igcm_h2o_ice) & * (zplev(1,l) - zplev(1,l+1)) / g & * rice(1,l) * (1.+nuice_ref) enddo rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight endif ! of if (scavenging) CALL WRITEDIAGFI(ngrid,'reffice', & 'reffice', & 'm',0,rave) !Alternative A. Pottier weighting rave2 = 0. totrave2 = 0. do l=1,nlayer rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l) totrave2 = totrave2 + zq(1,l,igcm_h2o_ice) end do rave2=max(rave2/max(totrave2,1.e-30),1.e-30) CALL WRITEDIAGFI(ngrid,'rmoym', & 'reffice', & 'm',0,rave2) do iq=1,nq call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s', & trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq)) end do call WRITEDIAGFI(ngrid,"watercap","Water ice thickness" & ,"kg.m-2",0,watercap) call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice', & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)) call WRITEDIAGFI(ngrid,'zdqcloud_vap','cloud vap', & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap)) call WRITEDIAGFI(ngrid,'zdqcloud','cloud ice', & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice) & +zdqcloud(1,:,igcm_h2o_vap)) IF (hdo) THEN call WRITEDIAGFI(ngrid,'zdqcloud_iceD','cloud ice hdo', & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_ice)) call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo', & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap)) ENDIF ! hdo call WRITEDIAGFI(ngrid,"rice","ice radius","m",1, & rice) if (CLFvarying) then call WRITEDIAGFI(ngrid,'totcloudfrac', & 'Total cloud fraction', & ' ',0,totcloudfrac) endif !clfvarying ENDIF ! of IF (water) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g do l=2,nlayer-1 tmean=zt(1,l) if(zt(1,l).ne.zt(1,l-1)) & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1)) zlocal(l)= zlocal(l-1) & -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g enddo zlocal(nlayer)= zlocal(nlayer-1)- & log(zplay(1,nlayer)/zplay(1,nlayer-1))* & rnew(1,nlayer)*tmean/g #endif 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 + co2ice(ig) 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 + co2ice(ig) enddo endif ! of if (igcm_co2_ice.ne.0) call WRITEDIAGFI(ngrid,'co2conservation', & 'Total CO2 mass conservation in physic', & '%',0,[(co2totA-co2totB)/co2totA]) ! 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("ls",zls*180./pi) 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) CALL send_xios_field("phisinit",phisfi) CALL send_xios_field("ps",ps) CALL send_xios_field("area",cell_area) ! CALL send_xios_field("ISR",fluxtop_sw_tot) CALL send_xios_field("OLR",fluxtop_lw) CALL send_xios_field("tsurf",tsurf) ! CALL send_xios_field("inertiedat",inertiedat) CALL send_xios_field("tsoil",tsoil) CALL send_xios_field("co2ice",co2ice) ! CALL send_xios_field("temp",zt) CALL send_xios_field("u",zu) CALL send_xios_field("v",zv) ! CALL send_xios_field("rho",rho) ! Orographic Gravity waves tendencies ! if (calllott) then ! CALL send_xios_field("dugw",zdugw/ptimestep) ! CALL send_xios_field("dvgw",zdvgw/ptimestep) ! CALL send_xios_field("dtgw",zdtgw/ptimestep) ! endif !CREATE IF CO2CYCLE ! if (tracer.and.(igcm_co2.ne.0)) then ! CALL send_xios_field("co2",zq(:,:,igcm_co2)) ! endif ! Water cycle ! if (water) then ! CALL send_xios_field("watercap",watercap) !CALL send_xios_field("watercaptag",watercaptag) ! CALL send_xios_field("mtot",mtot) ! CALL send_xios_field("icetot",icetot) ! if (igcm_h2o_vap.ne.0 .and. igcm_h2o_ice.ne.0) then ! CALL send_xios_field("h2o_vap",zq(:,:,igcm_h2o_vap)) ! CALL send_xios_field("h2o_ice",zq(:,:,igcm_h2o_ice)) ! endif ! endif ! if (.not.activice) then ! CALL send_xios_field("tauTESap",tauTES) ! else ! CALL send_xios_field("tauTES",taucloudtes) ! endif ! CALL send_xios_field("h2o_ice_s",qsurf(:,igcm_h2o_ice)) 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