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 watercommon_h, only : RLVTT, Psat_water,epsi,su_watercycle, RV, T_h2o_ice_liq use generic_cloud_common_h, only : epsi_generic, Psat_generic use thermcell_mod, only: init_thermcell_mod use gases_h, only: gnom, gfrac, ngasmx use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV use suaer_corrk_mod, only: suaer_corrk use radii_mod, only: h2o_reffrad, co2_reffrad use aerosol_mod, only: iniaerosol, iaero_co2, iaero_h2o 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 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, & alpha_lift, alpha_devil, qextrhor, & igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, & igcm_co2_ice, nesp, is_chim, is_condensable,constants_epsi_generic use time_phylmdz_mod, only: ecritphy, iphysiq, nday use phyetat0_mod, only: phyetat0,tab_cntrl_mod use surfini_mod, only: surfini use wstats_mod, only: callstats, wstats, mkstats use phyredem, only: physdem0, physdem1 !! use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, & !! noceanmx !! use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, & !! capcalocean, capcalseaice,capcalsno, nslay, & !! ocean_slab_final !! use surf_heat_transp_mod,only :init_masquv use ocean_slab_mod, only :ocean_slab_init, ocean_slab_noice, ocean_slab_ice, & ocean_slab_frac, ocean_slab_get_vars, & capcalocean, capcalseaice,capcalsno, nslay, knon, & ocean_slab_final 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 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, & calllott_nonoro, callrad, callsoil, nosurf, & calltherm, CLFvarying, co2cond, corrk, diagdtau, & diurnal, enertest, fat1au, flatten, j2, & hydrology, intheat, iradia, kastprof, & lwrite, mass_redistrib, massplanet, meanOLR, & nearco2cond, newtonian, noseason_day, oblate, & ok_slab_ocean, photochem, rings_shadow, rmean, & season, sedimentation,generic_condensation, & specOLR, callthermos, callvolcano, & startphy_file, testradtimes, tlocked, & tracer, UseTurbDiff, water, watercond, & waterrain, generic_rain, global1d, szangle, & moistadjustment, moistadjustment_generic, varspec, varspec_data, nvarlayer use generic_tracer_index_mod, only: generic_tracer_index use nonoro_gwd_ran_mod, only: nonoro_gwd_ran 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 conduction_mod, only: conduction use volcano_mod, only: volcano use pindex_mod, only: pindex use convadj_mod, only: convadj use vdifc_mod, only: vdifc use turbdiff_mod, only: turbdiff use turb_mod, only : q2,sensibFlux,turb_resolved use callsedim_mod, only: callsedim use mass_redistribution_mod, only: mass_redistribution use condensation_generic_mod, only: condensation_generic use datafile_mod, only: datadir use newton_cooling_hotJ, only: newtcool_MOCHA ! LT, adding for MOCHA protocol #ifndef MESOSCALE USE vertical_layers_mod, ONLY: ap,bp,aps,bps,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_SURFRAIN,comm_REEVAP,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_H2OICE_REFF, & 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 use write_output_mod, only: write_output 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.c 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 CO2). ! ! 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 l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice logical call_ice_vap_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 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 real f(ngrid) ! Mass flux norm real fm(ngrid,nlayer+1) ! Mass flux real fm_bis(ngrid,nlayer) ! Recasted fm real entr(ngrid,nlayer) ! Entrainment real detr(ngrid,nlayer) ! Detrainment real dqevap(ngrid,nlayer,nq) ! water tracer mass mixing ratio variations due to evaporation real dtevap(ngrid,nlayer) ! temperature variation due to evaporation real zqtherm(ngrid,nlayer,nq) ! vapor mass mixing ratio after evaporation real zttherm(ngrid,nlayer) ! temperature after evaporation real fraca(ngrid,nlayer+1) ! Fraction of the surface that plumes occupies real zw2(ngrid,nlayer+1) ! Vertical speed real zw2_bis(ngrid,nlayer) ! Recasted zw2 ! 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_co2 routine. real zdtsdif(ngrid) ! Turbdiff/vdifc routines. real zdtsurf_hyd(ngrid) ! Hydrol routine. ! For Thermosphere : (K/s) real zdtconduc(ngrid,nlayer) ! 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_co2 routine. real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines. real zdttherm(ngrid,nlayer) ! Calltherm routine. real zdtmr(ngrid,nlayer) ! Mass_redistribution routine. real zdtrain(ngrid,nlayer) ! Rain routine. real zdtrain_generic(ngrid,nlayer) ! Rain_generic routine. real dtmoist(ngrid,nlayer) ! Moistadj routine. real dt_ekman(ngrid,nslay), dt_hdiff(ngrid,nslay), dt_gm(ngrid,nslay) ! Slab_ocean 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_co2 routine. real zdqsdif(ngrid,nq) ! Turbdiff/vdifc routines. real zdqssed(ngrid,nq) ! Callsedim routine. real zdqsurfmr(ngrid,nq) ! Mass_redistribution routine. real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine. real zdqsrain_generic(ngrid,nq), zdqssnow_generic(ngrid,nq) ! Rain_generic routine. real dqs_hyd(ngrid,nq) ! Hydrol routine. real reevap_precip(ngrid) ! re-evaporation flux of precipitation (integrated over the atmospheric column) real reevap_precip_generic(ngrid,nq) ! For Tracers : (kg/kg_of_air/s) real zdqc(ngrid,nlayer,nq) ! Condense_co2 routine. real zdqadj(ngrid,nlayer,nq) ! Convadj routine. real zdqdif(ngrid,nlayer,nq) ! Turbdiff/vdifc routines. real zdqevap(ngrid,nlayer) ! Turbdiff routine. real zdqtherm(ngrid,nlayer,nq) ! Calltherm routine. real zdqsed(ngrid,nlayer,nq) ! Callsedim routine. real zdqmr(ngrid,nlayer,nq) ! Mass_redistribution routine. real zdqrain(ngrid,nlayer,nq) ! Rain routine. real dqmoist(ngrid,nlayer,nq) ! Moistadj routine. real dqvaplscale(ngrid,nlayer) ! Largescale routine. real dqcldlscale(ngrid,nlayer) ! Largescale routine. real dqvaplscale_generic(ngrid,nlayer,nq) ! condensation_generic routine. real dqcldlscale_generic(ngrid,nlayer,nq) ! condensation_generic routine. real dq_rain_generic_vap(ngrid,nlayer,nq) ! rain_generic routine real dq_rain_generic_cld(ngrid,nlayer,nq) ! rain_generic routine REAL,allocatable,save :: zdqchim(:,:,:) ! Calchim_asis routine REAL,allocatable,save :: zdqschim(:,:) ! Calchim_asis routine !$OMP THREADPRIVATE(zdqchim,zdqschim) real zdqvolc(ngrid,nlayer,nq) ! injection by volcanoes (kg/kg_of_air/s) REAL array_zero1(ngrid) REAL array_zero2(ngrid,nlayer) ! For Winds : (m/s/s) real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer) ! Convadj routine. real zdutherm(ngrid,nlayer), zdvtherm(ngrid,nlayer) ! Calltherm 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. ! 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 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 tice(ngrid) ! Sea ice temperature 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). ! 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 h2otot ! Total Amount of water. For diagnostic. real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic. real dWtot, dWtot_tmp, dWtots, dWtots_tmp logical,save :: watertest !$OMP THREADPRIVATE(watertest) real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air). real RH(ngrid,nlayer) ! Relative humidity. real psat_tmp 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 logical clearsky ! For double radiative transfer call. By BC ! For Clear Sky Case. real fluxsurfabs_sw1(ngrid) ! For SW/LW flux. real tau_col1(ngrid) ! For aerosol optical depth diagnostic. real OLR_nu1(ngrid,L_NSPECTI) ! Clear sky TOA LW radiation in each IR band real OSR_nu1(ngrid,L_NSPECTV) ! Clear sky TOA SW radiation in each VI band real GSR_nu1(ngrid,L_NSPECTV) ! Clear sky Surface SW radiation in each VI band real int_dtaui1(ngrid,nlayer,L_NSPECTI),int_dtauv1(ngrid,nlayer,L_NSPECTV) ! For optical thickness diagnostics. real tf, ntf real net_fluxsurf_lw(ngrid) ! net longwave flux at the surface (for diagnostics only) real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW ! For fixed variable molar mass real, dimension(:),allocatable,save :: p_var real, dimension(:),allocatable,save :: mu_var real, dimension(:,:),allocatable,save :: frac_var !$OMP THREADPRIVATE(p_var,mu_var,frac_var) real :: mu_vari(nlayer) real :: muvari(ngrid,nlayer) integer ios,k character*100 dt_file 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 ! 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 zdtconduc(:,:) = 0.0 ! Initialize fixed variable mu ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(varspec) then IF (.NOT.ALLOCATED(p_var)) ALLOCATE(p_var(nvarlayer)) IF (.NOT.ALLOCATED(mu_var)) ALLOCATE(mu_var(nvarlayer)) IF (.NOT.ALLOCATED(frac_var)) ALLOCATE(frac_var(nvarlayer,ngasmx)) p_var = 0.0 mu_var = 0.0 frac_var = 0.0 muvari = 0.0 dt_file=TRIM(varspec_data) open(33,file=dt_file,form='formatted',status='old',iostat=ios) if (ios.ne.0) then ! file not found write(*,*) 'Error from varspec' write(*,*) 'Data file ',trim(varspec_data),' not found.' write(*,*) 'Check that the data is in your run repository.' call abort_physic !a verifier else do k=1,nvarlayer read(33,*) p_var(k), mu_var(k),frac_var(k,1:ngasmx) !The order of columns in frac_var must correspond to order of molecules gases.def !The format of your file must be: ! pressure(k) molar_mass(k), molar_fraction_of_gas_1(k), molar_fraction_of_gas_2(k),..., molar_fraction_of_gas_ngasmx(k) enddo endif close(33) else IF (.NOT.ALLOCATED(p_var)) ALLOCATE(p_var(nlayer)) IF (.NOT.ALLOCATED(mu_var)) ALLOCATE(mu_var(nlayer)) IF (.NOT.ALLOCATED(frac_var)) ALLOCATE(frac_var(nlayer,ngasmx)) endif ! 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 call ini_conc_mod(ngrid,nlayer) IF (.NOT.ALLOCATED(zdqchim)) ALLOCATE(zdqchim(ngrid,nlayer,nesp)) IF (.NOT.ALLOCATED(zdqschim)) ALLOCATE(zdqschim(ngrid,nesp)) endif 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,albedo, & q2,qsurf,cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice) #else day_ini = pday #endif #ifndef MESOSCALE if (.not.startphy_file) then ! additional "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. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,& albedo_snow_SPECTV,albedo_co2_ice_SPECTV) ! Initialize orbital calculation. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) if(tlocked)then print*,'Planet is tidally locked at resonance n=',nres print*,'Make sure you have the right rotation rate!!!' endif ! Initialize oceanic variables. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (ok_slab_ocean)then !! call ocean_slab_init(ngrid,ptimestep, tslab, & !! sea_ice, pctsrf_sic) !! call ini_surf_heat_transp_mod() knindex(:) = 0 knon = 0 do ig=1,ngrid zmasq(ig)=1 ! knindex(ig) = ig if (nint(rnat(ig)).eq.0) then ! rnat = 0 (ocean), 1 (continent) zmasq(ig)=0 knon=knon+1 knindex(knon)=ig endif enddo call ocean_slab_init(ptimestep, pctsrf_sic, tslab, tice, sea_ice, zmasq) !! CALL init_masquv(ngrid,zmasq) ! ensure that "slab" variables match "physiq" ones call ocean_slab_get_vars(ngrid, tslab, tice, sea_ice, flux_g, & dt_hdiff, dt_ekman, dt_gm) endif ! end of 'ok_slab_ocean'. ! Initialize soil. ! ~~~~~~~~~~~~~~~~ if (callsoil) then call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, & ptimestep,tsurf,tsoil,capcal,fluxgrd) if (ok_slab_ocean) then ! mimic what is done during a regular time step do ig=1,ngrid fluxgrdocean(ig)=fluxgrd(ig) if (nint(rnat(ig)).eq.0) then capcal(ig)=capcalocean fluxgrd(ig)=0. ! Dividing by cell area to have flux in W/m2 fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1)+dt_gm(ig,1))/cell_area(ig) do iq=1,nsoilmx tsoil(ig,iq)=tsurf(ig) enddo if (pctsrf_sic(ig).gt.0.5) then capcal(ig)=capcalseaice if (qsurf(ig,igcm_h2o_ice).gt.0.) then capcal(ig)=capcalsno endif endif ! of if (pctsrf_sic(ig).gt.0.5) endif ! of if (nint(rnat(ig)).eq.0) enddo ! of do ig=1,ngrid endif ! end of 'ok_slab_ocean'. 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 ! Here is defined the type of the surface : Continent or Ocean. ! BC2014 : This is now already done in newstart. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (.not.ok_slab_ocean) then rnat(:)=1. if (.not.nosurf) then ! inertiedat only defined if there is a surface do ig=1,ngrid if(inertiedat(ig,1).gt.1.E4)then rnat(ig)=0 endif enddo print*,'WARNING! Surface type currently decided by surface inertia' print*,'This should be improved e.g. in newstart.F' endif endif ! end of 'ok_slab_ocean'. ! Initialize surface history variable. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ qsurf_hist(:,:)=qsurf(:,:) ! Initialize variable for dynamical heating and zonal wind tendency diagnostic ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ztprevious(:,:)=pt(:,:) zuprevious(:,:)=pu(:,:) ! Set temperature just above condensation temperature (for Early Mars) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(nearco2cond) 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 watertest=.false. if(water)then ! Initialize variables for water cycle if(enertest)then watertest = .true. endif endif ! 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 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (calltherm) then CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV) endif call su_watercycle ! even if we don't have a water cycle, we might ! need epsi for the wvp definitions in callcorrk.F ! or RETV, RLvCp for the thermal plume model #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 #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 print*, "physiq_mod: Correlated-k data base folder:",trim(datadir) call getin_p("corrkdir",corrkdir) 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. call suaer_corrk ! Set up aerosol optical properties. endif ! 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 write(*,*) "physiq: end of firstcall" endif ! end of '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 ! ------------------------------------------------------ ! 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.nearco2cond ) 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 net_fluxsurf_lw(1:ngrid) = 0.0 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 call orbite(zls,dist_star,declin,right_ascen) if (tlocked) then zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi) elseif (diurnal) then zlss=-2.*pi*(zday-.5) else if(diurnal .eqv. .false.) then zlss=9999. endif ! Compute variations of g with latitude (oblate case). ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (oblate .eqv. .false.) then glat(:) = g else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)' call abort else gmplanet = MassPlanet*grav*1e24 do ig=1,ngrid glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig)))) end do endif ! 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 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 ! ---------------------------------------------------------------- ! Compute mean mass, cp, and R ! -------------------------------- #ifndef MESOSCALE if(photochem) then call concentrations(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep) endif #endif !--------------------------------- ! II. Compute radiative tendencies !--------------------------------- ! Compute local stellar zenith angles ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (tlocked) then ! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb ztim1=SIN(declin) ztim2=COS(declin)*COS(zlss) ztim3=COS(declin)*SIN(zlss) call stelang(ngrid,sinlon,coslon,sinlat,coslat, & ztim1,ztim2,ztim3,mu0,fract, flatten) elseif (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, flatten) else if(diurnal .eqv. .false.) then call mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad,flatten) ! 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) endif endif if (callrad) then if( mod(icount-1,iradia).eq.0) then ! Eclipse incoming sunlight (e.g. Saturn ring shadowing). if(rings_shadow) then call call_rings(ngrid, ptime, pday, diurnal) endif 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 ! take into account water effect on mean molecular weight muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) elseif(generic_condensation) then ! take into account generic condensable specie (GCS) effect on mean molecular weight do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) if (call_ice_vap_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_vap)) muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap)) endif end do ! do iq=1,nq loop on tracers elseif(varspec) then !take into account fixed variable mean molecular weight do ig=1,ngrid call pindex(p_var,mu_var,pplay(ig,:),nvarlayer,nlayer,mu_vari) muvari(ig,:) = mu_vari enddo muvar(1:ngrid,1:nlayer) = muvari(1:ngrid,1:nlayer) muvar(1:ngrid,nlayer+1) = muvari(1:ngrid,nlayer) else muvar(1:ngrid,1:nlayer+1)=mugaz endif if(ok_slab_ocean) then tsurf2(:)=tsurf(:) do ig=1,ngrid if (nint(rnat(ig))==0) then tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25 endif enddo endif !(ok_slab_ocean) ! standard callcorrk clearsky=.false. call callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & 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, & clearsky,p_var,frac_var,firstcall,lastcall) !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 ! Option to call scheme once more for clear regions if(CLFvarying)then ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ... clearsky=.true. call callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt, & tsurf,fract,dist_star,aerosol,muvar, & zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1, & fluxsurfabs_sw1,fluxtop_lw1, & fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,GSR_nu1, & int_dtaui1,int_dtauv1, & tau_col1,cloudfrac,totcloudfrac, & clearsky,p_var,frac_var,firstcall,lastcall) clearsky = .false. ! just in case. ! Sum the fluxes and heating rates from cloudy/clear cases do ig=1,ngrid tf=totcloudfrac(ig) ntf=1.-tf fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig) fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig) albedo_equivalent(ig) = ntf*albedo_equivalent1(ig) + tf*albedo_equivalent(ig) fluxsurfabs_sw(ig) = ntf*fluxsurfabs_sw1(ig) + tf*fluxsurfabs_sw(ig) fluxtop_lw(ig) = ntf*fluxtop_lw1(ig) + tf*fluxtop_lw(ig) fluxabs_sw(ig) = ntf*fluxabs_sw1(ig) + tf*fluxabs_sw(ig) tau_col(ig) = ntf*tau_col1(ig) + tf*tau_col(ig) zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer) zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer) OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV) GSR_nu(ig,1:L_NSPECTV) = ntf*GSR_nu1(ig,1:L_NSPECTV) + tf*GSR_nu(ig,1:L_NSPECTV) OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI) if (diagdtau) then int_dtauv(ig,:,1:L_NSPECTV) = ntf*int_dtauv1(ig,:,1:L_NSPECTV) + tf*int_dtauv(ig,:,1:L_NSPECTV) int_dtaui(ig,:,1:L_NSPECTI) = ntf*int_dtaui1(ig,:,1:L_NSPECTI) + tf*int_dtaui(ig,:,1:L_NSPECTI) endif enddo endif ! end of CLFvarying. if(ok_slab_ocean) then tsurf(:)=tsurf2(:) endif ! 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) net_fluxsurf_lw(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid) !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 elseif(newtonian)then ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! II.b Call Newtonian cooling scheme ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall) call newtcool_MOCHA(ngrid,nlayer,coslon,coslat,pt,pplay,firstcall,lastcall,dtrad) zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep ! e.g. surface becomes proxy for 1st atmospheric layer ? fluxrad_sky(1:ngrid) = 0. ! otherwise, not initialized! else ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! II.c 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' fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1)) 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) net_fluxsurf_lw(1:ngrid)=net_fluxsurf_lw(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) ! -------------------------------------------- ! III. Vertical diffusion (turbulent mixing) : ! -------------------------------------------- if (calldifv) then zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid) ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff. if (UseTurbDiff) then call turbdiff(ngrid,nlayer,nq,rnat, & 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 zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) call vdifc(ngrid,nlayer,nq,rnat,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 flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid)) endif 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) ! 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. if(watertest.and.water)then call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp) call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots_tmp) do ig = 1, ngrid vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap)) enddo call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*cell_area(:)*ptimestep/totarea_planet,dWtots) dWtot = dWtot + dWtot_tmp dWtots = dWtots + dWtots_tmp do ig = 1, ngrid vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice)) enddo call planetwide_maxval(vdifcncons(:),nconsMAX) if (is_master) then print*,'---------------------------------------------------------------' print*,'In difv atmospheric water change =',dWtot,' kg m-2' print*,'In difv surface water change =',dWtots,' kg m-2' print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2' print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1' endif endif ! end of 'watertest' !------------------------- else ! calldifv if(.not.newtonian)then zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid) endif endif ! end of 'calldifv' !------------------- ! IV. Convection : !------------------- ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ ! IV.a Thermal plume model : ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (calltherm) THEN ! AB: We need to evaporate ice before calling thermcell_main. IF (water) THEN CALL evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,zqtherm,zttherm) !ALS24 adapats for generic tracer ELSEIF (generic_condensation .and. .not. water) THEN DO iq=1,nq CALL generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) CALL evap_generic(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,igcm_generic_vap, & igcm_generic_ice,dqevap,dtevap,zqtherm,zttherm) ENDDO ELSE zttherm(:,:) = pt(:,:) + pdt(:,:) * ptimestep zqtherm(:,:,:) = pq(:,:,:) + pdq(:,:,:) * ptimestep ENDIF CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall, & pplay, pplev, pphi, zpopsk, & pu, pv, zttherm, zqtherm, & zdutherm, zdvtherm, zdttherm, zdqtherm, & fm, entr, detr, zw2, fraca) pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdutherm(1:ngrid,1:nlayer) pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvtherm(1:ngrid,1:nlayer) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer) IF (tracer) THEN pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq) ENDIF ENDIF ! end of 'calltherm' ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 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 if(watertest)then call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp) do ig = 1, ngrid cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap)) enddo call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) dWtot = dWtot + dWtot_tmp do ig = 1, ngrid cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice)) enddo call planetwide_maxval(cadjncons(:),nconsMAX) if (is_master) then print*,'In convadj atmospheric water change =',dWtot,' kg m-2' print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1' endif endif ! end of 'watertest' endif ! end of 'calladj' !---------------------------------------------- ! Non orographic Gravity Waves: !--------------------------------------------- IF (calllott_nonoro) THEN CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, & cpnew, rnew, & pplay, & 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) print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin) print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin) ENDIF ! of IF (calllott_nonoro) !----------------------------------------------- ! V. Carbon dioxide condensation-sublimation : !----------------------------------------------- if (co2cond) then if (.not.tracer) then print*,'We need a CO2 ice tracer to condense CO2' call abort endif call condense_co2(ngrid,nlayer,nq,ptimestep, & capcal,pplay,pplev,tsurf,pt, & pdt,zdtsurf,pq,pdq, & qsurf,zdqsurfc,albedo,emis, & albedo_bareground,albedo_co2_ice_SPECTV, & zdtc,zdtsurfc,pdpsrf,zdqc) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(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_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid) ! 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 co2cloud atmospheric energy change =',dEtot,' W m-2' print*,'In co2cloud surface energy change =',dEtots,' W m-2' endif endif endif ! end of 'co2cond' !--------------------------------------------- ! VI. Specific parameterizations for tracers !--------------------------------------------- if (tracer) then ! ----------------------------------------------------- ! VI.0. Volcanoes injecting tracers in the atmosphere ! ----------------------------------------------------- if (callvolcano) then call volcano(ngrid,nlayer,pplev,pu,pv,pt,zdqvolc,nq,massarea,& zday,ptimestep) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + & zdqvolc(1:ngrid,1:nlayer,1:nq) endif ! --------------------- ! VI.1. Water and ice ! --------------------- if (water) then ! Water ice condensation in the atmosphere if(watercond.and.(RLVTT.gt.1.e-8))then if ((.not.calltherm).and.moistadjustment) then dqmoist(1:ngrid,1:nlayer,1:nq)=0. dtmoist(1:ngrid,1:nlayer)=0. ! Moist Convective Adjustment. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man) pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) & + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) & + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer) ! Test energy conservation. if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot) call planetwide_maxval(dtmoist(:,:),dtmoist_max) call planetwide_minval(dtmoist(:,:),dtmoist_min) madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:) do ig=1,ngrid madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:)) enddo if (is_master) then print*,'In moistadj atmospheric energy change =',dEtot,' W m-2' print*,'In moistadj MAX atmospheric energy change =',dtmoist_max*ptimestep,'K/step' print*,'In moistadj MIN atmospheric energy change =',dtmoist_min*ptimestep,'K/step' endif call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+ & massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) if (is_master) print*,'In moistadj atmospheric water change =',dWtot,' kg m-2' endif ! end of 'enertest' else ! rneb_man, dqmoist are output of moistadj, but used later on ! so set them to zero if moistadj is not called rneb_man(:,:)=0 dqmoist(:,:,:)=0 endif ! end of '(.not.calltherm).and.moistadjustment' ! Large scale condensation/evaporation. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer) pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer) pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer) ! Test energy conservation. if(enertest)then lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:) do ig=1,ngrid lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:)) enddo call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot) if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2' ! Test water conservation. call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+ & SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot) if (is_master) print*,'In largescale atmospheric water change =',dWtot,' kg m-2' endif ! end of 'enertest' ! Compute cloud fraction. do l = 1, nlayer do ig=1,ngrid cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) enddo enddo endif ! end of 'watercond' ! Water ice / liquid precipitation. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0 !JL18 need to do that everytimestep if mass redis is on. if(waterrain)then zdqsrain(1:ngrid) = 0.0 zdqssnow(1:ngrid) = 0.0 call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq, & zdtrain,zdqrain,zdqsrain,zdqssnow,reevap_precip,cloudfrac) pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) & + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) & + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer) dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid) ! Test energy conservation. if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot) if (is_master) print*,'In rain atmospheric T energy change =',dEtot,' W m-2' call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp) call planetwide_sumval(cell_area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot) dItot = dItot + dItot_tmp call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp) call planetwide_sumval(cell_area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot) dVtot = dVtot + dVtot_tmp dEtot = dItot + dVtot if (is_master) then print*,'In rain dItot =',dItot,' W m-2' print*,'In rain dVtot =',dVtot,' W m-2' print*,'In rain atmospheric L energy change =',dEtot,' W m-2' endif ! Test water conservation call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+ & massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*cell_area(:)*ptimestep/totarea_planet,dWtots) if (is_master) then print*,'In rain atmospheric water change =',dWtot,' kg m-2' print*,'In rain surface water change =',dWtots,' kg m-2' print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2' endif endif ! end of 'enertest' end if ! enf of 'waterrain' end if ! end of 'water' ! ------------------------- ! VI.2. Photochemistry ! ------------------------- #ifndef MESOSCALE IF (photochem) then DO ig=1,ngrid array_zero1(ig)=0.0 DO l=1,nlayer array_zero2(ig,l)=0. ENDDO ENDDO call calchim_asis(ngrid,nlayer,nq, & ptimestep,pplay,pplev,pt,pdt,dist_star,mu0, & fract,zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, & array_zero1,array_zero1, & pu,pdu,pv,pdv,array_zero2,array_zero2,icount,zdtchim) ! increment values of tracers: iesp = 0 DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry ! tracers is zero anyways ! September 2020: flag is_chim to increment only on chemical species IF (is_chim(iq)==1) THEN iesp = iesp + 1 DO l=1,nlayer DO ig=1,ngrid pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iesp) ENDDO ENDDO ENDIF ENDDO ! of DO iq=1,nq ! 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 ! increment values of temperature: pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtchim(1:ngrid,1:nlayer) END IF ! of IF (photochem) #endif ! ------------------------------------------ ! VI.3.a Generic Condensable Species (GCS) ! ------------------------------------------ if (generic_condensation) then if ((.not.water).and.(.not.calltherm).and.moistadjustment_generic) then write(*,*) "moist adjustment for GCS" dqmoist(1:ngrid,1:nlayer,1:nq)=0. dtmoist(1:ngrid,1:nlayer)=0. ! Moist Convective Adjustment. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call moistadj_generic(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) & + dqmoist(1:ngrid,1:nlayer,1:nq) pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer) ! Test energy conservation. if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot) call planetwide_maxval(dtmoist(:,:),dtmoist_max) call planetwide_minval(dtmoist(:,:),dtmoist_min) madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:) do ig=1,ngrid madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:)) enddo if (is_master) then print*,'In moistadj_generic atmospheric energy change =',dEtot,' W m-2' print*,'In moistadj_generic MAX atmospheric energy change =',dtmoist_max*ptimestep,'K/step' print*,'In moistadj_generic MIN atmospheric energy change =',dtmoist_min*ptimestep,'K/step' endif ! igcm_generic_vap & igcm_generic_ice are not declared below call planetwide_sumval(massarea(:,:)*dqmoist(:,:,1)*ptimestep/totarea_planet+ & massarea(:,:)*dqmoist(:,:,2)*ptimestep/totarea_planet,dWtot) if (is_master) print*,'In moistadj_generic atmospheric GCS change =',dWtot,' kg m-2' endif ! end of 'enertest' else ! rneb_man, dqmoist are output of moistadj, but used later on ! so set them to zero if moistadj is not called rneb_man(:,:)=0 dqmoist(:,:,:)=0 endif ! end of '(.not.calltherm).and.moistadjustment_generic' 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 if (.not. water) then ! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction ! Water is the priority ! If you have set water and generic_condensation, then cloudfrac will be water cloudfrac ! ! If you have set generic_condensation (and not water) and you have set several GCS ! then cloudfrac will be only the cloudfrac of the last generic tracer ! (Because it is rewritten every tracer in the loop) ! ! Maybe one should create a cloudfrac_generic(ngrid,nlayer,nq) with 3 dimensions, the last one for tracers ! Let's loop on tracers cloudfrac(:,:)=0.0 do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) if (call_ice_vap_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 !Generic Rain if (generic_rain) then zdqsrain_generic(1:ngrid,1:nq) = 0.0 zdqssnow_generic(1:ngrid,1:nq) = 0.0 call rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq, & zdtrain_generic,dq_rain_generic_vap,dq_rain_generic_cld, & zdqsrain_generic,zdqssnow_generic,reevap_precip_generic,cloudfrac) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) & + dq_rain_generic_vap(1:ngrid,1:nlayer,1:nq) ! only the parts with indexes of generic vapor tracers are filled in dq_rain_generic_vap, other parts are 0. pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) & + dq_rain_generic_cld(1:ngrid,1:nlayer,1:nq) ! only the parts with indexes of generic ice(cloud) tracers are filled in dq_rain_generic_cld, other parts are 0. pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain_generic(1:ngrid,1:nlayer) dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq) ! Test energy conservation. if(enertest)then call planetwide_sumval(cpp*massarea(:,:)*zdtrain_generic(:,:)/totarea_planet,dEtot) if (is_master) print*,'In rain_generic atmospheric T energy change =',dEtot,' W m-2' ! Test conservationfor each generic condensable specie (GCS) tracer ! Let's loop on tracers do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic) if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer call planetwide_sumval(massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dItot_tmp) call planetwide_sumval(cell_area(:)*zdqssnow_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dItot) dItot = dItot + dItot_tmp call planetwide_sumval(massarea(:,:)*dq_rain_generic_vap(:,:,igcm_generic_vap)/totarea_planet*RLVTT/cpp,dVtot_tmp) call planetwide_sumval(cell_area(:)*zdqsrain_generic(:,igcm_generic_ice)/totarea_planet*RLVTT/cpp,dVtot) dVtot = dVtot + dVtot_tmp dEtot = dItot + dVtot if (is_master) then print*,'In rain_generic dItot =',dItot,' W m-2' print*,'In rain_generic dVtot =',dVtot,' W m-2' print*,'In rain_generic atmospheric L energy change =',dEtot,' W m-2' endif call planetwide_sumval(massarea(:,:)*dq_rain_generic_vap(:,:,igcm_generic_vap)*ptimestep/totarea_planet+ & massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)*ptimestep/totarea_planet,dWtot) call planetwide_sumval((zdqsrain_generic(:,igcm_generic_ice)+zdqssnow_generic(:,igcm_generic_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots) if (is_master) then print*,'In rain_generic atmospheric generic tracer change =',dWtot,' kg m-2' print*,'In rain_generic surface generic tracer change =',dWtots,' kg m-2' print*,'In rain_generic non-cons factor =',dWtot+dWtots,' kg m-2' endif endif end do ! do iq=1,nq loop on tracers endif ! end of 'enertest' endif !generic_rain ! ------------------------- ! VI.3.b. Sedimentation. ! ------------------------- if (sedimentation) then zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0 zdqssed(1:ngrid,1:nq) = 0.0 if(watertest)then iq=igcm_h2o_ice call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot) call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots) if (is_master) then print*,'Before sedim pq =',dWtot,' kg m-2' print*,'Before sedim pdq =',dWtots,' kg m-2' endif endif call callsedim(ngrid,nlayer,ptimestep, & pplev,zzlev,pt,pdt,pq,pdq, & zdqsed,zdqssed,nq) if(watertest)then iq=igcm_h2o_ice call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot) call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots) if (is_master) then print*,'After sedim pq =',dWtot,' kg m-2' print*,'After sedim pdq =',dWtots,' kg m-2' endif 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) ! Test water conservation if(watertest)then call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot) call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots) if (is_master) then print*,'In sedim atmospheric ice change =',dWtot,' kg m-2' print*,'In sedim surface ice change =',dWtots,' kg m-2' print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2' endif endif 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) * & ( zdqevap(1:ngrid,1:nlayer) & + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) & + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) & + 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, & rnat,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 ! ------------------ ! VI.5. Slab Ocean ! ------------------ if (ok_slab_ocean)then do ig=1,ngrid qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice) enddo !! These two commented functions are from BC2014 (only kept for record) !! !! call ocean_slab_ice(ptimestep, & !! ngrid, knindex, tsea_ice, fluxrad, & !! zdqssnow, qsurf(:,igcm_h2o_ice), & !! - zdqsdif(:,igcm_h2o_vap), & !! flux_sens_lat,tsea_ice, pctsrf_sic, & !! taux,tauy,icount) !! call ocean_slab_get_vars(ngrid,tslab, & !! sea_ice, flux_o, & !! flux_g, dt_hdiff, & !! dt_ekman) call ocean_slab_noice(icount, ptimestep, knon, knindex, & zdqssnow, tsea_ice, & fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, & tsea_ice, taux, tauy, zmasq) call ocean_slab_ice(icount, ptimestep, knon, knindex, & zdqssnow, tsea_ice, & fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, & tsea_ice, -zdqsdif(:,igcm_h2o_vap), taux, tauy, zmasq) call ocean_slab_frac(pctsrf_sic, zmasq) call ocean_slab_get_vars(ngrid, tslab, tice, sea_ice, flux_g, & dt_hdiff, dt_ekman, dt_gm) !!! call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, & !!! dt_hdiff, dt_ekman, dt_gm) !! sea_ice defines the sea ice thickness in kg/m2 !! pctsrf_sic defines the percentage of the oceanic grid that is covered by sea ice (between 0 and **almost** 1) !! qsurf(:,igcm_h2o_ice) defines the amount of snow that accumulates on top of the sea ice (in kg/m2) ; note that the snow fraction depends on the amount of snow. !! therefore, beware that qsurf(:,igcm_h2o_ice) has its own meaning for oceanic grid points (when ok_slab_ocean activated) do ig=1,ngrid if (nint(rnat(ig)).eq.1)then tslab(ig,1) = 0. tslab(ig,2) = 0. tsea_ice(ig) = 0. ! tice(ig) = 0. sea_ice(ig) = 0. pctsrf_sic(ig) = 0. qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice) endif enddo endif ! end of 'ok_slab_ocean' ! ----------------------------- ! VI.6. Surface Tracer Update ! ----------------------------- if(hydrology)then call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & capcal,albedo,albedo_bareground, & albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & sea_ice) zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid) dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq) qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq) ! Test energy conservation if(enertest)then call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots) if (is_master) print*,'In hydrol surface energy change =',dEtots,' W m-2' endif ! test water conservation if(watertest)then call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*cell_area(:)*ptimestep/totarea_planet,dWtots) if (is_master) print*,'In hydrol surface ice change =',dWtots,' kg m-2' call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots) if (is_master) then print*,'In hydrol surface water change =',dWtots,' kg m-2' print*,'---------------------------------------------------------------' endif endif else ! of if (hydrology) qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq) end if ! of if (hydrology) ! 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(:,:) endif! end of if 'tracer' !------------------------------------------------ ! VII. Surface and sub-surface soil temperature !------------------------------------------------ ! Increment surface temperature if(ok_slab_ocean)then do ig=1,ngrid if (nint(rnat(ig)).eq.0)then zdtsurf(ig)= (tslab(ig,1) & + pctsrf_sic(ig)*(tsea_ice(ig)-tslab(ig,1))-tsurf(ig))/ptimestep endif tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) enddo else tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) endif ! end of 'ok_slab_ocean' ! Compute soil temperatures and subsurface heat flux. if (callsoil) then call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, & ptimestep,tsurf,tsoil,capcal,fluxgrd) endif if (ok_slab_ocean) then do ig=1,ngrid fluxgrdocean(ig)=fluxgrd(ig) if (nint(rnat(ig)).eq.0) then capcal(ig)=capcalocean fluxgrd(ig)=0. !! fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1)) ! Dividing by cell area to have flux in W/m2 fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1)+dt_gm(ig,1))/cell_area(ig) ! fluxgrdocean(ig)=(flux_g(ig)-fluxabs_sw(ig))+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1)+dt_gm(ig,1))/cell_area(ig) do iq=1,nsoilmx tsoil(ig,iq)=tsurf(ig) enddo if (pctsrf_sic(ig).gt.0.5) then capcal(ig)=capcalseaice if (qsurf(ig,igcm_h2o_ice).gt.0.) then capcal(ig)=capcalsno endif endif endif enddo endif !end of 'ok_slab_ocean' ! 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 !--------------------------------------------------- ! 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 IF (calltherm) THEN DO ig=1,ngrid DO l=1,nlayer zw2_bis(ig,l) = zw2(ig,l) fm_bis(ig,l) = fm(ig,l) ENDDO ENDDO ENDIF ! 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 ! 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 ! Generalised for arbitrary aerosols now. By LK reffcol(1:ngrid,1:naerkind)=0.0 if(co2cond.and.(iaero_co2.ne.0))then call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2)) do ig=1,ngrid reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)*reffrad(ig,1:nlayer,iaero_co2)*mass(ig,1:nlayer)) enddo endif if(water.and.(iaero_h2o.ne.0))then call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, & reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o)) do ig=1,ngrid reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)*reffrad(ig,1:nlayer,iaero_h2o)*mass(ig,1:nlayer)) enddo endif endif ! end of 'tracer' ! Test for water conservation. if(water)then call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf) call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf) call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol) call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol) h2otot = icesrf + liqsrf + icecol + vapcol if (is_master) then print*,' Total water amount [kg m^-2]: ',h2otot print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] ' print*, icesrf,liqsrf,icecol,vapcol endif if(meanOLR .and. is_master)then if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then ! to record global water balance open(98,file="h2o_bal.out",form='formatted',position='append') write(98,*) zday,icesrf,liqsrf,icecol,vapcol close(98) endif endif endif ! Calculate RH (Relative Humidity) for diagnostic. if(water)then do l = 1, nlayer do ig=1,ngrid call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l)) RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l) enddo enddo ! Compute maximum possible H2O column amount (100% saturation). do ig=1,ngrid H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:)) enddo endif ! end of 'water' ! 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_vap,igcm_generic_ice,call_ice_vap_generic) if (call_ice_vap_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_vap) / qsat_generic(ig,l,iq) enddo enddo end if end do ! iq=1,nq endif ! end of 'generic_condensation' if (is_master) print*,'--> Ls =',zls*180./pi ! ------------------------- ! IX. Thermosphere ! ------------------------- if (callthermos) then call conduction(ngrid,nlayer,nq,ptimestep,pplay,pplev, & pt,tsurf,zzlev,zzlay,muvar,pq,firstcall,zdtconduc) pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtconduc(1:ngrid,1:nlayer) endif ! of if (callthermos) !---------------------------------------------------------------------- ! 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 writting 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 if(lastcall) then ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) #ifndef MESOSCALE if (ngrid.ne.1) then write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin if (ok_slab_ocean) then ! fetch "ocean variables" to ensure they are stored call ocean_slab_get_vars(ngrid, tslab, tice, sea_ice, flux_g, & dt_hdiff, dt_ekman, dt_gm) endif call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin, & tsurf,tsoil,emis,albedo,q2,qsurf_hist, & cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab,tsea_ice,tice,sea_ice) endif #endif if(ok_slab_ocean) then call ocean_slab_final!(tslab, seaice) end if 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,"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 if (water) then vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap) call wstats(ngrid,"vmr_h2ovapor", & "H2O vapour volume mixing ratio","mol/mol", & 3,vmr) endif endif ! end of 'tracer' if(watercond.and.CLFvarying)then call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man) call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc) call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac) call wstats(ngrid,"RH","relative humidity"," ",3,RH) endif if (ok_slab_ocean) then call wstats(ngrid,"dt_hdiff1","dt_hdiff1","W m-2",2,dt_hdiff(:,1)) call wstats(ngrid,"dt_hdiff2","dt_hdiff2","W m-2",2,dt_hdiff(:,2)) call wstats(ngrid,"dt_ekman1","dt_ekman1","W m-2",2,dt_ekman(:,1)) call wstats(ngrid,"dt_ekman2","dt_ekman2","W m-2",2,dt_ekman(:,2)) call wstats(ngrid,"dt_gm1","dt_gm1","W m-2",2,dt_gm(:,1)) call wstats(ngrid,"dt_gm2","dt_gm2","W m-2",2,dt_gm(:,2)) call wstats(ngrid,"tslab1","tslab1","K",2,tslab(:,1)) call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2)) call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic) call wstats(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice) ! call wstats(ngrid,"tice","sea ice temperature","K",2,tice) call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice) call wstats(ngrid,"rnat","nature of the surface","",2,rnat) endif 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 write_output("Ls","solar longitude","deg",zls*180./pi) call write_output("Lss","sub solar longitude","deg",zlss*180./pi) call write_output("RA","right ascension","deg",right_ascen*180./pi) call write_output("Declin","solar declination","deg",declin*180./pi) call write_output("dist_star","dist_star","AU",dist_star) call write_output("tsurf","Surface temperature","K",tsurf) call write_output("ps","Surface pressure","Pa",ps) call write_output("emis","Emissivity","",emis) call write_output("temperature","temperature","K",zt) call write_output("teta","potential temperature","K",zh) call write_output("u","Zonal wind","m.s-1",zu) call write_output("v","Meridional wind","m.s-1",zv) call write_output("w","Vertical wind","m.s-1",pw) call write_output("p","Pressure","Pa",pplay) call write_output("omega","omega","Pa/s",omega) ! Subsurface temperatures ! call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf) ! Oceanic layers if(ok_slab_ocean) then call write_output("pctsrf_sic","pct ice/sea","",pctsrf_sic) call write_output("tsea_ice","top layer temp, snow/ice","K",tsea_ice) ! call write_output("tice","sea ice temperature","K",tice) call write_output("sea_ice","sea ice","kg/m2",sea_ice) call write_output("tslab1","tslab1","K",tslab(:,1)) call write_output("tslab2","tslab2","K",tslab(:,2)) call write_output("dt_hdiff1","dt_hdiff1","W m-2",dt_hdiff(:,1)) call write_output("dt_hdiff2","dt_hdiff2","W m-2",dt_hdiff(:,2)) call write_output("dt_ekman1","dt_ekman1","W m-2",dt_ekman(:,1)) call write_output("dt_ekman2","dt_ekman2","W m-2",dt_ekman(:,2)) call write_output("dt_gm1","dt_gm1","W m-2",dt_gm(:,1)) call write_output("dt_gm2","dt_gm2","W m-2",dt_gm(:,2)) call write_output("rnat","nature of the surface","",rnat) call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux) call write_output("latentFlux","latent heat flux","w.m^-2",zdqsdif(:,igcm_h2o_vap)*RLVTT) endif ! Thermal plume model if (calltherm) then call write_output('entr','Entrainment','kg m$^{-2}$ s$^{-1}$', entr) call write_output('detr','Detrainment','kg m$^{-2}$ s$^{-1}$', detr) call write_output('fm','Mass flux','kg m$^{-2}$ s$^{-1}$', fm_bis) call write_output('w_plm','Squared vertical velocity','m s$^{-1}$', zw2_bis) call write_output('fraca','Updraft fraction','', fraca) call write_output("zdttherm","dt due to plumes","K-1 s-1",zdttherm) !ALS24 IF (tracer) THEN call write_output("zdqtherm","dq due to plumes, iq = 1","kg kg-1 s-1",& zdqtherm(1:ngrid,1:nlayer,1)) !ALS24 call write_output("zdqtherm2","dq due to plumes; iq = 2","kg kg-1 s-1",& zdqtherm(1:ngrid,1:nlayer,2)) !ALS24 ENDIF endif ! GW non-oro outputs if (calllott_nonoro) then call write_output("dugwno","GW non-oro dU","m/s2", d_u_hin) call write_output("dvgwno","GW non-oro dV","m/s2", d_v_hin) endif ! Total energy balance diagnostics if(callrad.and.(.not.newtonian))then call write_output("ALB","Surface albedo"," ",albedo_equivalent) call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw) call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn) call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw) call write_output("DYN","dynamical heat input","W m-2",fluxdyn) call write_output("shad","rings"," ", fract) ! call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1) ! call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1) ! 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,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1) ! call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1) ! call writediagfi(ngrid,"netfluxsurflw","lw net surface flux","W m-2",2,net_fluxsurf_lw) if(ok_slab_ocean) then call write_output("GND","heat flux from ground","W m-2",fluxgrdocean) else call write_output("GND","heat flux from ground","W m-2",fluxgrd) endif endif ! end of 'callrad' if(enertest) then if (calldifv) then call write_output("q2","turbulent kinetic energy","J.kg^-1",q2) call write_output("sensibFlux","sensible heat flux","w.m^-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 write_output("dEzradsw","radiative heating","w.m^-2",dEzradsw) call write_output("dEzradlw","radiative heating","w.m^-2",dEzradlw) endif if(watercond) then call write_output("lscaledE","heat from largescale","W m-2",lscaledE) if ((.not.calltherm).and.moistadjustment) then call write_output("madjdE","heat from moistadj","W m-2",madjdE) endif call write_output("qsatatm","atm qsat"," ",qsat) ! call write_output("lscaledEz","heat from largescale","W m-2",lscaledEz) ! call write_output("madjdEz","heat from moistadj","W m-2",madjdEz) ! call write_output("h2o_max_col","maximum H2O column amount","kg.m^-2",H2Omaxcol) endif if (generic_condensation) then call write_output("genericconddE","heat from generic condensation","W m-2",genericconddE) call write_output("dt_generic_condensation","heating from generic condensation","K s-1",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 write_output('dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',int_dtauv(:,nlayer:1:-1,nw)) enddo do nw=1,L_NSPECTI write(str2,'(i2.2)') nw call write_output('dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',int_dtaui(:,nlayer:1:-1,nw)) enddo endif ! Temporary inclusions for heating diagnostics. call write_output("zdtsw","SW heating","T s-1",zdtsw) call write_output("zdtlw","LW heating","T s-1",zdtlw) call write_output("dtrad","radiative heating","K s-1",dtrad) call write_output("zdtdyn","Dyn. heating","T s-1",zdtdyn) if (tracer .or. water) call write_output("dtmoistadj","moist adj heating","K s-1",dtmoist) if (calladj) call write_output("dtdryadj","dry adj heating","K s-1",zdtadj) ! For Debugging. !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat)) !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi) ! Output aerosols. if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & call write_output('CO2ice_reff','CO2ice_reff','m',reffrad(1,1,iaero_co2)) if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & call write_output('H2Oice_reff','H2Oice_reff','m',reffrad(:,:,iaero_h2o)) if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & call write_output('CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',reffcol(1,iaero_co2)) if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & call write_output('H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',reffcol(1,iaero_h2o)) ! Output tracers. if (tracer) then do iq=1,nq call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq)) call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & 'kg m^-2',qsurf_hist(1,iq) ) call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col', & 'kg m^-2',qcol(:,iq) ) ! call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & ! 'kg m^-2',qsurf(:,iq) ) if(watercond.or.CLFvarying)then call write_output("rneb_man","H2O cloud fraction (conv)"," ",rneb_man) call write_output("rneb_lsc","H2O cloud fraction (large scale)"," ",rneb_lsc) call write_output("CLF","H2O cloud fraction"," ",cloudfrac) call write_output("CLFt","H2O column cloud fraction"," ",totcloudfrac) call write_output("RH","relative humidity"," ",RH) call write_output("vteta","virtual potential temperature","K",zh * (1.+(1./epsi-1.) * & zq(1:ngrid,1:nlayer,1))) endif if(waterrain)then call write_output("rain","rainfall","kg m-2 s-1",zdqsrain) call write_output("snow","snowfall","kg m-2 s-1",zdqssnow) call write_output("reevap","reevaporation of precipitation","kg m-2 s-1",reevap_precip) endif if(generic_condensation)then call write_output("CLF","GCS cloud fraction"," ",cloudfrac) !AF24: TODO fix rneb_generic and RH_generic failure with write_output() ! call write_output("rneb_generic","GCS cloud fraction (generic condensation)"," ",rneb_generic) ! call write_output("RH_generic","GCS relative humidity"," ",RH_generic) call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic) call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic) call write_output("vteta","virtual potential temperature","K",zh * (1.+(1./epsi_generic-1.) * & zq(1:ngrid,1:nlayer,1) )) endif if(generic_rain)then call write_output("rain","generic rainfall","kg m-2 s-1",zdqsrain_generic) call write_output("snow","generic snowfall","kg m-2 s-1",zdqssnow_generic) call write_output("reevap","generic reevaporation of precipitation","kg m-2 s-1",reevap_precip_generic) endif if((hydrology).and.(.not.ok_slab_ocean))then call write_output("hice","oceanic ice height","m",hice) endif call write_output("tau_col","Total aerosol optical depth","[]",tau_col) enddo ! end of 'nq' loop 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.(water)) then comm_CLOUDFRAC(1:ngrid,1:nlayer)=cloudfrac(1:ngrid,1:nlayer) comm_TOTCLOUDFRAC(1:ngrid)=totcloudfrac(1:ngrid) comm_SURFRAIN(1:ngrid)=zdqsrain(1:ngrid) comm_DQVAP(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_vap) comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_ice) comm_H2OICE_REFF(1:ngrid,1:nlayer)=reffrad(1:ngrid,1:nlayer,iaero_h2o) comm_REEVAP(1:ngrid)=reevap_precip(1:ngrid) comm_DTRAIN(1:ngrid,1:nlayer)=zdtrain(1:ngrid,1:nlayer) comm_DTLSC(1:ngrid,1:nlayer)=dtlscale(1:ngrid,1:nlayer) comm_RH(1:ngrid,1:nlayer)=RH(1:ngrid,1:nlayer) else if ((tracer).and.(generic_condensation).and.(.not. water)) 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_vap,igcm_generic_ice,call_ice_vap_generic) if (call_ice_vap_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_vap) 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) ! if (turb_resolved) then ! open(17,file='lsf.txt',form='formatted',status='old') ! rewind(17) ! DO l=1,nlayer ! read(17,*) lsf_dt(l),lsf_dq(l) ! ENDDO ! close(17) ! do ig=1,ngrid ! if ((tracer).and.(water)) then ! pdq(ig,:,igcm_h2o_vap) = pdq(ig,:,igcm_h2o_vap) + lsf_dq(:) ! endif ! pdt(ig,:) = pdt(ig,:) + lsf_dt(:) ! comm_HR_DYN(ig,:) = lsf_dt(:) ! enddo ! endif #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("controle",tab_cntrl_mod,1) CALL send_xios_field("ap",ap,1) CALL send_xios_field("bp",bp,1) CALL send_xios_field("aps",aps,1) CALL send_xios_field("bps",bps,1) IF (water) THEN CALL send_xios_field('h2o_layer1',zq(:,1,igcm_h2o_vap)) CALL send_xios_field('co2_layer1',zq(:,1,igcm_co2_ice)) ENDIF 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