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 write_field_phy, only: Writefield_phy !! use ioipsl_getin_p_mod, only: getin_p use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir use generic_cloud_common_h, only : epsi_generic, Psat_generic use gases_h, only: gnom, gfrac use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV use suaer_corrk_mod, only: suaer_corrk use radii_mod, only: su_aer_radii,haze_reffrad_fix use aerosol_mod, only: iaero_haze, i_haze, haze_prof use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, & dryness use comdiurn_h, only: coslat, sinlat, coslon, sinlon use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa use geometry_mod, only: latitude, longitude, cell_area USE comgeomfi_h, only: totarea, totarea_planet USE tracer_h, only: noms, mmol, radius, rho_q, qext, & igcm_n2,igcm_ch4_gas,igcm_ch4_ice,igcm_haze,& igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,lw_ch4,lw_co,& alpha_lift, alpha_devil, qextrhor, & nesp, is_chim, is_condensable,constants_epsi_generic use time_phylmdz_mod, only: ecritphy, iphysiq, nday use phyetat0_mod, only: phyetat0 use wstats_mod, only: callstats, wstats, mkstats use phyredem, only: physdem0, physdem1 use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval use mod_phys_lmdz_para, only : is_master use planete_mod, only: apoastr, periastr, year_day, peri_day, & obliquit, nres, z0, adjust, tpal use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp use time_phylmdz_mod, only: daysec use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, & callrad, callsoil, nosurf, & callconduct,callmolvis,callmoldiff, & aerohaze, corrk, diagdtau, & diurnal, enertest, fat1au, & icetstep, intheat, iradia, kastprof, & lwrite, mass_redistrib, meanOLR, & fast,fasthaze,haze,metcloud,monoxcloud, & n2cond,nearn2cond,noseason_day,conservn2, & convergeps,kbo,triton,paleo,paleoyears,glaflow, & carbox, methane,condmetsurf,condcosurf, & oldplutovdifc,oldplutocorrk,oldplutosedim, & aerohaze,haze_proffix,haze_radproffix, & source_haze, tsurfmax, albmin_ch4, & season, sedimentation,generic_condensation, & specOLR, & startphy_file, testradtimes, & tracer, UseTurbDiff, & global1d, szangle use generic_tracer_index_mod, only: generic_tracer_index use check_fields_mod, only: check_physics_fields use conc_mod, only: rnew, cpnew, ini_conc_mod use phys_state_var_mod use callcorrk_mod, only: callcorrk use callcorrk_pluto_mod, only: callcorrk_pluto use vdifc_mod, only: vdifc use vdifc_pluto_mod, only: vdifc_pluto use turbdiff_mod, only: turbdiff use turb_mod, only : q2,sensibFlux,turb_resolved use mass_redistribution_mod, only: mass_redistribution use condensation_generic_mod, only: condensation_generic use datafile_mod, only: datadir #ifndef MESOSCALE use vertical_layers_mod, only: presnivs, pseudoalt use mod_phys_lmdz_omp_data, ONLY: is_omp_master #else use comm_wrf, only : comm_HR_SW, comm_HR_LW, & comm_CLOUDFRAC,comm_TOTCLOUDFRAC, comm_RH, & comm_HR_DYN, & comm_DQICE,comm_DQVAP,comm_ALBEQ, & comm_FLUXTOP_DN,comm_FLUXABS_SW, & comm_FLUXTOP_LW,comm_FLUXSURF_SW, & comm_FLUXSURF_LW,comm_FLXGRD, & comm_DTRAIN,comm_DTLSC, & comm_LATENT_HF #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 implicit none !================================================================== ! ! Purpose ! ------- ! Central subroutine for all the physics parameterisations in the ! universal model. Originally adapted from the Mars LMDZ model. ! ! The model can be run without or with tracer transport ! depending on the value of "tracer" in file "callphys.def". ! ! ! It includes: ! ! I. Initialization : ! I.1 Firstcall initializations. ! I.2 Initialization for every call to physiq. ! ! II. Compute radiative transfer tendencies (longwave and shortwave) : ! II.a Option 1 : Call correlated-k radiative transfer scheme. ! II.b Option 2 : Call Newtonian cooling scheme. ! II.! Option 3 : Atmosphere has no radiative effect. ! ! III. Vertical diffusion (turbulent mixing) : ! ! IV. Convection : ! IV.a Thermal plume model ! IV.b Dry convective adjusment ! ! V. Condensation and sublimation of gases (currently just N2). ! ! VI. Tracers ! VI.1. Water and water ice. ! VI.2 Photochemistry ! VI.3. Aerosols and particles. ! VI.4. Updates (pressure variations, surface budget). ! VI.5. Slab Ocean. ! VI.6. Surface Tracer Update. ! ! VII. Surface and sub-surface soil temperature. ! ! VIII. Perform diagnostics and write output files. ! ! ! arguments ! --------- ! ! INPUT ! ----- ! ! ngrid Size of the horizontal grid. ! nlayer Number of vertical layers. ! nq Number of advected fields. ! ! firstcall True at the first call. ! lastcall True at the last call. ! ! pday Number of days counted from the North. Spring equinoxe. ! ptime Universal time (00 WHEN DOWNWARDS !!) real omega(ngrid,nlayer) ! omega velocity (Pa/s, >0 when downward) integer i,l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_gas, igcm_generic_ice logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer real zls ! Solar longitude (radians). real zlss ! Sub solar point longitude (radians). real zday ! Date (time since Ls=0, calculated in sols). REAL,save :: saveday ! saved date REAL,save :: savedeclin ! saved declin real zzlay(ngrid,nlayer) ! Altitude at the middle of the atmospheric layers. real zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries. ! VARIABLES for the thermal plume model (AF24: deleted) ! TENDENCIES due to various processes : ! For Surface Temperature : (K/s) real zdtsurf(ngrid) ! Cumulated tendencies. real zdtsurfmr(ngrid) ! Mass_redistribution routine. real zdtsurfc(ngrid) ! Condense_n2 routine. real zdtsdif(ngrid) ! Turbdiff/vdifc routines. ! real zdtsurf_hyd(ngrid) ! Hydrol routine. ! For Atmospheric Temperatures : (K/s) real dtlscale(ngrid,nlayer) ! Largescale routine. real dt_generic_condensation(ngrid,nlayer) ! condensation_generic routine. real zdtc(ngrid,nlayer) ! Condense_n2 routine. real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines. real zdtmr(ngrid,nlayer) ! Mass_redistribution routine. real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine. real zdtchim(ngrid,nlayer) ! Calchim routine. ! For Surface Tracers : (kg/m2/s) real dqsurf(ngrid,nq) ! Cumulated tendencies. !real zdqsurfc(ngrid) ! Condense_n2 routine. REAL zdqsc(ngrid,nq) ! Condense_n2 routine. real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines. real zdqssed(ngrid,nq) ! Callsedim routine. real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine. ! For Tracers : (kg/kg_of_air/s) real zdqc(ngrid,nlayer,nq) ! Condense_n2 routine. real zdqadj(ngrid,nlayer,nq) ! Convadj routine. real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines. real zdqevap(ngrid,nlayer) ! Turbdiff routine. real zdqsed(ngrid,nlayer,nq) ! Callsedim routine. real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine. real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine. real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine. REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine !$OMP THREADPRIVATE(zdqchim,zdqschim) !! PLUTO variables REAL zdqch4cloud(ngrid,nlayer,nq) REAL zdqsch4cloud(ngrid,nq) REAL zdtch4cloud(ngrid,nlayer) REAL zdqcocloud(ngrid,nlayer,nq) REAL zdqscocloud(ngrid,nq) REAL zdtcocloud(ngrid,nlayer) REAL rice_ch4(ngrid,nlayer) ! Methane ice geometric mean radius (m) REAL rice_co(ngrid,nlayer) ! CO ice geometric mean radius (m) REAL zdqsch4fast(ngrid) ! used only for fast model nogcm REAL zdqch4fast(ngrid) ! used only for fast model nogcm REAL zdqscofast(ngrid) ! used only for fast model nogcm REAL zdqcofast(ngrid) ! used only for fast model nogcm REAL zdqflow(ngrid,nq) REAL zdtconduc(ngrid,nlayer) ! (K/s) REAL zdumolvis(ngrid,nlayer) REAL zdvmolvis(ngrid,nlayer) real zdqmoldiff(ngrid,nlayer,nq) ! Haze relatated tendancies REAL zdqhaze(ngrid,nlayer,nq) REAL zdqprodhaze(ngrid,nq) REAL zdqsprodhaze(ngrid) REAL zdqhaze_col(ngrid) REAL zdqphot_prec(ngrid,nlayer) REAL zdqphot_ch4(ngrid,nlayer) REAL zdqconv_prec(ngrid,nlayer) REAL zdq_source(ngrid,nlayer,nq) ! Fast Haze relatated tendancies REAL fluxbot(ngrid) REAL gradflux(ngrid) REAL fluxlym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface REAL fluxlym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface REAL flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1 REAL flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 REAL zfluxuv ! Lyman alpha flux at 1AU REAL array_zero1(ngrid) REAL array_zero2(ngrid,nlayer) ! For Winds : (m/s/s) real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine. real zdumr(ngrid,nlayer), zdvmr(ngrid,nlayer) ! Mass_redistribution routine. real zdvdif(ngrid,nlayer), zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines. real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines. real zdhadj(ngrid,nlayer) ! Convadj routine. REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine. ! For Pressure and Mass : real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s). real zdmassmr_col(ngrid) ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s). real zdpsrfmr(ngrid) ! Pressure tendency for mass_redistribution routine (Pa/s). ! Local variables for LOCAL CALCULATIONS: ! --------------------------------------- real zflubid(ngrid) real zplanck(ngrid),zpopsk(ngrid,nlayer) REAL zdum1(ngrid,nlayer) REAL zdum2(ngrid,nlayer) real ztim1,ztim2,ztim3, z1,z2 real ztime_fin real zdh(ngrid,nlayer) real gmplanet real taux(ngrid),tauy(ngrid) ! local variables for DIAGNOSTICS : (diagfi & stat) ! ------------------------------------------------- real ps(ngrid) ! Surface Pressure. real zt(ngrid,nlayer) ! Atmospheric Temperature. real zu(ngrid,nlayer),zv(ngrid,nlayer) ! Zonal and Meridional Winds. real zq(ngrid,nlayer,nq) ! Atmospheric Tracers. real zdtadj(ngrid,nlayer) ! Convadj Diagnostic. real zdtdyn(ngrid,nlayer) ! Dynamical Heating (K/s). real zdudyn(ngrid,nlayer) ! Dynamical Zonal Wind tendency (m.s-2). real reff(ngrid,nlayer) ! Effective dust radius (used if doubleq=T). real vmr(ngrid,nlayer) ! volume mixing ratio real time_phys real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic. real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2). ! Pluto adding variables real vmr_ch4(ngrid) ! vmr ch4 real vmr_co(ngrid) ! vmr co real rho(ngrid,nlayer) ! density real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3 real zrho_co(ngrid,nlayer) ! density CO kg.m-3 real zrho_haze(ngrid,nlayer) ! density haze kg.m-3 real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1 real zq1temp_ch4(ngrid) ! real qsat_ch4(ngrid) ! real qsat_ch4_l1(ngrid) ! ! CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix real sensiblehf1(ngrid) ! sensible heat flux real sensiblehf2(ngrid) ! sensible heat flux ! included by RW for H2O Manabe scheme real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj). real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale). ! to test energy conservation (RW+JL) real mass(ngrid,nlayer),massarea(ngrid,nlayer) real dEtot, dEtots, AtmToSurf_TurbFlux real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW !$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW) !JL12 conservation test for mean flow kinetic energy has been disabled temporarily real dtmoist_max,dtmoist_min real dItot, dItot_tmp, dVtot, dVtot_tmp real dWtot, dWtot_tmp, dWtots, dWtots_tmp real psat_tmp ! AF24: to remove? real qsat_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) specific concentration at saturation (kg/kg_of_air). real RH_generic(ngrid,nlayer,nq) ! generic condensable tracers (GCS) Relative humidity. real rneb_generic(ngrid,nlayer,nq) ! GCS cloud fraction (generic condensation). real psat_tmp_generic real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic !$OMP THREADPRIVATE(metallicity) real reffrad_generic_zeros_for_wrf(ngrid,nlayer) ! !!! this is temporary, it is only a list of zeros, it will be replaced when a generic aerosol will be implemented ! For Clear Sky Case. (AF24: deleted) real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW real,save,allocatable :: reffcol(:,:) !$OMP THREADPRIVATE(reffcol) ! Non-oro GW tendencies REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer) REAL d_t_hin(ngrid,nlayer) ! Diagnostics 2D of gw_nonoro REAL zustrhi(ngrid), zvstrhi(ngrid) real :: tsurf2(ngrid) !! real :: flux_o(ngrid),flux_g(ngrid) real :: flux_g(ngrid) real :: flux_sens_lat(ngrid) real :: qsurfint(ngrid,nq) #ifdef MESOSCALE REAL :: lsf_dt(nlayer) REAL :: lsf_dq(nlayer) #endif ! local variables for skin depth check real :: inertia_min,inertia_max real :: diurnal_skin ! diurnal skin depth (m) real :: annual_skin ! anuual skin depth (m) ! 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. !$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs) ! Misc character*2 :: str2 character(len=10) :: tmp1 character(len=10) :: tmp2 !================================================================================================== ! ----------------- ! I. INITIALISATION ! ----------------- ! -------------------------------- ! I.1 First Call Initialisation. ! -------------------------------- if (firstcall) then call getin_p("check_physics_inputs", check_physics_inputs) call getin_p("check_physics_outputs", check_physics_outputs) ! Allocate saved arrays (except for 1D model, where this has already ! been done) #ifndef MESOSCALE if (ngrid>1) call phys_state_var_init(nq) #endif ! Variables set to 0 ! ~~~~~~~~~~~~~~~~~~ dtrad(:,:) = 0.0 fluxrad(:) = 0.0 tau_col(:) = 0.0 zdtsw(:,:) = 0.0 zdtlw(:,:) = 0.0 ! Initialize tracer names, indexes and properties. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) if (tracer) then call initracer(ngrid,nq) ! if(photochem) then !AF24: removed endif ! Initialize aerosol indexes. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! call iniaerosol ! allocate related local arrays ! (need be allocated instead of automatic because of "naerkind") allocate(aerosol(ngrid,nlayer,naerkind)) allocate(reffcol(ngrid,naerkind)) #ifdef CPP_XIOS ! Initialize XIOS context write(*,*) "physiq: call wxios_context_init" CALL wxios_context_init #endif ! Read 'startfi.nc' file. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #ifndef MESOSCALE call phyetat0(startphy_file, & ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq, & day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,inertiedat) #else day_ini = pday #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 if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!" tsurf(:)=pt(:,1) if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!" do isoil=1,nsoilmx tsoil(1:ngrid,isoil)=tsurf(1:ngrid) enddo if (is_master) write(*,*) "Physiq: initializing day_ini to pday !" day_ini=pday endif #endif if (pday.ne.day_ini) then write(*,*) "ERROR in physiq.F90:" write(*,*) "bad synchronization between physics and dynamics" write(*,*) "dynamics day: ",pday write(*,*) "physics day: ",day_ini stop endif write (*,*) 'In physiq day_ini =', day_ini ! Initialize albedo calculation. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ albedo(:,:)=0.0 albedo_bareground(:)=0.0 albedo_snow_SPECTV(:)=0.0 albedo_n2_ice_SPECTV(:)=0.0 ptime0=ptime write (*,*) 'In physiq ptime0 =', ptime call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV) ! Initialize orbital calculation. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) savedeclin=0. saveday=pday adjust=0. ! albedo adjustment for convergeps ! Initialize soil. ! ~~~~~~~~~~~~~~~~ if (callsoil) then call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, & ptimestep,tsurf,tsoil,capcal,fluxgrd) else ! else of 'callsoil'. print*,'WARNING! Thermal conduction in the soil turned off' capcal(:)=1.e6 fluxgrd(:)=intheat print*,'Flux from ground = ',intheat,' W m^-2' endif ! end of 'callsoil'. icount=1 ! Initialize surface history variable. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ qsurf_hist(:,:)=qsurf(:,:) !! call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) ! Initialize variable for dynamical heating and zonal wind tendency diagnostic ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ztprevious(:,:)=pt(:,:) zuprevious(:,:)=pu(:,:) ! Set temperature just above condensation temperature (for Early Mars) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(nearn2cond) then write(*,*)' WARNING! Starting at Tcond+1K' do l=1, nlayer do ig=1,ngrid pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4 & -pt(ig,l)) / ptimestep enddo enddo endif if(meanOLR)then call system('rm -f rad_bal.out') ! to record global radiative balance. call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures. call system('rm -f h2o_bal.out') ! to record global hydrological balance. endif !! Initialize variables for water cycle ! AF24: removed ! Set metallicity for GCS ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic ! Set some parameters for the thermal plume model !AF24: removed ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #ifndef MESOSCALE if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d. call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & ptimestep,pday+nday,time_phys,cell_area, & albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) endif !! call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) #endif if (corrk) then ! We initialise the spectral grid here instead of ! at firstcall of callcorrk so we can output XspecIR, XspecVI ! when using Dynamico if (is_master) print*, "physiq_mod: Correlated-k data base folder:",trim(datadir) call getin_p("corrkdir",corrkdir) if (is_master) print*,"corrkdir = ", corrkdir write (tmp1, '(i4)') L_NSPECTI write (tmp2, '(i4)') L_NSPECTV banddir=trim(trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))) banddir=trim(trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))) call setspi !Basic infrared properties. call setspv ! Basic visible properties. call sugas_corrk ! Set up gaseous absorption properties. if (aerohaze) then call suaer_corrk ! Set up aerosol optical properties. endif endif !! call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) ! XIOS outputs #ifdef CPP_XIOS write(*,*) "physiq: call initialize_xios_output" call initialize_xios_output(pday,ptime,ptimestep,daysec, & year_day,presnivs,pseudoalt,WNOI,WNOV) #endif !! call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) write(*,*) "physiq: end of firstcall" endif ! end of 'firstcall' !! call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_gas),1) !! call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) 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 ! call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat) ! call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal) ! ------------------------------------------------------ ! I.2 Initializations done at every physical timestep: ! ------------------------------------------------------ #ifdef CPP_XIOS ! update XIOS time/calendar call update_xios_timestep #endif ! Initialize various variables ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if ( .not.nearn2cond ) then pdt(1:ngrid,1:nlayer) = 0.0 endif zdtsurf(1:ngrid) = 0.0 pdq(1:ngrid,1:nlayer,1:nq) = 0.0 dqsurf(1:ngrid,1:nq)= 0.0 pdu(1:ngrid,1:nlayer) = 0.0 pdv(1:ngrid,1:nlayer) = 0.0 pdpsrf(1:ngrid) = 0.0 zflubid(1:ngrid) = 0.0 flux_sens_lat(1:ngrid) = 0.0 taux(1:ngrid) = 0.0 tauy(1:ngrid) = 0.0 if (conservn2) then write(*,*) 'conservn2 iniloop' call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) endif zday=pday+ptime ! Compute time, in sols (and fraction thereof). ! Compute Stellar Longitude (Ls), and orbital parameters. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (season) then call stellarlong(zday,zls) else call stellarlong(noseason_day,zls) end if ! Get Lyman alpha flux at specific Ls if (haze) then call lymalpha(zls,zfluxuv) print*, 'Haze lyman-alpha zls,zfluxuv=',zls,zfluxuv end if IF (triton) then CALL orbitetriton(zls,zday,dist_star,declin) ELSE call orbite(zls,dist_star,declin,right_ascen) ENDIF if (diurnal) then zlss=-2.*pi*(zday-.5) else if(diurnal .eqv. .false.) then zlss=9999. endif glat(:) = g !AF24: removed oblateness ! Compute geopotential between layers. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer) do l=1,nlayer zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid) enddo zzlev(1:ngrid,1)=0. do l=2,nlayer do ig=1,ngrid z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) enddo enddo !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22 zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1) ! Compute potential temperature ! Note : Potential temperature calculation may not be the same in physiq and dynamic... ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do l=1,nlayer do ig=1,ngrid zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp zh(ig,l)=pt(ig,l)/zpopsk(ig,l) mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig) massarea(ig,l)=mass(ig,l)*cell_area(ig) enddo enddo ! 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 if (.not.fast) then 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)) enddo ! omega in Pa/s do l=1,nlayer-1 omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1)) enddo omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0 do l=1,nlayer omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid) enddo endif if (conservn2) then write(*,*) 'conservn2 thermo' call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) endif ! -------------------------------------------------------- ! 1.3 thermosphere ! -------------------------------------------------------- ! ajout de la conduction depuis la thermosphere IF (callconduct) THEN call conduction (ngrid,nlayer,ptimestep, & pplay,pt,zzlay,zzlev,zdtconduc,tsurf) DO l=1,nlayer DO ig=1,ngrid pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l) ENDDO ENDDO ENDIF ! ajout de la viscosite moleculaire IF (callmolvis) THEN call molvis(ngrid,nlayer,ptimestep, & pplay,pt,zzlay,zzlev, & zdtconduc,pu,tsurf,zdumolvis) call molvis(ngrid,nlayer,ptimestep, & pplay,pt,zzlay,zzlev, & zdtconduc,pv,tsurf,zdvmolvis) DO l=1,nlayer DO ig=1,ngrid ! pdt(ig,l)=pdt(ig,l)+ zdtconduc(ig,l) pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) ENDDO ENDDO ENDIF IF (callmoldiff) THEN call moldiff_red(ngrid,nlayer,nq, & pplay,pplev,pt,pdt,pq,pdq,ptimestep, & zzlay,zdtconduc,zdqmoldiff) DO l=1,nlayer DO ig=1,ngrid DO iq=1, nq pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) ENDDO ENDDO ENDDO ENDIF if (conservn2) then write(*,*) 'conservn2 thermo' call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) endif !--------------------------------- ! II. Compute radiative tendencies !--------------------------------- ! Saving qsurf to compute paleo flux condensation/sublimation DO iq=1, nq DO ig=1,ngrid IF (qsurf(ig,iq).lt.0.) then qsurf(ig,iq)=0. ENDIF qsurf1(ig,iq)=qsurf(ig,iq) ENDDO ENDDO ! Compute local stellar zenith angles ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ fract = 0 if (diurnal) then ztim1=SIN(declin) ztim2=COS(declin)*COS(2.*pi*(zday-.5)) ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) call stelang(ngrid,sinlon,coslon,sinlat,coslat, & ztim1,ztim2,ztim3,mu0,fract) else if(diurnal .eqv. .false.) then call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad) ! WARNING: this function appears not to work in 1D if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight. mu0 = cos(pi*szangle/180.0) fract= 1/(4*mu0) ! AF24: from pluto.old endif endif ! Pluto albedo /IT changes depending on surface ices (only in 3D) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (ngrid.ne.1) then !! Specific to change albedo of N2 so that Psurf !! converges toward 1.4 Pa in "1989" seasons for Triton !! converges toward 1.1 Pa in "2015" seasons for Pluto if (convergeps) then if (triton) then ! 1989 declination if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45. & .and.zday.gt.saveday+1000. & .and.declin.lt.savedeclin) then call globalaverage2d(ngrid,pplev(:,1),globave) if (globave.gt.1.5) then adjust=adjust+0.005 else if (globave.lt.1.3) then adjust=adjust-0.005 endif saveday=zday endif else ! Pluto : 2015 declination current epoch if (declin*180./pi.gt.50.and.declin*180./pi.lt.51. & .and.zday.gt.saveday+10000. & .and.declin.gt.savedeclin) then call globalaverage2d(ngrid,pplev(:,1),globave) if (globave.gt.1.2) then adjust=adjust+0.005 else if (globave.lt.1.) then adjust=adjust-0.005 endif saveday=zday endif endif end if end if ! if ngrid ne 1 call surfprop(ngrid,nq,fract,qsurf,tsurf, & capcal,adjust,dist_star,flusurfold,ptimestep,zls,& albedo,emis,inertiedat) ! do i=2,ngrid ! albedo(i,:) = albedo(1,:) ! enddo if (firstcall.and.callsoil) then ! AF24 Originally in soil.F, but inertiedat is modified by surfprop ! Additional checks: is the vertical discretization sufficient ! to resolve diurnal and annual waves? do ig=1,ngrid ! extreme inertia for this column inertia_min=minval(inertiedat(ig,:)) inertia_max=maxval(inertiedat(ig,:)) print*, "inertia min max" , inertia_min, inertia_max ! diurnal and annual skin depth diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi) annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi) if (0.5*diurnal_skinlayer(nsoilmx)) then ! one should have the full soil be at least twice the diurnal skin depth write(*,*) "soil Error: grid point ig=",ig write(*,*) " longitude=",longitude(ig)*(180./pi) write(*,*) " latitude=",latitude(ig)*(180./pi) write(*,*) " total soil layer depth ",layer(nsoilmx) write(*,*) " not large enough for an annual skin depth of ", & annual_skin write(*,*) " change soil layer distribution (comsoil_h.F90)" call abort_physic("soil","change soil layer distribution (comsoil_h.F90)",1) endif enddo ! of do ig=1,ngrid end if ! callsoil if (callrad) then if( mod(icount-1,iradia).eq.0.or.lastcall) then ! Eclipse incoming sunlight !AF24: removed !! call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) !! call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) if (corrk) then ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! II.a Call correlated-k radiative transfer scheme ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(kastprof)then print*,'kastprof should not = true here' call abort endif ! if(water) then !AF24: removed if(generic_condensation) then do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer epsi_generic=constants_epsi_generic(iq) muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_gas)) muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_gas)) endif end do ! do iq=1,nq loop on tracers ! take into account generic condensable specie (GCS) effect on mean molecular weight else muvar(1:ngrid,1:nlayer+1)=mugaz endif ! if(ok_slab_ocean) then !AF24: removed ! standard callcorrk if (oldplutocorrk) then call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf, & albedo(:,1),emis,mu0,pplev,pplay,pt, & zzlay,tsurf,fract,dist_star,aerosol, & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw, & fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, & cloudfrac,totcloudfrac,.false., & firstcall,lastcall) albedo_equivalent(1:ngrid)=albedo(1:ngrid,1) fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ & fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid,1)) else call callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & zzlay,tsurf,fract,dist_star,aerosol,muvar, & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & fluxsurfabs_sw,fluxtop_lw, & fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, & int_dtaui,int_dtauv, & tau_col,cloudfrac,totcloudfrac, & .false.,firstcall,lastcall) ! Radiative flux from the sky absorbed by the surface (W.m-2). GSR=0.0 fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+ & fluxsurfabs_sw(1:ngrid) endif ! oldplutocorrk !GG (feb2021): Option to "artificially" decrease the raditive time scale in !the deep atmosphere press > 0.1 bar. Suggested by MT. !! COEFF_RAD to be "tuned" to facilitate convergence of tendency !coeff_rad=0. ! 0 values, it doesn't accelerate the convergence !coeff_rad=0.5 !coeff_rad=1. !do l=1, nlayer ! do ig=1,ngrid ! if(pplay(ig,l).ge.1.d4) then ! zdtsw(ig,l)=zdtsw(ig,l)*(pplay(ig,l)/1.d4)**coeff_rad ! zdtlw(ig,l)=zdtlw(ig,l)*(pplay(ig,l)/1.d4)**coeff_rad ! endif ! enddo !enddo ! AF24: removed CLFvarying Option !if(noradsurf)then ! no lower surface; SW flux just disappears ! GSR = SUM(fluxsurf_sw(1:ngrid)*cell_area(1:ngrid))/totarea ! fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid) ! print*,'SW lost in deep atmosphere = ',GSR,' W m^-2' !endif ! Net atmospheric radiative heating rate (K.s-1) dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer) ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! if (firstcall .and. albedo_spectral_mode) then call spectral_albedo_calc(albedo_snow_SPECTV) endif ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! II.b Call Newtonian cooling scheme !AF24: removed ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! II.! Atmosphere has no radiative effect ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ fluxtop_dn(1:ngrid) = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2 if(ngrid.eq.1)then ! / by 4 globally in 1D case... fluxtop_dn(1) = fract(1)*Fat1AU/dist_star**2/2.0 endif fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid) print*,'------------WARNING---WARNING------------' ! by MT2015. print*,'You are in corrk=false mode, ' print*,'and the surface albedo is taken equal to the first visible spectral value' albedo_equivalent(1:ngrid)=albedo(1:ngrid,1) fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1)) fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid) fluxrad_sky(1:ngrid) = fluxsurfabs_sw(1:ngrid) fluxtop_lw(1:ngrid) = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4 dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating endif ! end of corrk endif ! of if(mod(icount-1,iradia).eq.0) ! Transformation of the radiative tendencies ! ------------------------------------------ zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid) zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid) fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid) pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer) ! Test of energy conservation !---------------------------- if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW) call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW) !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk call planetwide_sumval(fluxsurfabs_sw(:)*cell_area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*cell_area(:)/totarea_planet,dEtotsLW) dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:) dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:) if (is_master) then print*,'---------------------------------------------------------------' print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2' print*,'In corrk LW atmospheric heating =',dEtotLW,' W m-2' print*,'atmospheric net rad heating (SW+LW) =',dEtotLW+dEtotSW,' W m-2' print*,'In corrk SW surface heating =',dEtotsSW,' W m-2' print*,'In corrk LW surface heating =',dEtotsLW,' W m-2' print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2' endif endif ! end of 'enertest' endif ! of if (callrad) !! call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) !! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) if (conservn2) then write(*,*) 'conservn2 radiat' call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) endif ! -------------------------------------------- ! III. Vertical diffusion (turbulent mixing) : ! -------------------------------------------- if (calldifv) then zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid) if (oldplutovdifc) then zdum1(:,:) = 0 zdum2(:,:) = 0 zdh(:,:)=pdt(:,:)/zpopsk(:,:) ! Calling vdif (Martian version WITH N2 condensation) CALL vdifc_pluto(ngrid,nlayer,nq,zpopsk, & ptimestep,capcal,lwrite, & pplay,pplev,zzlay,zzlev,z0, & pu,pv,zh,pq,pt,tsurf,emis,qsurf, & zdum1,zdum2,zdh,pdq,pdt,zflubid, & zdudif,zdvdif,zdhdif,zdtsdif,q2, & zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4) zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only bcond=1./tcond1p4Pa acond=r/lw_n2 !------------------------------------------------------------------ ! test methane conservation ! if(methane)then ! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice, ! & ptimestep,pplev,zdqdif,zdqsdif,'CH4',' vdifc ') ! endif ! methane !------------------------------------------------------------------ ! test CO conservation ! if(carbox)then ! call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice, ! & ptimestep,pplev,zdqdif,zdqsdif,'CO ',' vdifc ') ! endif ! carbox !------------------------------------------------------------------ ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff. else if (UseTurbDiff) then call turbdiff(ngrid,nlayer,nq, & ptimestep,capcal, & pplay,pplev,zzlay,zzlev,z0, & pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf, & pdt,pdq,zflubid, & zdudif,zdvdif,zdtdif,zdtsdif, & sensibFlux,q2,zdqdif,zdqevap,zdqsdif, & taux,tauy) else ! if (oldplutovdifc) .and. (UseTurbDiff) zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) call vdifc(ngrid,nlayer,nq,zpopsk, & ptimestep,capcal,lwrite, & pplay,pplev,zzlay,zzlev,z0, & pu,pv,zh,pq,tsurf,emis,qsurf, & zdh,pdq,zflubid, & zdudif,zdvdif,zdhdif,zdtsdif, & sensibFlux,q2,zdqdif,zdqsdif) zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only zdqevap(1:ngrid,1:nlayer)=0. end if !end of 'UseTurbDiff' zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid) !!! this is always done, except for turbulence-resolving simulations if (.not. turb_resolved) then pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer) pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer) pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer) endif ! if(ok_slab_ocean)then !AF24: removed ! flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid)) ! endif !! call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_gas)) if (tracer) then pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq) dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq) end if ! of if (tracer) !! call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) !! call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) ! test energy conservation !------------------------- if(enertest)then dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:) do ig = 1, ngrid dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground enddo call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot) dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:) call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots) call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux) if (is_master) then if (UseTurbDiff) then print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2' print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2' else print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2' print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2' end if endif ! end of 'is_master' ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere. endif ! end of 'enertest' ! ! Test water conservation. !AF24: removed else ! calldifv ! if(.not.newtonian)then zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid) ! ------------------------------------------------------------------ ! Methane surface sublimation and condensation in fast model (nogcm) ! ------------------------------------------------------------------ if ((methane).and.(fast).and.condmetsurf) THEN call ch4surf(ngrid,nlayer,nq,ptimestep, & tsurf,zdtsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, & zdqch4fast,zdqsch4fast) dqsurf(1:ngrid,igcm_ch4_ice)= dqsurf(1:ngrid,igcm_ch4_ice) + & zdqsch4fast(1:ngrid) pdq(1:ngrid,1,igcm_ch4_gas)= pdq(1:ngrid,1,igcm_ch4_gas) + & zdqch4fast(1:ngrid) zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_ch4*zdqsch4fast(1:ngrid)/capcal(1:ngrid) end if ! ------------------------------------------------------------------ ! CO surface sublimation and condensation in fast model (nogcm) ! ------------------------------------------------------------------ if ((carbox).and.(fast).and.condcosurf) THEN call cosurf(ngrid,nlayer,nq,ptimestep, & tsurf,pplev,pdpsrf,pq,pdq,qsurf,dqsurf, & zdqcofast,zdqscofast) dqsurf(1:ngrid,igcm_co_ice)= dqsurf(1:ngrid,igcm_co_ice) + & zdqscofast(1:ngrid) pdq(1:ngrid,1,igcm_co_gas)= pdq(1:ngrid,1,igcm_co_gas) + & zdqcofast(1:ngrid) zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+lw_co*zdqscofast(1:ngrid)/capcal(1:ngrid) end if endif ! end of 'calldifv' if (conservn2) then write(*,*) 'conservn2 diff' call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ & dqsurf(:,1)*ptimestep) endif !------------------- ! IV. Convection : !------------------- ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ ! IV.a Thermal plume model : AF24: removed ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! IV.b Dry convective adjustment : ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(calladj) then zdh(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) zduadj(1:ngrid,1:nlayer)=0.0 zdvadj(1:ngrid,1:nlayer)=0.0 zdhadj(1:ngrid,1:nlayer)=0.0 call convadj(ngrid,nlayer,nq,ptimestep, & pplay,pplev,zpopsk, & pu,pv,zh,pq, & pdu,pdv,zdh,pdq, & zduadj,zdvadj,zdhadj, & zdqadj) pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer) pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvadj(1:ngrid,1:nlayer) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only if(tracer) then pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) end if ! Test energy conservation if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot) if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2' endif ! ! Test water conservation !AF24: removed endif ! end of 'calladj' !---------------------------------------------- ! Non orographic Gravity Waves: AF24: removed !--------------------------------------------- !----------------------------------------------- ! V. Nitrogen condensation-sublimation : !----------------------------------------------- if (n2cond) then if (.not.tracer) then print*,'We need a N2 ice tracer to condense N2' call abort endif zdqc(:,:,:)=0. zdqsc(:,:)=0. call condense_n2(ngrid,nlayer,nq,ptimestep, & capcal,pplay,pplev,tsurf,pt, & pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, & qsurf(1,igcm_n2),albedo,emis, & zdtc,zdtsurfc,pdpsrf,zduc,zdvc, & zdqc,zdqsc(1,igcm_n2)) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer) pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer) pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer) zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq) dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2) !! call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) !! call writediagfi(ngrid,"condense_n2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) ! test energy conservation if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot) call planetwide_sumval(capcal(:)*zdtsurfc(:)*cell_area(:)/totarea_planet,dEtots) if (is_master) then print*,'In n2cloud atmospheric energy change =',dEtot,' W m-2' print*,'In n2cloud surface energy change =',dEtots,' W m-2' endif endif endif ! end of 'n2cond' if (conservn2) then write(*,*) 'conservn2 n2cond' call testconservmass(ngrid,nlayer,pplev(:,1)+ & pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep) endif !--------------------------------------------- ! VI. Specific parameterizations for tracers !--------------------------------------------- if (tracer) then ! 7a. Methane, CO, and ice ! --------------------------------------- ! Methane ice condensation in the atmosphere ! ---------------------------------------- rice_ch4(:,:)=0 ! initialization needed for callsedim zdqch4cloud(:,:,:)=0. if ((methane).and.(metcloud).and.(.not.fast)) THEN call ch4cloud(ngrid,nlayer,naerkind,ptimestep, & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, & pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud, & nq,rice_ch4) DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,igcm_ch4_gas)=pdq(ig,l,igcm_ch4_gas)+ & zdqch4cloud(ig,l,igcm_ch4_gas) pdq(ig,l,igcm_ch4_ice)=pdq(ig,l,igcm_ch4_ice)+ & zdqch4cloud(ig,l,igcm_ch4_ice) ENDDO ENDDO ! Increment methane ice surface tracer tendency DO ig=1,ngrid dqsurf(ig,igcm_ch4_ice)=dqsurf(ig,igcm_ch4_ice)+ & zdqsch4cloud(ig,igcm_ch4_ice) ENDDO ! update temperature tendancy DO ig=1,ngrid DO l=1,nlayer pdt(ig,l)=pdt(ig,l)+zdtch4cloud(ig,l) ENDDO ENDDO end if ! --------------------------------------- ! CO ice condensation in the atmosphere ! ---------------------------------------- zdqcocloud(:,:,:)=0. IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN call cocloud(ngrid,nlayer,naerkind,ptimestep, & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, & pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud, & nq,rice_co,qsurf(1,igcm_n2),dqsurf(1,igcm_n2)) DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,igcm_co_gas)=pdq(ig,l,igcm_co_gas)+ & zdqcocloud(ig,l,igcm_co_gas) pdq(ig,l,igcm_co_ice)=pdq(ig,l,igcm_co_ice)+ & zdqcocloud(ig,l,igcm_co_ice) ENDDO ENDDO ! Increment CO ice surface tracer tendency DO ig=1,ngrid dqsurf(ig,igcm_co_ice)=dqsurf(ig,igcm_co_ice)+ & zdqscocloud(ig,igcm_co_ice) ENDDO ! update temperature tendancy DO ig=1,ngrid DO l=1,nlayer pdt(ig,l)=pdt(ig,l)+zdtcocloud(ig,l) ENDDO ENDDO ELSE rice_co(:,:)=0 ! initialization needed for callsedim END IF ! of IF (carbox) !------------------------------------------------------------------ ! test methane conservation ! if(methane)then ! call testconserv(ngrid,nlayer,nq,igcm_ch4_gas,igcm_ch4_ice, ! & ptimestep,pplev,zdqch4cloud,zdqsch4cloud,'CH4','ch4clou') ! endif ! methane !------------------------------------------------------------------ ! test CO conservation ! if(carbox)then ! call testconserv(ngrid,nlayer,nq,igcm_co_gas,igcm_co_ice, ! & ptimestep,pplev,zdqcocloud,zdqscocloud,'CO ','cocloud') ! endif ! carbox !------------------------------------------------------------------ ! 7b. Haze particle production ! ------------------- IF (haze) THEN zdqphot_prec(:,:)=0. zdqphot_ch4(:,:)=0. zdqhaze(:,:,:)=0 ! Forcing to a fixed haze profile if haze_proffix if (haze_proffix.and.i_haze.gt.0.) then call haze_prof(ngrid,nlayer,zzlay,pplay,pt, & reffrad,profmmr) zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze)) & /ptimestep else call hazecloud(ngrid,nlayer,nq,ptimestep, & pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze, & zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin) endif DO iq=1, nq ! should be updated DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq) ENDDO ENDDO ENDDO ENDIF IF (fast.and.fasthaze) THEN call prodhaze(ngrid,nlayer,nq,ptimestep,pplev,pq,pdq,dist_star, & mu0,declin,zdqprodhaze,zdqsprodhaze,gradflux,fluxbot, & fluxlym_sol_bot,fluxlym_ipm_bot,flym_sol,flym_ipm) DO ig=1,ngrid pdq(ig,1,igcm_ch4_gas)=pdq(ig,1,igcm_ch4_gas)+ & zdqprodhaze(ig,igcm_ch4_gas) pdq(ig,1,igcm_prec_haze)=pdq(ig,1,igcm_prec_haze)+ & zdqprodhaze(ig,igcm_prec_haze) pdq(ig,1,igcm_haze)=abs(pdq(ig,1,igcm_haze)+ & zdqprodhaze(ig,igcm_haze)) qsurf(ig,igcm_haze)= qsurf(ig,igcm_haze)+ & zdqsprodhaze(ig)*ptimestep ENDDO ENDIF ! ------------------------- ! VI.3. Aerosol particles ! ------------------------- !Generic Condensation if (generic_condensation) then call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay, & pt,pq,pdt,pdq,dt_generic_condensation, & dqvaplscale_generic,dqcldlscale_generic,rneb_generic) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dt_generic_condensation(1:ngrid,1:nlayer) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqvaplscale_generic(1:ngrid,1:nlayer,1:nq) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq)+dqcldlscale_generic(1:ngrid,1:nlayer,1:nq) if(enertest)then do ig=1,ngrid genericconddE(ig) = cpp*SUM(mass(:,:)*dt_generic_condensation(:,:)) enddo call planetwide_sumval(cpp*massarea(:,:)*dt_generic_condensation(:,:)/totarea_planet,dEtot) if (is_master) print*,'In generic condensation atmospheric energy change =',dEtot,' W m-2' end if ! Let's loop on tracers cloudfrac(:,:)=0.0 do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer do l = 1, nlayer do ig=1,ngrid cloudfrac(ig,l)=rneb_generic(ig,l,iq) enddo enddo endif end do ! do iq=1,nq loop on tracers ! endif ! .not. water endif !generic_condensation ! Sedimentation. if (sedimentation) then zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0 zdqssed(1:ngrid,1:nq) = 0.0 if (oldplutosedim)then call callsedim_pluto(ngrid,nlayer,ptimestep, & pplev,zzlev,pt,pdt,rice_ch4,rice_co, & pq,pdq,zdqsed,zdqssed,nq,pphi) else call callsedim(ngrid,nlayer,ptimestep, & pplev,zzlev,pt,pdt,pq,pdq, & zdqsed,zdqssed,nq,pphi) endif ! Whether it falls as rain or snow depends only on the surface temperature pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq) dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq) end if ! end of 'sedimentation' ! --------------- ! VI.4. Updates ! --------------- ! Updating Atmospheric Mass and Tracers budgets. if(mass_redistrib) then zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * 0 ! ( zdqevap(1:ngrid,1:nlayer) & ! ! + zdqrain(1:ngrid,1:nlayer,igcm_h2o_gas) & ! ! + dqmoist(1:ngrid,1:nlayer,igcm_h2o_gas) & ! + dqvaplscale(1:ngrid,1:nlayer) ) do ig = 1, ngrid zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer)) enddo ! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr) ! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col) call writediagfi(ngrid,"mass","mass","kg/m2",3,mass) call mass_redistribution(ngrid,nlayer,nq,ptimestep, & capcal,pplay,pplev,pt,tsurf,pq,qsurf, & pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, & zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq) dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer) pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer) pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer) pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid) zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid) endif ! call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) !! call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_gas)) !! call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) ! ------------------ ! VI.5. Slab Ocean !AF24: removed ! ------------------ ! ----------------------------- ! VI.6. Surface Tracer Update ! ----------------------------- ! AF24: deleted hydrology qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq) ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain. qsurf_hist(:,:) = qsurf(:,:) ! if(ice_update)then ! ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice)) ! endif endif! end of if 'tracer' if (conservn2) then write(*,*) 'conservn2 tracer' call testconservmass(ngrid,nlayer,pplev(:,1)+ & pdpsrf(:)*ptimestep,qsurf(:,1)) endif DO ig=1,ngrid flusurf(ig,igcm_n2)=(qsurf(ig,igcm_n2)- & qsurf1(ig,igcm_n2))/ptimestep flusurfold(ig,igcm_n2)=flusurf(ig,igcm_n2) if (methane) then flusurf(ig,igcm_ch4_ice)=(qsurf(ig,igcm_ch4_ice)- & qsurf1(ig,igcm_ch4_ice))/ptimestep flusurfold(ig,igcm_ch4_ice)=flusurf(ig,igcm_ch4_ice) endif if (carbox) then flusurf(ig,igcm_co_ice)=(qsurf(ig,igcm_co_ice)- & qsurf1(ig,igcm_co_ice))/ptimestep !flusurfold(ig,igcm_co_ice)=flusurf(ig,igcm_co_ice) endif ENDDO !! Special source of haze particle ! ! todo: should be placed in haze and use tendency of n2 instead of flusurf IF (source_haze) THEN write(*,*) "Source haze not supported yet." stop ! call hazesource(ngrid,nlayer,nq,ptimestep, & ! pplev,flusurf,mu0,zdq_source) DO iq=1, nq DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+zdq_source(ig,l,iq) ENDDO ENDDO ENDDO ENDIF !------------------------------------------------ ! VII. Surface and sub-surface soil temperature !------------------------------------------------ ! For diagnostic ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid rho(ig,1) = pplay(ig,1)/(r*pt(ig,1)) sensiblehf1(ig)=rho(ig,1)*cpp*(0.4/log(zzlay(ig,1)/z0))**2* & (pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))**0.5* & (tsurf(ig)-pt(ig,1)) if (calldifv) then sensiblehf2(ig)=zflubid(ig)-capcal(ig)*zdtsdif(ig) end if ENDDO ! VII.1 Increment surface temperature ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature ! Lellouch et al., 2000,2011 IF (tsurfmax) THEN DO ig=1,ngrid if (albedo_equivalent(ig).gt.albmin_ch4.and. & qsurf(ig,igcm_n2).eq.0.) then tsurf(ig)=min(tsurf(ig),54.) endif ENDDO ENDIF ! VII.2 Compute soil temperatures and subsurface heat flux. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (callsoil) then call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, & ptimestep,tsurf,tsoil,capcal,fluxgrd) endif ! ! For output : ! tidat_out(:,:)=0. ! DO l=1,min(nlayermx,nsoilmx) ! tidat_out(:,l)=inertiedat(:,l) ! ENDDO ! Test energy conservation if(enertest)then call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots) if (is_master) print*,'Surface energy change =',dEtots,' W m-2' endif ! VII.3 multiply tendencies of cond/subli for paleo loop only in the ! last Pluto year of the simulation ! Year day must be adapted in the startfi for each object ! Paleo uses year_day to calculate the annual mean tendancies ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (paleo) then if (zday.gt.day_ini+ptime0+nday-year_day) then DO iq=1,nq DO ig=1,ngrid qsurfyear(ig,iq)=qsurfyear(ig,iq)+ & (qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep ENDDO ENDDO endif endif ! VII.4 Glacial flow at each timestep glastep or at lastcall ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (fast.and.glaflow) THEN if ((mod(zday-day_ini-ptime0,glastep)).lt.1. & .or.lastcall) then IF (lastcall) then dstep=mod(zday-day_ini-ptime0,glastep)*daysec else dstep=glastep*daysec endif zdqflow(:,:)=qsurf(:,:) IF (paleo) then call spreadglacier_paleo(ngrid,nq,qsurf, & phisfinew,dstep,tsurf) else call spreadglacier_simple(ngrid,nq,qsurf,dstep) endif zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep if (conservn2) then write(*,*) 'conservn2 glaflow' call testconservmass(ngrid,nlayer,pplev(:,1)+ & pdpsrf(:)*ptimestep,qsurf(:,1)) endif endif ENDIF !--------------------------------------------------- ! VIII. Perform diagnostics and write output files !--------------------------------------------------- ! Note : For output only: the actual model integration is performed in the dynamics. ! Temperature, zonal and meridional winds. zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep !! Recast thermal plume vertical velocity array for outputs !! AF24: removed ! Diagnostic. zdtdyn(1:ngrid,1:nlayer) = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer) zdudyn(1:ngrid,1:nlayer) = (pu(1:ngrid,1:nlayer)-zuprevious(1:ngrid,1:nlayer)) / ptimestep zuprevious(1:ngrid,1:nlayer) = zu(1:ngrid,1:nlayer) if(firstcall)then zdtdyn(1:ngrid,1:nlayer)=0.0 zdudyn(1:ngrid,1:nlayer)=0.0 endif ! Dynamical heating diagnostic. do ig=1,ngrid fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp enddo ! Tracers. zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep ! Surface pressure. ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep call globalaverage2d(ngrid,ps,globave) ! pressure density !pluto specific IF (.not.fast) then ! do ig=1,ngrid do l=1,nlayer zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig) zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig) rho(ig,l) = zplay(ig,l)/(r*zt(ig,l)) enddo zplev(ig,nlayer+1)=pplev(ig,nlayer+1)/pplev(ig,1)*ps(ig) enddo ENDIF ! Surface and soil temperature information call planetwide_sumval(cell_area(:)*tsurf(:)/totarea_planet,Ts1) call planetwide_minval(tsurf(:),Ts2) call planetwide_maxval(tsurf(:),Ts3) if(callsoil)then TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea ! mean temperature at bottom soil layer if (is_master) then print*,' ave[Tsurf] min[Tsurf] max[Tsurf] ave[Tdeep]' print*,Ts1,Ts2,Ts3,TsS end if else if (is_master) then print*,' ave[Tsurf] min[Tsurf] max[Tsurf]' print*,Ts1,Ts2,Ts3 endif end if ! Check the energy balance of the simulation during the run if(corrk)then call planetwide_sumval(cell_area(:)*fluxtop_dn(:)/totarea_planet,ISR) call planetwide_sumval(cell_area(:)*fluxabs_sw(:)/totarea_planet,ASR) call planetwide_sumval(cell_area(:)*fluxtop_lw(:)/totarea_planet,OLR) call planetwide_sumval(cell_area(:)*fluxgrd(:)/totarea_planet,GND) call planetwide_sumval(cell_area(:)*fluxdyn(:)/totarea_planet,DYN) do ig=1,ngrid if(fluxtop_dn(ig).lt.0.0)then print*,'fluxtop_dn has gone crazy' print*,'fluxtop_dn=',fluxtop_dn(ig) print*,'tau_col=',tau_col(ig) print*,'aerosol=',aerosol(ig,:,:) print*,'temp= ',pt(ig,:) print*,'pplay= ',pplay(ig,:) call abort endif end do if(ngrid.eq.1)then DYN=0.0 endif if (is_master) then print*,' ISR ASR OLR GND DYN [W m^-2]' print*, ISR,ASR,OLR,GND,DYN endif if(enertest .and. is_master)then print*,'SW flux/heating difference SW++ - ASR = ',dEtotSW+dEtotsSW-ASR,' W m-2' print*,'LW flux/heating difference LW++ - OLR = ',dEtotLW+dEtotsLW+OLR,' W m-2' print*,'LW energy balance LW++ + ASR = ',dEtotLW+dEtotsLW+ASR,' W m-2' endif if(meanOLR .and. is_master)then if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then ! to record global radiative balance open(92,file="rad_bal.out",form='formatted',position='append') write(92,*) zday,ISR,ASR,OLR close(92) open(93,file="tem_bal.out",form='formatted',position='append') if(callsoil)then write(93,*) zday,Ts1,Ts2,Ts3,TsS else write(93,*) zday,Ts1,Ts2,Ts3 endif close(93) endif endif endif ! end of 'corrk' ! Diagnostic to test radiative-convective timescales in code. if(testradtimes)then open(38,file="tau_phys.out",form='formatted',position='append') ig=1 do l=1,nlayer write(38,*) -1./pdt(ig,l),pt(ig,l),pplay(ig,l) enddo close(38) print*,'As testradtimes enabled,' print*,'exiting physics on first call' call abort endif ! Compute column amounts (kg m-2) if tracers are enabled. if(tracer)then qcol(1:ngrid,1:nq)=0.0 do iq=1,nq do ig=1,ngrid qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer)) enddo enddo endif ! end of 'tracer' ! Calculate RH_generic (Generic Relative Humidity) for diagnostic. if(generic_condensation)then RH_generic(:,:,:)=0.0 do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer do l = 1, nlayer do ig=1,ngrid call Psat_generic(zt(ig,l),pplay(ig,l),metallicity,psat_tmp_generic,qsat_generic(ig,l,iq)) RH_generic(ig,l,iq) = zq(ig,l,igcm_generic_gas) / qsat_generic(ig,l,iq) enddo enddo end if end do ! iq=1,nq endif ! end of 'generic_condensation' if (methane) then IF (fast) then ! zq is the mixing ratio supposingly mixed in all atmosphere DO ig=1,ngrid vmr_ch4(ig)=zq(ig,1,igcm_ch4_gas)* & mmol(igcm_n2)/mmol(igcm_ch4_gas)*100. ENDDO ELSE DO ig=1,ngrid ! compute vmr methane vmr_ch4(ig)=qcol(ig,igcm_ch4_gas)* & g/ps(ig)*mmol(igcm_n2)/mmol(igcm_ch4_gas)*100. ! compute density methane DO l=1,nlayer zrho_ch4(ig,l)=zq(ig,l,igcm_ch4_gas)*rho(ig,l) ENDDO ENDDO ENDIF endif if (carbox) then IF (fast) then DO ig=1,ngrid vmr_co(ig)=zq(ig,1,igcm_co_gas)* & mmol(igcm_n2)/mmol(igcm_co_gas)*100. ENDDO ELSE DO ig=1,ngrid ! compute vmr CO vmr_co(ig)=qcol(ig,igcm_co_gas)* & g/ps(ig)*mmol(igcm_n2)/mmol(igcm_co_gas)*100. ! compute density CO DO l=1,nlayer zrho_co(ig,l)=zq(ig,l,igcm_co_gas)*rho(ig,l) ENDDO ENDDO ENDIF endif zrho_haze(:,:)=0. zdqrho_photprec(:,:)=0. IF (haze.and.aerohaze) then DO ig=1,ngrid DO l=1,nlayer zrho_haze(ig,l)=zq(ig,l,igcm_haze)*rho(ig,l) zdqrho_photprec(ig,l)=zdqphot_prec(ig,l)*rho(ig,l) ENDDO ENDDO ENDIF IF (fasthaze) then DO ig=1,ngrid qcol(ig,igcm_haze)=zq(ig,1,igcm_haze)*pplev(ig,1)/g qcol(ig,igcm_prec_haze)=zq(ig,1,igcm_prec_haze)*pplev(ig,1)/g ENDDO ENDIF ! Info about Ls, declin... IF (fast) THEN if (is_master) write(*,*),'Ls=',zls*180./pi,' dec=',declin*180./pi if (is_master) write(*,*),'zday=',zday,' ps=',globave IF (lastcall) then if (is_master) write(*,*),'lastcall' ENDIF ELSE if (is_master) write(*,*),'Ls=',zls*180./pi,'decli=',declin*180./pi,'zday=',zday ENDIF lecttsoil=0 ! default value for lecttsoil call getin_p("lecttsoil",lecttsoil) IF (lastcall.and.(ngrid.EQ.1).and.(lecttsoil.eq.1)) THEN ! save tsoil temperature profile for 1D profile OPEN(13,file='proftsoil.out',form='formatted') DO i=1,nsoilmx write(13,*) tsoil(1,i) ENDDO CLOSE(13) ENDIF if (is_master) print*,'--> Ls =',zls*180./pi if(lastcall) then ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) !! Update surface ice distribution to iterate to steady state if requested !! AF24: removed ! endif if (paleo) then ! time range for tendencies of ice flux qsurfyear zdt_tot=year_day ! Last year of simulation masslost(:)=0. massacc(:)=0. DO ig=1,ngrid ! update new reservoir of ice on the surface DO iq=1,nq ! kg/m2 to be sublimed or condensed during paleoyears qsurfyear(ig,iq)=qsurfyear(ig,iq)* & paleoyears*365.25/(zdt_tot*daysec/86400.) ! special case if we sublime the entire reservoir !! AF: TODO : fix following lines (real_area), using line below: ! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost) ! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN ! masslost(iq)=masslost(iq)+real_area(ig)* & ! (-qsurfyear(ig,iq)-qsurf(ig,iq)) ! qsurfyear(ig,iq)=-qsurf(ig,iq) ! ENDIF ! IF (qsurfyear(ig,iq).gt.0.) THEN ! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq) ! ENDIF ENDDO ENDDO DO ig=1,ngrid DO iq=1,nq qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq) IF (qsurfyear(ig,iq).gt.0.) THEN qsurfpal(ig,iq)=qsurfpal(ig,iq)- & qsurfyear(ig,iq)*masslost(iq)/massacc(iq) ENDIF ENDDO ENDDO ! Finally ensure conservation of qsurf DO iq=1,nq call globalaverage2d(ngrid,qsurf(:,iq),globaveice(iq)) call globalaverage2d(ngrid,qsurfpal(:,iq), & globavenewice(iq)) IF (globavenewice(iq).gt.0.) THEN qsurfpal(:,iq)=qsurfpal(:,iq)* & globaveice(iq)/globavenewice(iq) ENDIF ENDDO ! update new geopotential depending on the ice reservoir phisfipal(:)=phisfinew(:)+qsurfpal(:,igcm_n2)*g/1000. !phisfipal(ig)=phisfi(ig) if (kbo.or.triton) then ! case of Triton : we do not change the orbital parameters pdaypal=pday ! no increment of pdaypal to keep same evolution of the subsolar point eccpal=1.-periastr/((periastr+apoastr)/2.) !no change of ecc peri_daypal=peri_day ! no change oblipal=obliquit ! no change tpalnew=tpal adjustnew=adjust else ! Pluto ! update new pday and tpal (Myr) to be set in startfi controle pdaypal=int(day_ini+paleoyears*365.25/6.3872) tpalnew=tpal+paleoyears*1.e-6 ! Myrs ! update new N2 ice adjustment (not tested yet on Pluto) adjustnew=adjust ! update milankovitch parameters : obliquity,Lsp,ecc call calcmilank(tpalnew,oblipal,peri_daypal,eccpal) !peri_daypal=peri_day !eccpal=0.009 endif if (is_master) write(*,*) "Paleo peri=",peri_daypal," tpal=",tpalnew if (is_master) write(*,*) "Paleo eccpal=",eccpal," tpal=",tpalnew !---------------------------------------------------------------------- ! Writing NetCDF file "RESTARTFI" at the end of the run !---------------------------------------------------------------------- ! Note: 'restartfi' is stored just before dynamics are stored ! in 'restart'. Between now and the writing of 'restart', ! there will have been the itau=itau+1 instruction and ! a reset of 'time' (lastacll = .true. when itau+1= itaufin) ! thus we store for time=time+dtvr #ifndef MESOSCALE ! create restartfi if (ngrid.ne.1) then print*, "physdem1pal not yet implemented" stop !TODO: import this routine from pluto.old ! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, & ! ptimestep,pdaypal, & ! ztime_fin,tsurf,tsoil,emis,q2,qsurfpal, & ! cell_area,albedodat,inertiedat,zmea,zstd,zsig, & ! zgam,zthe,oblipal,eccpal,tpalnew,adjustnew,phisfipal, & ! peri_daypal) endif else ! 'paleo' if (ngrid.ne.1) then write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin, & tsurf,tsoil,emis,q2,qsurf_hist) endif #endif endif ! end of 'paleo' endif ! end of 'lastcall' ! ----------------------------------------------------------------- ! WSTATS: Saving statistics ! ----------------------------------------------------------------- ! ("stats" stores and accumulates key variables in file "stats.nc" ! which can later be used to make the statistic files of the run: ! 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,"fluxsurf_lw", & "Thermal IR radiative flux to surface","W.m-2",2, & fluxsurf_lw) call wstats(ngrid,"fluxtop_lw", & "Thermal IR radiative flux to space","W.m-2",2, & fluxtop_lw) ! call wstats(ngrid,"fluxsurf_sw", & ! "Solar radiative flux to surface","W.m-2",2, & ! fluxsurf_sw_tot) ! call wstats(ngrid,"fluxtop_sw", & ! "Solar radiative flux to space","W.m-2",2, & ! fluxtop_sw_tot) call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) ! call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) ! call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) call wstats(ngrid,"p","Pressure","Pa",3,pplay) call wstats(ngrid,"emis","Emissivity","",2,emis) 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,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2) if (tracer) then do iq=1,nq call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',& 'kg m^-2',2,qsurf(1,iq) ) call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 'kg m^-2',2,qcol(1,iq) ) call wstats(ngrid,trim(noms(iq))//'_reff',trim(noms(iq))//'_reff',& 'm',3,reffrad(1,1,iq)) end do endif ! end of 'tracer' if(lastcall.and.callstats) then write (*,*) "Writing stats..." call mkstats(ierr) endif #ifndef MESOSCALE !------------------------------------------------------------------------------ ! OUTPUT in netcdf file "DIAGFI.NC", ! containing any variable for diagnostic ! ! Note 1 : output with period "ecritphy", set in "run.def" ! Note 2 : writediagfi can also be called from any other subroutine ! for any variable, but its preferable to keep all the ! calls in one place ... !------------------------------------------------------------------------------ call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi) call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi) call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi) call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi) call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf) call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps) call writediagfi(ngrid,"emis","Emissivity","",2,emis) !! Pluto outputs call writediagfi(ngrid,"dist_star","dist_star","AU",0,dist_star) if (fast) then call writediagfi(ngrid,"globave","surf press","Pa",0,globave) call writediagfi(ngrid,"fluxrad","fluxrad","W m-2",2,fluxrad) call writediagfi(ngrid,"fluxgrd","fluxgrd","W m-2",2,fluxgrd) call writediagfi(ngrid,"capcal","capcal","W.s m-2 K-1",2,capcal) ! call writediagfi(ngrid,"dplanck","dplanck","W.s m-2 K-1",2,dplanck) call writediagfi(ngrid,"tsoil","tsoil","K",3,tsoil) else if (check_physics_outputs) then ! Check the validity of updated fields at the end of the physics step call check_physics_fields("HERE physiq:", zt, zu, zv, pplev, zq) endif call writediagfi(ngrid,"temp","temperature","K",3,zt) call writediagfi(ngrid,"teta","potential temperature","K",3,zh) 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,"p","Pressure","Pa",3,pplay) endif ! Subsurface temperatures ! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf) ! call writediagsoil(ngrid,"temp","temperature","K",3,tsoil) ! Total energy balance diagnostics if(callrad)then call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent) call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1)) call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw) call writediagfi(ngrid,"shad","rings"," ", 2, fract) ! call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw) ! call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw) call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd) ! endif call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn) endif ! end of 'callrad' if(enertest) then if (calldifv) then call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2) call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux) ! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff) ! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff) ! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs) endif if (corrk) then call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw) call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw) endif if (generic_condensation) then call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE) call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation) endif endif ! end of 'enertest' ! Diagnostics of optical thickness ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19 if (diagdtau) then do nw=1,L_NSPECTV write(str2,'(i2.2)') nw call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw)) enddo do nw=1,L_NSPECTI write(str2,'(i2.2)') nw call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw)) enddo endif ! Temporary inclusions for heating diagnostics. call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw) call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw) call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad) call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn) ! For Debugging. !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat)) !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi) ! Output tracers. if (tracer) then ! call writediagfi(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc) do iq=1,nq call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & ! 'kg m^-2',2,qsurf_hist(1,iq) ) call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 'kg m^-2',2,qcol(1,iq) ) call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 'kg m^-2',2,qsurf(1,iq) ) if(generic_condensation)then call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic) call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac) call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic) endif ! call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col) enddo ! end of 'nq' loop !Pluto specific call writediagfi(ngrid,'n2_iceflux','n2_iceflux',"kg m^-2 s^-1",2,flusurf(1,igcm_n2) ) if (haze_radproffix)then call writediagfi(ngrid,'haze_reff','haze_reff','m',3,reffrad(1,1,1)) end if if (methane) then ! call writediagfi(ngrid,"rice_ch4","ch4 ice mass mean radius","m",3,rice_ch4) ! call writediagfi(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4) ! call writediagfi(ngrid,"qsat_ch4"," "," ",2,qsat_ch4) ! call writediagfi(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1) call writediagfi(ngrid,'ch4_iceflux','ch4_iceflux',& "kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) ) call writediagfi(ngrid,"vmr_ch4","vmr_ch4","%",2,vmr_ch4) if (.not.fast) then call writediagfi(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",3,zrho_ch4(:,:)) endif ! Tendancies call writediagfi(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",& 3,zdqc(:,:,igcm_ch4_gas)) call writediagfi(ngrid,"zdqdif_ch4","zdqdif ch4","",& 3,zdqdif(:,:,igcm_ch4_gas)) call writediagfi(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",& 2,zdqsdif(:,igcm_ch4_ice)) call writediagfi(ngrid,"zdqadj_ch4","zdqadj ch4","",& 3,zdqadj(:,:,igcm_ch4_gas)) if (sedimentation) then call writediagfi(ngrid,"zdqsed_ch4","zdqsed ch4","",& 3,zdqsed(:,:,igcm_ch4_gas)) call writediagfi(ngrid,"zdqssed_ch4","zdqssed ch4","",& 2,zdqssed(:,igcm_ch4_gas)) endif if (metcloud.and.(.not.fast)) then call writediagfi(ngrid,"zdtch4cloud","ch4 cloud","T s-1",& 3,zdtch4cloud) call writediagfi(ngrid,"zdqch4cloud","ch4 cloud","T s-1",& 3,zdqch4cloud(1,1,igcm_ch4_gas)) endif endif if (carbox) then ! call writediagfi(ngrid,"zdtcocloud","tendancy T cocloud","K",3,zdtcocloud) call writediagfi(ngrid,'co_iceflux','co_iceflux',& "kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) ) call writediagfi(ngrid,"vmr_co","vmr_co","%",2,vmr_co) if (.not.fast) THEN call writediagfi(ngrid,"zrho_co","zrho_co","kg.m-3",3,zrho_co(:,:)) endif endif if (haze) then ! call writediagfi(ngrid,"zrho_haze","zrho_haze","kg.m-3",3,zrho_haze(:,:)) call writediagfi(ngrid,"zdqrho_photprec","zdqrho_photprec",& "kg.m-3.s-1",3,zdqrho_photprec(:,:)) call writediagfi(ngrid,"zdqphot_prec","zdqphot_prec","",& 3,zdqphot_prec(:,:)) call writediagfi(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",& 3,zdqhaze(:,:,igcm_ch4_gas)) call writediagfi(ngrid,"zdqhaze_prec","zdqhaze_prec","",& 3,zdqhaze(:,:,igcm_prec_haze)) if (igcm_haze.ne.0) then call writediagfi(ngrid,"zdqhaze_haze","zdqhaze_haze","",& 3,zdqhaze(:,:,igcm_haze)) if (sedimentation) then call writediagfi(ngrid,"zdqssed_haze","zdqssed haze",& "kg/m2/s",2,zdqssed(:,igcm_haze)) endif endif call writediagfi(ngrid,"zdqphot_ch4","zdqphot_ch4","",& 3,zdqphot_ch4(:,:)) call writediagfi(ngrid,"zdqconv_prec","zdqconv_prec","",& 3,zdqconv_prec(:,:)) ! call writediagfi(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s", ! & 2,zdqhaze_col(:)) endif if (aerohaze) then call writediagfi(ngrid,"tau_col",& "Total aerosol optical depth","opacity",2,tau_col) endif endif ! end of 'tracer' ! Output spectrum. if(specOLR.and.corrk)then call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu) call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu) call writediagspecVI(ngrid,"GSR3D","GSR(lon,lat,band)","W/m^2/cm^-1",3,GSR_nu) endif #else comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer) comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer) comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid) if (.not.calldifv) comm_LATENT_HF(:)=0.0 if ((tracer).and.(generic_condensation)) then ! If you have set generic_condensation (and not water) and you have set several GCS ! then the outputs given to WRF will be only the ones for the last generic tracer ! (Because it is rewritten every tracer in the loop) ! WRF can take only one moist tracer do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer reffrad_generic_zeros_for_wrf(:,:) = 1. comm_CLOUDFRAC(1:ngrid,1:nlayer) = cloudfrac(1:ngrid,1:nlayer) comm_TOTCLOUDFRAC(1:ngrid) = totcloudfrac(1:ngrid) !?????? comm_SURFRAIN(1:ngrid) = zdqsrain_generic(1:ngrid,iq) comm_DQVAP(1:ngrid,1:nlayer) = pdq(1:ngrid,1:nlayer,igcm_generic_gas) comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_generic_ice) ! comm_H2OICE_REFF(1:ngrid,1:nlayer) = reffrad_generic_zeros_for_wrf(1:ngrid,1:nlayer) ! for now zeros ! !comm_H2OICE_REFF(1:ngrid,1:nlayer) = 0*zdtrain_generic(1:ngrid,1:nlayer) ! for now zeros ! comm_REEVAP(1:ngrid) = reevap_precip_generic(1:ngrid,iq) comm_DTRAIN(1:ngrid,1:nlayer) = zdtrain_generic(1:ngrid,1:nlayer) comm_DTLSC(1:ngrid,1:nlayer) = dt_generic_condensation(1:ngrid,1:nlayer) comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq) endif end do ! do iq=1,nq loop on tracers else comm_CLOUDFRAC(1:ngrid,1:nlayer)=0. comm_TOTCLOUDFRAC(1:ngrid)=0. comm_SURFRAIN(1:ngrid)=0. comm_DQVAP(1:ngrid,1:nlayer)=0. comm_DQICE(1:ngrid,1:nlayer)=0. ! comm_H2OICE_REFF(1:ngrid,1:nlayer)=0. comm_REEVAP(1:ngrid)=0. comm_DTRAIN(1:ngrid,1:nlayer)=0. comm_DTLSC(1:ngrid,1:nlayer)=0. comm_RH(1:ngrid,1:nlayer)=0. endif ! if water, if generic_condensation, else comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid) comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid) comm_FLUXTOP_LW(1:ngrid)=fluxtop_lw(1:ngrid) comm_FLUXSURF_SW(1:ngrid)=fluxsurf_sw(1:ngrid) comm_FLUXSURF_LW(1:ngrid)=fluxsurf_lw(1:ngrid) comm_FLXGRD(1:ngrid)=fluxgrd(1:ngrid) sensibFlux(1:ngrid) = zflubid(1:ngrid) - capcal(1:ngrid)*zdtsdif(1:ngrid) !!! ???? comm_HR_DYN(1:ngrid,1:nlayer) = zdtdyn(1:ngrid,1:nlayer) #endif ! 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) CALL send_xios_field("ps",ps) CALL send_xios_field("area",cell_area) CALL send_xios_field("p",pplay) CALL send_xios_field("temperature",zt) CALL send_xios_field("u",zu) CALL send_xios_field("v",zv) call send_xios_field("w",pw) CALL send_xios_field("ISR",fluxtop_dn) CALL send_xios_field("OLR",fluxtop_lw) CALL send_xios_field("ASR",fluxabs_sw) if (specOLR .and. corrk) then call send_xios_field("OLR3D",OLR_nu) call send_xios_field("IR_Bandwidth",DWNI) call send_xios_field("VI_Bandwidth",DWNV) call send_xios_field("OSR3D",OSR_nu) call send_xios_field("GSR3D",GSR_nu) endif if (lastcall.and.is_omp_master) then write(*,*) "physiq: 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, pplev, zq) endif icount=icount+1 end subroutine physiq end module physiq_mod