Index: LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 (revision 5197) +++ LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 (revision 5198) @@ -1,53 +1,11 @@ +! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $ MODULE physiqex_mod + IMPLICIT NONE + CONTAINS -! ------------------------------------------------------------------------------------------------ -! -! Main authors : F. Hourdin, S. Riette, Q. Libois, A. Idelkadi -! -------------- -! -! Object : A version of the LMDZ physics that include PHYEX -! ------- Tu be used in particular to be coupled to non hydrostatic version of the dynamical -! core dynamico or its limited area version. -! -! -! Development steps : -! ------------------- -! -! physiqex_mod.F90 has been introducued in phylmd as an empty driver for the physics (with only -! a newtonian relaxation for temperature and drag at the surface). -! This driver has been used to introduce the PHYEX package, an externailzed physics package coming -! from MesoNH, representing turbulence, convection, clouds and precipitations. -! This work has been done for a large part in Dephy sessions (Oleron, 2022, Fréjus 2023, Lège 2024). -! -! Fréjus 2023 : first version runing in 1D for Arm cumulus -! Lège 2024 : first versions with a generic soil model (developped at this occasion) and ECrad -! -! -! To do list -! --------- -! To do list pour un branchement plus propre -! * PHYEX considère toutes les espèces microphysiques et la tke comme pronostiques, des termes de tendances -! sont calculés. L'avance temporelle est faite en fin de pas de temps mais pourrait être déplacée au dessus si -! les tendances de ces variables passaient par l'interface -! * PHYEX a besoin de variables avec mémoire d'un pas de temps à l'autre (en plus des variables pronostiques -! décrites juste au dessus). Ce sont les tableaux ALLOCATABLE. Ces variables pourraient devenir -! des traceurs. -! * La variable ZDZMIN est ici calculée avec un MINVAL qui conduira à des résultats différents -! lorsque la distribution (noeuds/procs) sera différente. Cette variable est utile pour déterminer -! un critère CFL pour la sédimentation des précipitations. Il suffit donc d'avoir une valeur -! approchante. -! * Certains allocatable sont en klev, d'autres en klev+2. Il est possible de changer ceci mais il faut gérer -! les recopies pour que les params voient des tableaux en klev+2 (en effectuant des recopies des -! niveaux extrêmes comme fait pour le vent par exemple) -! * L'eau tombée en surface (précipitations, sédimentation du nuage, terme de dépôt) se trouve dans -! les variables ZINPRC, ZINPRR, ZINPRS, ZINPRG -! -!------------------------------------------------------------------------------------------------- - - SUBROUTINE physiqex (nlon,nlev, & - & debut,lafin,pdtphys_, & + & debut,lafin,pdtphys, & & paprs,pplay,pphi,pphis,presnivs, & & u,v,rot,t,qx, & @@ -57,40 +15,10 @@ USE dimphy, only : klon,klev USE infotrac_phy, only : nqtot - USE geometry_mod, only : latitude, cell_area, latitude_deg, longitude_deg + USE geometry_mod, only : latitude ! USE comcstphy, only : rg USE ioipsl, only : ymds2ju - !USE phys_state_var_mod, only : phys_state_var_init - USE phys_state_var_mod + USE phys_state_var_mod, only : phys_state_var_init USE phyetat0_mod, only: phyetat0 USE output_physiqex_mod, ONLY: output_physiqex -!!!!!! AI 04-2024 - USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time - USE limit_read_mod, ONLY : init_limit_read - USE conf_phys_m, only: conf_phys - USE ozonecm_m, only: ozonecm - USE aero_mod - USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, & - year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour, calend - USE radlwsw_m, only: radlwsw - USE MODD_CST - USE phys_output_var_mod, ONLY : phys_output_var_init, cloud_cover_sw - USE write_field_phy - USE phyaqua_mod, only: zenang_an -!!!!!! - ! PHYEX internal modules - USE MODE_INIT_PHYEX, ONLY: INIT_PHYEX, FILL_DIMPHYEX - USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t - USE MODD_PHYEX, ONLY: PHYEX_t - USE MODI_INI_PHYEX, ONLY: INI_PHYEX - USE MODI_ICE_ADJUST, ONLY: ICE_ADJUST - USE MODI_RAIN_ICE, ONLY: RAIN_ICE - USE MODI_TURB - USE MODI_SHALLOW_MF - - USE MODD_CST - USE ioipsl_getin_p_mod, ONLY : getin_p - USE generic_soil_ini, ONLY : soil_ini - USE generic_soil_conduct, ONLY : soil_conduct - USE generic_soil_fluxes, ONLY : soil_fluxes IMPLICIT none @@ -98,9 +26,10 @@ ! Routine argument: ! + integer,intent(in) :: nlon ! number of atmospheric colums integer,intent(in) :: nlev ! number of vertical levels (should be =klev) logical,intent(in) :: debut ! signals first call to physics logical,intent(in) :: lafin ! signals last call to physics - real,intent(in) :: pdtphys_ ! physics time step (s) + real,intent(in) :: pdtphys ! physics time step (s) real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) @@ -119,157 +48,6 @@ real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure - real :: d_qr(klon, klev), d_qs(klon, klev), d_qg(klon, klev) ! tendency for rain, snow, graupel - real :: d_tke(klon, klev) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! AI 04-2024 a nettoyer - ! - CHARACTER (LEN=20) :: modname='physiq_mod' - CHARACTER*80 abort_message - ! Declarations pas de temps - INTEGER, SAVE :: itap=0 ! compteur pour la physique - !$OMP THREADPRIVATE(itap) - INTEGER lmt_pas - SAVE lmt_pas ! frequence de mise a jour - !$OMP THREADPRIVATE(lmt_pas) - REAL zdtime, zdtime1, zdtime2, zlongi - ! Declarations pour appel de conf_phys - LOGICAL ok_hf - SAVE ok_hf - !$OMP THREADPRIVATE(ok_hf) - LOGICAL ok_journe ! sortir le fichier journalier - SAVE ok_journe - !$OMP THREADPRIVATE(ok_journe) - LOGICAL ok_mensuel ! sortir le fichier mensuel - SAVE ok_mensuel - !$OMP THREADPRIVATE(ok_mensuel) - LOGICAL ok_instan ! sortir le fichier instantane - SAVE ok_instan - !$OMP THREADPRIVATE(ok_instan) - LOGICAL ok_LES ! sortir le fichier LES - SAVE ok_LES - !$OMP THREADPRIVATE(ok_LES) - LOGICAL callstats ! sortir le fichier stats - SAVE callstats - !$OMP THREADPRIVATE(callstats) - LOGICAL ok_region ! sortir le fichier regional - PARAMETER (ok_region=.FALSE.) - REAL seuil_inversion - SAVE seuil_inversion - !$OMP THREADPRIVATE(seuil_inversion) - ! Parametres lies au nouveau schema de nuages (SB, PDF) - REAL, SAVE :: fact_cldcon - REAL, SAVE :: facttemps - !$OMP THREADPRIVATE(fact_cldcon,facttemps) - LOGICAL, SAVE :: ok_newmicro - !$OMP THREADPRIVATE(ok_newmicro) - INTEGER, SAVE :: iflag_cld_th - !$OMP THREADPRIVATE(iflag_cld_th) - REAL ratqsbas,ratqshaut,tau_ratqs - SAVE ratqsbas,ratqshaut,tau_ratqs - !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs) - ! ozone - INTEGER,SAVE :: read_climoz ! Read ozone climatology - ! (let it keep the default OpenMP shared attribute) - ! Allowed values are 0, 1 and 2 - ! 0: do not read an ozone climatology - ! 1: read a single ozone climatology that will be used day and night - ! 2: read two ozone climatologies, the average day and night - ! climatology and the daylight climatology - INTEGER,SAVE :: ncid_climoz - REAL zzzo - REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for - ! the ozone fields, old method. - ! - ! aerosols - ! params et flags aerosols - LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not - LOGICAL ok_alw ! Apply aerosol LW effect or not - LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013) - REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995) - SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1 - !$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1) - LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA - ! false : lecture des aerosol dans un fichier - INTEGER, SAVE :: flag_aerosol - !$OMP THREADPRIVATE(flag_aerosol) - LOGICAL, SAVE :: flag_bc_internal_mixture - !$OMP THREADPRIVATE(flag_bc_internal_mixture) - ! - !--STRAT AEROSOL - INTEGER, SAVE :: flag_aerosol_strat - !$OMP THREADPRIVATE(flag_aerosol_strat) - ! - !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION - LOGICAL, SAVE :: flag_aer_feedback - !$OMP THREADPRIVATE(flag_aer_feedback) - ! - LOGICAL, SAVE :: chemistry_couple ! true : use INCA chemistry O3 - ! false : use offline chemistry O3 - !$OMP THREADPRIVATE(chemistry_couple) - LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques - !$OMP THREADPRIVATE(ok_volcan) - INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato - !$OMP THREADPRIVATE(flag_volc_surfstrat) - ! - INTEGER, SAVE :: iflag_radia=0 ! active ou non le rayonnement (MPL) - !$OMP THREADPRIVATE(iflag_radia) - ! - REAL, SAVE :: alp_offset - !$OMP THREADPRIVATE(alp_offset) -! Declarations pour le rayonnement -! Le rayonnement n'est pas calcule tous les pas, il faut donc -! sauvegarder les sorties du rayonnement - ! Surface - REAL zxtsol(KLON) - INTEGER itaprad - SAVE itaprad - !$OMP THREADPRIVATE(itaprad) - ! Date equinoxe de printemps - !LOGICAL, parameter :: ok_rad=.true. ! a nettoyer - REAL,SAVE :: solarlong0 - !$OMP THREADPRIVATE(solarlong0) - ! 6 bandes SW : argument non utilise dans radlwsw a nettoyer ? - REAL,DIMENSION(6), parameter :: SFRWL = (/ 1.28432794E-03, 0.12304168, 0.33106142, & - & 0.32870591, 0.18568763, 3.02191470E-02 /) - LOGICAL, parameter :: new_orbit = .TRUE. - INTEGER, parameter :: mth_eq=3, day_eq=21 - REAL :: day_since_equinox, jD_eq - !REAL zdtime, zdtime1, zdtime2, zlongi - REAL dist, rmu0(klon), fract(klon), zrmu0(klon), zfract(klon) - CHARACTER(len=512) :: namelist_ecrad_file - ! - REAL, dimension(klon, klev) :: zqsat - REAL, dimension(klon, klev) :: ref_liq, ref_ice, ref_liq_pi, ref_ice_pi ! rayons effectifs des gout - REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique - REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique - REAL, dimension(klon, klev) :: cldemirad ! emissivite pour - REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse - ! aerosols AI attention a les declarer ailleurs dans phylmd : phys_output ou phy_stats et non phys_local - REAL, dimension(klon,klev+1) :: ZLWFT0_i, ZSWFT0_i,ZFLDN0, ZFLUP0, ZFSDN0, ZFSUP0 - REAL, dimension(klon) :: oplwad_aero, sollwad_aero,& - toplwai_aero, sollwai_aero, & - toplwad0_aero, sollwad0_aero - REAL, dimension(klon,9) :: solsw_aero, solsw0_aero,topsw_aero, topsw0_aero - - REAL, dimension(klon) :: & - topswcf_aero, solswcf_aero, & - solswad_aero, & - solswad0_aero, & - topswai_aero, solswai_aero, & - toplwad_aero, topswad0_aero, & - topswad_aero - REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer - - REAL,dimension(klon,klev) :: d_t_swr, d_t_sw0, d_t_lwr, d_t_lw0 - ! - EXTERNAL suphel ! AI attention - ! - include "YOMCST.h" - include "clesphys.h" -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Fin AI - - include "flux_arp.h" - +! include "clesphys.h" INTEGER length PARAMETER ( length = 100 ) @@ -279,143 +57,16 @@ !$OMP THREADPRIVATE(clesphy0) -! Saved variables -REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PSIGS !variance of s -REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PCF_MF, PRC_MF, PRI_MF !shallow convection cloud -REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ZQR, ZQS, ZQG !rain, snow, graupel specifiq contents -REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PTKEM ! TKE -TYPE(DIMPHYEX_t),SAVE :: D -TYPE(PHYEX_t), SAVE :: PHYEX -! -INTEGER, PARAMETER :: KRR=6, KRRL=2, KRRI=3, KSV=0 -INTEGER :: JRR -! Time-State variables and Sources variables -REAL, DIMENSION(klon,klev+2) :: ZUT, ZVT ! wind component on klev+2 -REAL, DIMENSION(klon,klev+2) :: ZPABST ! absolute pressure -REAL, DIMENSION(klon,klev+2,KRR) :: ZRX,ZRXS,ZRXS0 ! qx and source of q from LMDZ to rt -REAL, DIMENSION(klon,klev+2) :: ZRHOD ! rho_dry -REAL, DIMENSION(klon,klev+2,0) :: ZSVT ! passive scal. var. -REAL, DIMENSION(klon,klev+2) :: PWT ! vertical wind velocity (used only for diagnostic) -REAL, DIMENSION(klon,klev+2) :: ZRUS,ZRVS,ZRWS,ZRTHS,ZRTKES ! sources of momentum, conservative potential temperature, Turb. Kin. Energy -REAL, DIMENSION(KLON,KLEV+2) :: ZTHETAS !tendency -REAL, DIMENSION(KLON,KLEV+2) :: ZTHETAS0 -REAL, DIMENSION(KLON,KLEV+2) :: ZTKES -REAL, DIMENSION(KLON,KLEV+2) :: ZTKES0 -! Source terms for all water kinds, PRRS(:,:,:,1) is used for the conservative ! mixing ratio -REAL, DIMENSION(klon,klev+2,KRR) :: ZRRS -REAL, DIMENSION(klon,klev+2) :: PTHVREF ! Virtual Potential Temperature of the reference state -! Adjustment variables -REAL, DIMENSION(KLON,KLEV+2) :: zqdm !1-qt=1/(1+rt) -REAL, DIMENSION(KLON,KLEV+2) :: zqt !mixing ratio and total specifiq content -REAL, DIMENSION(KLON,KLEV+2) :: ZTHETA, ZEXN !theta and exner function -REAL, DIMENSION(KLON,KLEV+2) :: zz_mass !altitude above ground of mass points -REAL, DIMENSION(KLON,KLEV+2) :: zz_flux !altitude above ground of flux points -REAL, DIMENSION(KLON,KLEV+2) :: zdzm !distance between two consecutive mass points (expressed on flux points) -REAL, DIMENSION(KLON,KLEV+2) :: zdzf !distance between two consecutive flux points (expressed on mass points) -REAL, DIMENSION(KLON) :: zs !surface orography -REAL, DIMENSION(KLON) :: zsigqsat !coeff for the extra term for s variance -REAL, DIMENSION(0) :: ZMFCONV -REAL, DIMENSION(KLON,KLEV+2) :: ZICLDFR, ZWCLDFR, ZSSIO, ZSSIU, ZIFR, ZICE_CLD_WGT !used only in HIRLAM config -REAL, DIMENSION(KLON,KLEV+2) :: ZSRC ! bar(s'rc') -REAL, DIMENSION(KLON,KLEV+2) :: ZCLDFR !cloud fraction -REAL, DIMENSION(KLON,KLEV+2) :: ZHLC_HRC, ZHLC_HCF, ZHLI_HRI, ZHLI_HCF !subgrid autoconversion -! Rain_ice variables -real :: ZDZMIN -real, dimension(klon, klev+2) :: ZCIT !ice concentration -real, dimension(klon) :: ZINPRC, ZINPRR, ZINPRS, ZINPRG !precipitation flux at ground -real, dimension(klon, klev+2) :: ZEVAP3D !evaporation (diag) -real, dimension(klon, klev+2) :: ZRAINFR -real, dimension(klon) :: ZINDEP !deposition flux, already contained in ZINPRC -real, dimension(klon) :: ZSEA, ZTOWN !sea and town fractions in the frid cell -! Turbulence variables -REAL, DIMENSION(klon,klev+2) :: PDXX,PDYY,PDZX,PDZY ! metric coefficients -REAL, DIMENSION(klon,klev+2) :: PZZ ! physical distance between 2 succesive grid points along the K direction -REAL, DIMENSION(klon) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW ! Director Cosinus along x, y and z directions at surface w-point -REAL, DIMENSION(klon) :: PCOSSLOPE ! cosinus of the angle between i and the slope vector -REAL, DIMENSION(klon) :: PSINSLOPE ! sinus of the angle between i and the slope vector -REAL, DIMENSION(klon,klev+2) :: PRHODJ ! dry density * Grid size -REAL, DIMENSION(0,0) :: MFMOIST ! moist mass flux dual scheme -REAL, DIMENSION(0,0,0) :: PHGRAD ! horizontal gradients -REAL, DIMENSION(klon) :: PSFTH,PSFRV,PSFU,PSFV ! normal surface fluxes of theta, Rv, (u,v) parallel to the orography -REAL, DIMENSION(klon,0) :: PSFSV ! normal surface fluxes of Scalar var. KSV=0 -REAL, DIMENSION(klon) :: ZBL_DEPTH ! BL height for TOMS -REAL, DIMENSION(klon) :: ZSBL_DEPTH ! SBL depth for RMC01 -REAL, DIMENSION(klon,klev+2) :: ZCEI ! Cloud Entrainment instability index to emphasize localy turbulent fluxes -REAL, DIMENSION(klon,klev+2,0):: ZRSVS ! Source terms for all passive scalar variables -REAL, DIMENSION(klon,klev+2) :: ZFLXZTHVMF ! MF contribution for vert. turb. transport used in the buoy. prod. of TKE -REAL, DIMENSION(klon,klev+2) :: ZWTH ! heat flux -REAL, DIMENSION(klon,klev+2) :: ZWRC ! cloud water flux -REAL, DIMENSION(klon,klev+2,0):: ZWSV ! scalar flux -REAL, DIMENSION(klon,klev+2) :: ZTP ! Thermal TKE production (MassFlux + turb) -REAL, DIMENSION(klon,klev+2) :: ZDP ! Dynamic TKE production -REAL, DIMENSION(klon,klev+2) :: ZTDIFF ! Diffusion TKE term -REAL, DIMENSION(klon,klev+2) :: ZTDISS ! Dissipation TKE term -REAL, DIMENSION(0,0) :: PLENGTHM, PLENGTHH ! length scale from vdfexcu (HARMONIE-AROME) -REAL :: ZTHVREFZIKB ! for electricity scheme interface -REAL, DIMENSION(klon,klev+2) :: ZDXX,ZDYY,ZDZX,ZDZY,ZZZ -REAL, DIMENSION(klon) :: ZDIRCOSXW,ZDIRCOSYW,ZDIRCOSZW,ZCOSSLOPE,ZSINSLOPE -! Shallow variables -REAL, DIMENSION(klon,klev+2) :: PDUDT_MF ! tendency of U by massflux scheme -REAL, DIMENSION(klon,klev+2) :: PDVDT_MF ! tendency of V by massflux scheme -REAL, DIMENSION(klon,klev+2) :: PDTHLDT_MF ! tendency of thl by massflux scheme -REAL, DIMENSION(klon,klev+2) :: PDRTDT_MF ! tendency of rt by massflux scheme -REAL, DIMENSION(klon,klev+2,KSV):: PDSVDT_MF ! tendency of Sv by massflux scheme -REAL, DIMENSION(klon,klev+2) :: PSIGMF -REAL, DIMENSION(klon,klev+2) :: ZFLXZTHMF -REAL, DIMENSION(klon,klev+2) :: ZFLXZRMF -REAL, DIMENSION(klon,klev+2) :: ZFLXZUMF -REAL, DIMENSION(klon,klev+2) :: ZFLXZVMF -REAL, DIMENSION(klon,klev+2) :: PTHL_UP ! Thl updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PRT_UP ! Rt updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PRV_UP ! Vapor updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PU_UP ! U wind updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PV_UP ! V wind updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PRC_UP ! cloud content updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PRI_UP ! ice content updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PTHV_UP ! Thv updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PW_UP ! vertical speed updraft characteristics -REAL, DIMENSION(klon,klev+2) :: PFRAC_UP ! updraft fraction -REAL, DIMENSION(klon,klev+2) :: PEMF ! updraft mass flux -REAL, DIMENSION(klon,klev+2) :: PDETR ! updraft detrainment -REAL, DIMENSION(klon,klev+2) :: PENTR ! updraft entrainment -INTEGER,DIMENSION(klon) ::IKLCL,IKETL,IKCTL ! level of LCL,ETL and CTL -! Values after saturation adjustement -REAL, DIMENSION(klon,klev) :: t_adj ! Adjusted temperature -REAL, DIMENSION(klon,klev) :: qv_adj, ql_adj, qi_adj, qr_adj, qs_adj, qg_adj !specific contents after adjustement -! + real :: temp_newton(klon,klev) integer :: k logical, save :: first=.true. !$OMP THREADPRIVATE(first) -!real,save :: rg=9.81 -!!$OMP THREADPRIVATE(rg) + +real,save :: rg=9.81 +!$OMP THREADPRIVATE(rg) ! For I/Os integer :: itau0 - -real, dimension(nlon) :: swdnsrf,lwdnsrf,lwupsrf,sensible,latent,qsats,capcal,fluxgrd, fluxu,fluxv,fluxt,fluxq - -integer, parameter :: nsoil=11 -real, dimension(nsoil) :: zout -real, dimension(nlon,nlev) :: tsoil_out - -integer :: isoil -real,dimension(:), allocatable, save :: tsrf -real, dimension(:,:), allocatable, save :: tsoil -!$OMP THREADPRIVATE(tsrf,tsoil) - -real :: thermal_inertia0,z0m0,z0h0,z0q0,betasoil - -!AI -!real, save :: rpi, stefan, karman, fsw0, timesec -!!$OMP THREADPRIVATE(rpi, stefan, karman, fsw0, timesec) -real, save :: stefan, fsw0, timesec -!$OMP THREADPRIVATE(stefan, fsw0, timesec) -integer, save :: iecri=1,iflag_surface=0 -!$OMP THREADPRIVATE(iecri,iflag_surface) -real,save :: zjulian -!$OMP THREADPRIVATE(zjulian) -! AI -!integer, save :: itap=0 -!!$OMP THREADPRIVATE(iecri,iflag_surface,rpi, stefan, karman, fsw0, timesec,tsrf, tsoil,itap) +real :: zjulian @@ -425,124 +76,28 @@ print*,'Debut physiqex',debut -! AI -pdtphys=pdtphys_ -CALL update_time(pdtphys) -phys_tstep=NINT(pdtphys) - ! initializations if (debut) then ! Things to do only for the first call to physics - print*,'Debut physiqex IN' - print*,'NOUVEAU PHYSIQEX 222222222222222222222' - CALL suphel ! initialiser constantes et parametres phys. - ! Read parameters and flags in files .def - CALL conf_phys(ok_journe, ok_mensuel, & - ok_instan, ok_hf, & - ok_LES, & - callstats, & - solarlong0,seuil_inversion, & - fact_cldcon, facttemps,ok_newmicro,iflag_radia, & - iflag_cld_th,ratqsbas,ratqshaut,tau_ratqs, & - ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, & - chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, & - flag_bc_internal_mixture, bl95_b0, bl95_b1, & - read_climoz, alp_offset) +print*,'Debut physiqex IN' - ! IF (.NOT. create_etat0_limit) - CALL init_limit_read(days_elapsed) +! load initial conditions for physics (including the grid) + call phys_state_var_init(1) ! some initializations, required before calling phyetat0 + call phyetat0("startphy.nc", clesphy0, tabcntr0) - ! load initial conditions for physics (including the grid) - ! AI - !call phys_state_var_init(1) ! some initializations, required before calling phyetat0 - call phys_state_var_init(read_climoz) ! some initializations, required before calling phyetat0 - call phyetat0("startphy.nc", clesphy0, tabcntr0) - CALL phys_output_var_init - - ! AI - ! Init increments et pas - !itap = 0 - itaprad = 0 - lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours - WRITE(*,*)'La frequence de lecture surface est de ', & - lmt_pas - IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN - radpas = NINT( 86400./phys_tstep)/nbapp_rad - ELSE - WRITE(*,*) 'le nombre de pas de temps physique doit etre un ', & - 'multiple de nbapp_rad' - WRITE(*,*) 'changer nbapp_rad ou alors commenter ce test ', & - 'mais 1+1<>2' - abort_message='nbre de pas de temps physique n est pas multiple ' & - // 'de nbapp_rad' - CALL abort_physic(modname,abort_message,1) - ENDIF +! Initialize outputs: + itau0=0 + ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 + !CALL ymds2ju(annee0, month, dayref, hour, zjulian) + call ymds2ju(1979, 1, 1, 0.0, zjulian) - ! Initialize outputs: - itau0=0 - ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 - !CALL ymds2ju(annee0, month, dayref, hour, zjulian) - !call ymds2ju(1979, 1, 1, 0.0, zjulian) - !print*,'zjulian=',zjulian - - ZDZMIN=MINVAL((pphi(:,2:) - pphi(:,1:klev-1))/9.81) - CALL INIT_PHYEX(pdtphys, ZDZMIN, PHYEX) - CALL FILL_DIMPHYEX(KLON, KLEV, D) - - ! - ! Variables saved - ALLOCATE(PTKEM(klon,klev+2)) - ALLOCATE(PSIGS(klon,klev+2)) - ALLOCATE(PCF_MF(klon,klev+2)) - ALLOCATE(PRC_MF(klon,klev+2)) - ALLOCATE(PRI_MF(klon,klev+2)) - ALLOCATE(ZQR(klon, klev)) - ALLOCATE(ZQS(klon, klev)) - ALLOCATE(ZQG(klon, klev)) - PSIGS=0. - PCF_MF=0. - PRC_MF=0. - PRI_MF=0. - ZQR=0. - ZQS=0. - ZQG=0. - PTKEM(:,:) = PHYEX%TURBN%XTKEMIN ! TODO: init from TKE at stationnary state - - ! ---------------------------------------------------------------- - ! Initialisation du sol generique - ! ---------------------------------------------------------------- - allocate(tsrf(nlon)) - allocate(tsoil(nlon,nsoil)) - thermal_inertia0=2000. - z0m0=0.05 - z0h0=0.05 - z0q0=0.05 - betasoil=0.2 - call soil_ini(klon,nsoil,PHYEX%CST%XCPD,PHYEX%CST%XLVTT,thermal_inertia0,z0m0,z0h0,z0q0,betasoil) -! #ifndef CPP_IOIPSL_NO_OUTPUT ! Initialize IOIPSL output file #endif -! AI Partie Surface depplacee -! compute tendencies to return to the dynamics: -! "friction" on the first layer -! tsrf(1:nlon)=t(1:nlon,1) -! do isoil=1,nsoil -! tsoil(1:nlon,isoil)=tsrf(1:nlon) -! !tsoil(1:nlon,isoil)=302. -! enddo -! rpi=2.*asin(1.) -! timesec=5*3600 -! stefan=5.67e-8 -! fsw0=900. -! call getin_p('iecri',iecri) -! call getin_p('iflag_surface',iflag_surface) endif ! of if (debut) -! AI -! Increment -!itap = itap + 1 !------------------------------------------------------------ ! Initialisations a chaque pas de temps !------------------------------------------------------------ + ! set all tendencies to zero @@ -551,672 +106,28 @@ d_t(1:klon,1:klev)=0. d_qx(1:klon,1:klev,1:nqtot)=0. -d_qr(1:klon,1:klev)=0. -d_qs(1:klon,1:klev)=0. -d_qg(1:klon,1:klev)=0. d_ps(1:klon)=0. -d_tke(1:klon,1:klev)=0. -! -! AI attention -d_t_swr = 0. -d_t_sw0 = 0. -d_t_lwr = 0. -d_t_lw0 = 0. -! -ZDXX(:,:) = 0. -ZDYY(:,:) = 0. -ZDZX(:,:) = 0. -ZDZY(:,:) = 0. -ZDIRCOSXW(:) = 1. -ZDIRCOSYW(:) = 1. -ZDIRCOSZW(:) = 1. -ZCOSSLOPE(:) = 0. -ZSINSLOPE(:) = 1. -PHGRAD(:,:,:) = 0. -ZBL_DEPTH(:) = 0. ! needed only with LRMC01 key (correction in the surface boundary layer) -ZSBL_DEPTH(:) = 0. -ZCEI(:,:) = 0. ! needed only if HTURBLEN_CL /= 'NONE' modification of mixing lengh inside clouds -ZSVT(:,:,:) = 0. -PWT(:,:) = 0. -ZUT(:,2:klev+1) = u(:,:) -ZVT(:,2:klev+1) = v(:,:) -! -!------------------------------------------------------------ -! Conversions and extra levels -!------------------------------------------------------------ -!TODO check in Meso-NH how values are extrapolated outside of the physical domain -zqt(:,2:klev+1) = qx(:,:,1) + qx(:,:,2) + qx(:,:,3) -zqt(:,1)=0. -zqt(:,klev+2)=0. -zqdm(:,:)=1.-zqt(:,:) !equal to 1/(1+rt) -ZRX(:,2:klev+1,1) = qx(:,:,1) / zqdm(:,2:klev+1) -ZRX(:,2:klev+1,2) = qx(:,:,2) / zqdm(:,2:klev+1) -ZRX(:,2:klev+1,4) = qx(:,:,3) / zqdm(:,2:klev+1) -ZRX(:,2:klev+1,3) = ZQR(:,:) -ZRX(:,2:klev+1,5) = ZQS(:,:) -ZRX(:,2:klev+1,6) = ZQG(:,:) -DO JRR=1,KRR - CALL VERTICAL_EXTEND(ZRX(:,:,JRR),klev) -END DO -! -ZEXN(:,2:klev+1) = (pplay(:,:) / PHYEX%CST%XP00) ** (PHYEX%CST%XRD/PHYEX%CST%XCPD) -ZTHETA(:,2:klev+1) = t(:,:) / ZEXN(:,2:klev+1) -CALL VERTICAL_EXTEND(ZEXN,klev) -CALL VERTICAL_EXTEND(ZTHETA,klev) - -!TODO check in Meso-NH how zz_mass and zz_flux are initialized outside of the physical domain -zs(:) = pphis(:)/PHYEX%CST%XG -zz_mass(:,2:klev+1) = pphi(:,:) / PHYEX%CST%XG -zz_mass(:,1) = 2*zs-zz_mass(:,2) -zz_mass(:,klev+2)=2.*zz_mass(:,klev+1)-zz_mass(:,klev) - -do k=2, klev+2 - zz_flux(:,k)=(zz_mass(:,k-1)+zz_mass(:,k))/2. -enddo -zz_flux(:,1)=2*zz_mass(:,1)-zz_flux(:,2) - -!zdzf is the distance between two consecutive flux points (expressed on mass points) -do k=1,klev+1 - zdzf(:,k)=zz_flux(:,k+1)-zz_flux(:,k) -enddo -zdzf(:,klev+2)=(zz_mass(:,klev+2)-zz_flux(:,klev+2))*2. - -!zdzm distance between two consecutive mass points (expressed on flux points) -do k=2,klev+2 - zdzm(:,k)=zz_mass(:,k)-zz_mass(:,k-1) -enddo -zdzm(:,1)=(zz_mass(:,1)-zz_flux(:,1))*2. - -ZPABST(:,2:klev+1) = pplay(:,:) -ZRHOD(:,2:klev+1)=ZPABST(:,2:klev+1)/(t*(PHYEX%CST%XRD+ZRX(:,2:klev+1,1)*PHYEX%CST%XRV)) -DO k=2,klev+1 - PRHODJ(:,k) = ZRHOD(:,k) * (zdzf(:,k)*cell_area(:)) -END DO -PTHVREF(:,:) = ZTHETA(:,:) * (1. + PHYEX%CST%XRV/PHYEX%CST%XRD * ZRX(:,:,1)) * ZQDM(:,:) - -CALL VERTICAL_EXTEND(ZPABST,klev) -CALL VERTICAL_EXTEND(PRHODJ,klev) -CALL VERTICAL_EXTEND(ZRHOD,klev) -CALL VERTICAL_EXTEND(ZUT,klev) -CALL VERTICAL_EXTEND(ZVT,klev) - !------------------------------------------------------------ -! Tendencies +! Calculs !------------------------------------------------------------ -!For Meso-NH, initialia values for the tendencies are filled with -!a pseudo-tendecy computed by dividing the state variable by the time step -!This mechanism enables the possibility for the different parametrisations -!to guess the value at the end of the time step. -!For the wind components, we could do the same way but it is not needed -!as the parametrisations don't use the S varaible to guess the futur value of the wind. -ZRXS(:,:,:) = ZRX(:,:,:)/pdtphys -ZTHETAS(:,:)=ZTHETA(:,:)/pdtphys -ZTKES(:,:)=PTKEM(:,:)/pdtphys -!To compute the actual tendency, we save the initial values of these variables -ZRXS0(:,:,:) = ZRXS(:,:,:) -ZTHETAS0(:,:)=ZTHETAS(:,:) -ZTKES0(:,:)=ZTKES(:,:) -!------------------------------------------------------------ -! Adjustment -!------------------------------------------------------------ -! -ZSRC(:,:) = 0. -ZSIGQSAT=PHYEX%NEBN%VSIGQSAT -CALL ICE_ADJUST (D, PHYEX%CST, PHYEX%RAIN_ICE_PARAMN, PHYEX%NEBN, PHYEX%TURBN, PHYEX%PARAM_ICEN, & - &PHYEX%MISC%TBUCONF, KRR, & - &'ADJU', & - &pdtphys, ZSIGQSAT, & - &PRHODJ, ZEXN, ZRHOD, PSIGS, .FALSE., zmfconv, & - &ZPABST, ZZ_MASS, & - &ZEXN, PCF_MF, PRC_MF, PRI_MF, & - &ZICLDFR, ZWCLDFR, ZSSIO, ZSSIU, ZIFR, & - &ZRX(:,:,1), ZRX(:,:,2), ZRXS(:,:,1), ZRXS(:,:,2), ZTHETA, ZTHETAS, & - &PHYEX%MISC%COMPUTE_SRC, ZSRC, ZCLDFR, & - &ZRX(:,:,3), ZRX(:,:,4), ZRXS(:,:,4), ZRX(:,:,5), ZRX(:,:,6), & - &PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET, & - &ZICE_CLD_WGT, & - &PHLC_HRC=ZHLC_HRC, PHLC_HCF=ZHLC_HCF, PHLI_HRI=ZHLI_HRI, PHLI_HCF=ZHLI_HCF ) -! -!Variables are updated with their adjusted values (to be used by the other parametrisations) -ZTHETA(:,:)=ZTHETAS(:,:)*pdtphys -ZRX(:,:,:)=ZRXS(:,:,:)*pdtphys -t_adj=ZTHETA(:,2:klev+1)*ZEXN(:,2:klev+1) -qv_adj=ZRX(:,2:klev+1,1)*zqdm(:,2:klev+1) -ql_adj=ZRX(:,2:klev+1,2)*zqdm(:,2:klev+1) -qi_adj=ZRX(:,2:klev+1,4)*zqdm(:,2:klev+1) -qr_adj=ZRX(:,2:klev+1,3)*zqdm(:,2:klev+1) -qs_adj=ZRX(:,2:klev+1,5)*zqdm(:,2:klev+1) -qg_adj=ZRX(:,2:klev+1,6)*zqdm(:,2:klev+1) -ZRHOD(:,2:klev+1)=ZPABST(:,2:klev+1)/(t*(PHYEX%CST%XRD+ZRX(:,2:klev+1,1)*PHYEX%CST%XRV)) -CALL VERTICAL_EXTEND(ZRHOD,klev) -! -!------------------------------------------------------------ -! Radiation -!------------------------------------------------------------ -! -!AI a nettoyer apres -! calcul de qsat Fred mai 2024 -zqsat(:,:) = EXP( XALPW - XBETAW / t(:,:) - XGAMW * ALOG( t(:,:) ) ) -zqsat(:,:) = XRD / XRV * zqsat(:,:) / ( pplay(:,:) - zqsat(:,:)) -! -if (iflag_radia.ge.1) then -! -!ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors -!Cloud fraction is available in the ZCLDFR variable -! -! AI attention : dans un 1er temps read_climoz=0 -! Update ozone if day change - !print*,'solarlong0, itap, lmt_pas, days_elapsed : ', solarlong0, itap, lmt_pas, days_elapsed - !stop - IF (MOD(itap-1,lmt_pas) == 0) THEN -! IF (read_climoz <= 0) THEN - ! Once per day, update ozone from Royer: - IF (solarlong0<-999.) then - ! Generic case with evolvoing season - zzzo=real(days_elapsed+1) - ELSE IF (abs(solarlong0-1000.)<1.e-4) then - ! Particular case with annual mean insolation - zzzo=real(90) ! could be revisited - IF (read_climoz/=-1) THEN - abort_message ='read_climoz=-1 is recommended when ' & - // 'solarlong0=1000.' - CALL abort_physic (modname,abort_message,1) - ENDIF - ELSE - ! Case where the season is imposed with solarlong0 - zzzo=real(90) ! could be revisited - ENDIF - zzzo=real(days_elapsed+1) - wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzzo) - !print*,'wo = ', wo -! ELSE -! !--- ro3i = elapsed days number since current year 1st january, 0h -! ro3i=days_elapsed+jh_cur-jh_1jan -! !--- scaling for old style files (360 records) -! IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len -! IF(adjust_tropopause) THEN -! CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), & -! ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), & -! time_climoz , longitude_deg, latitude_deg, & -! dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)) -! ELSE -! CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), & -! ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), & -! time_climoz ) -! ENDIF -! ! Convert from mole fraction of ozone to column density of ozone in a -! ! cell, in kDU: -! FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd & -! * zmasse / dobson_u / 1e3 -! ! (By regridding ozone values for LMDZ only once a day, we -! ! have already neglected the variation of pressure in one -! ! day. So do not recompute "wo" at each time step even if -! ! "zmasse" changes a little.) -! ENDIF - ENDIF -!========================================================================= -! Calculs de l'orbite. -! Necessaires pour le rayonnement et la surface (calcul de l'albedo). -! doit donc etre plac\'e avant radlwsw et pbl_surface - !CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq) - CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq) - print*,'mth_eq, day_eq, jD_eq, jH_cur =', mth_eq, day_eq, jD_eq, jH_cur - day_since_equinox = (jD_cur + jH_cur) - jD_eq - !choix entre calcul de la longitude solaire vraie ou valeur fixee = solarlong0 - IF (solarlong0<-999.) THEN - IF (new_orbit) THEN - ! calcul selon la routine utilisee pour les planetes - CALL solarlong(day_since_equinox, zlongi, dist) - print*,'day_since_equinox, zlongi, dist : ', day_since_equinox, zlongi, dist - ELSE - ! calcul selon la routine utilisee pour l'AR4 - CALL orbite(REAL(days_elapsed+1),zlongi,dist) - ENDIF - ELSE - zlongi=solarlong0 ! longitude solaire vraie - dist=1. ! distance au soleil / moyenne - ENDIF - !IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist -! AI a refaire en definissant les itap, itaprad - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! ============================ - ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur - ! l'annee a partir d'une formule analytique. - ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et - ! non nul aux poles. - IF (abs(solarlong0-1000.)<1.e-4) THEN - CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, & - latitude_deg,longitude_deg,rmu0,fract) - swradcorr(:) = 1.0 - ! JrNt(:) = 1.0 - zrmu0(:) = rmu0(:) - ELSE - ! recode par Olivier Boucher en sept 2015 - SELECT CASE (iflag_cycle_diurne) - CASE(0) - ! Sans cycle diurne - CALL angle(zlongi, latitude_deg, fract, rmu0) - swradcorr = 1.0 - ! JrNt = 1.0 - zrmu0 = rmu0 - CASE(1) - ! Avec cycle diurne sans application des poids - ! bit comparable a l ancienne formulation cycle_diurne=true - ! on integre entre gmtime et gmtime+radpas - zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s) - CALL zenang(zlongi,jH_cur,0.0,zdtime, & - latitude_deg,longitude_deg,rmu0,fract) - zrmu0 = rmu0 - swradcorr = 1.0 - ! Calcul du flag jour-nuit - ! JrNt = 0.0 - !WHERE (fract.GT.0.0) JrNt = 1.0 - CASE(2) - ! Avec cycle diurne sans application des poids - ! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1) - ! Comme cette routine est appele a tous les pas de temps de - ! la physique meme si le rayonnement n'est pas appele je - ! remonte en arriere les radpas-1 pas de temps - ! suivant. Petite ruse avec MOD pour prendre en compte le - ! premier pas de temps de la physique pendant lequel - ! itaprad=0 - zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1) - zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1) - CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, & - latitude_deg,longitude_deg,rmu0,fract) - ! - ! Calcul des poids - ! - zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le - zdtime2=0.0 !--pas de temps de la physique qui se termine - CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, & - latitude_deg,longitude_deg,zrmu0,zfract) - swradcorr = 0.0 - WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) & - swradcorr=zfract/fract*zrmu0/rmu0 - ! Calcul du flag jour-nuit - !JrNt = 0.0 - !WHERE (zfract.GT.0.0) JrNt = 1.0 - END SELECT - ENDIF - !print*,'fract = ',fract - !sza_o = ACOS (rmu0) *180./pi - - ! Calcul de l'ensoleillement : - IF (MOD(itaprad,radpas).EQ.0) THEN - namelist_ecrad_file='namelist_ecrad' - ! - ! Q0 gestion des pas de temps : - ! itap, itaprad ????? - ! - ! Q1 champs lies a la surface ? sorties appel au modele de surface - ! zxtsol temp au sol : dvrait etre recuperee de la partie surface - ! albsol_dir, albsol_dif aussi - albsol_dir = 1. ! AI attention - albsol_dif = 1. ! AI attention - zxtsol(:) = t(:,1) ! AI attention - ! - ! Q2 ozone ? lmdz ok a completer - ! wo : pour l'instant appel ozonecm - ! - ! Q3 orbite ? lmdz ok a comleter - ! dist, rmu0 et fract - ! pour l'instant cas solarlong0<-999. - ! comment introduire les pas de temps rayonnement - ! - ! Q3 proprietes optiques des nuages ? call cloud_optics_prop.F90 ? non necessaires pour ecrad - ! epaisseur optique, emissivite : cldfrarad, cldemirad, cldtaurad , non necessaire pour ecrad - ! rayons effectifs : ref_liq, ref_ice ????? - ! cloud liq/Ice water mixing contenent (kg/kg) ???? - ! les rayons effectifs seront calcules en utilisant la routine de l'IFS - ! en fonct des contenus en eau, ..... - cldemirad = 1. - cldtaurad = 5. - cldtaupirad = 5. - ref_liq = 1. - ref_ice = 1. - ref_liq_pi = 1. - ref_ice_pi = 1. - ! - ! Q4 aerosols ? - ! simplifier l'ecriture : mettre tous les flags et parametres dans aero_mod.F90 ? - ! ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, & - ! flag_aerosol, flag_aerosol_strat, flag_aer_feedback, & - ! tau_aero, piz_aero, cg_aero, & - ! tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & - ! tau_aero_lw_rrtm, & - ! - ! Q5 thermo : zqsat - ! zqsat specific humidity at liq saturation => - ! calcul de Fred ou a calculer dans ecrad enf fct de t et p - ! flwc, fiwc : qx( 2) et qx( 3) - - print*,'Entree radlwsw a itap : ', itap - CALL radlwsw & - (debut, dist, rmu0, fract, & - & paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, & - & t,qx(:,:,1),wo, & - & ZCLDFR, cldemirad, cldtaurad, & - & ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, & - & flag_aerosol, flag_aerosol_strat, flag_aer_feedback, & - & tau_aero, piz_aero, cg_aero, & - & tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & - ! Rajoute par OB pour RRTM - & tau_aero_lw_rrtm, & - & cldtaupirad, m_allaer, & - ! - ! & zqsat, flwc, fiwc, & - & zqsat, qx(:,:,2), qx(:,:,2), & - & ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, & - & namelist_ecrad_file, & - & heat,heat0,cool,cool0,albpla, & - & heat_volc,cool_volc, & - & topsw,toplw,solsw,solswfdiff,sollw, & - & sollwdown, & - & topsw0,toplw0,solsw0,sollw0, & - & lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, & - & swdnc0, swdn0, swdn, swupc0, swup0, swup, & - & topswad_aero, solswad_aero, & - & topswai_aero, solswai_aero, & - & topswad0_aero, solswad0_aero, & - & topsw_aero, topsw0_aero, & - & solsw_aero, solsw0_aero, & - & topswcf_aero, solswcf_aero, & - !-C. Kleinschmitt for LW diagnostics - & toplwad_aero, sollwad_aero,& - & toplwai_aero, sollwai_aero, & - & toplwad0_aero, sollwad0_aero,& - !-end - & ZLWFT0_i, ZFLDN0, ZFLUP0, & - & ZSWFT0_i, ZFSDN0, ZFSUP0, & - & cloud_cover_sw) - itaprad = 0 - ENDIF ! IF (MOD(itaprad,radpas).EQ.0) - itaprad = itaprad + 1 - - IF (iflag_radia.eq.0) THEN - heat=0. - cool=0. - sollw=0. ! MPL 01032011 - solsw=0. - radsol=0. - swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars - swup0=0. - lwup=0. - lwup0=0. - lwdn=0. - lwdn0=0. - ENDIF - ! Ajouter la tendance des rayonnements (tous les pas) - ! avec une correction pour le cycle diurne dans le SW - ! - DO k=1, klev - ! AI attention - !d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY - !d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY - !d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY - !d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY - d_t_swr(:,k)=swradcorr(:)*heat(:,k)/RDAY - d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)/RDAY - d_t_lwr(:,k)=-cool(:,k)/RDAY - d_t_lw0(:,k)=-cool0(:,k)/RDAY - ENDDO - - !t_seri(:,:) = t_seri(:,:) + d_t_swr(:,:) + d_t_lwr(:,:) - d_t = d_t + d_t_swr + d_t_lwr - - !CALL writefield_phy('d_t_swr',d_t_swr,klev) - !CALL writefield_phy('d_t_lwr',d_t_lwr,klev) - !CALL writefield_phy('temp',t,klev) - -endif ! iflag_radia a nettoyer -! -! AI passer peut-etre la surface avant radiation -! albedo, temp surf, ... pour le rayonnement -! -!ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors -!Cloud fraction is available in the ZCLDFR variable - -! -!------------------------------------------------------------ -! Surface -!------------------------------------------------------------ -! -! AI ! compute tendencies to return to the dynamics: ! "friction" on the first layer -if (debut) then - tsrf(1:nlon)=t(1:nlon,1) - do isoil=1,nsoil - tsoil(1:nlon,isoil)=tsrf(1:nlon) - !tsoil(1:nlon,isoil)=302. - enddo - rpi=2.*asin(1.) - timesec=5*3600 - stefan=5.67e-8 - fsw0=900. - call getin_p('iecri',iecri) - call getin_p('iflag_surface',iflag_surface) -endif +d_u(1:klon,1)=-u(1:klon,1)/86400. +d_v(1:klon,1)=-v(1:klon,1)/86400. +! newtonian relaxation towards temp_newton() +do k=1,klev + temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 + d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5 +enddo -timesec=timesec+pdtphys -! Computation of qsats. Done even it is not used for ouptuts -! AI -!qsats(:) = EXP( XALPW - XBETAW / t(:,1) - XGAMW * ALOG( t(:,1) ) ) -!qsats(:) = XRD / XRV * qsats(:) / ( pplay(:,1) - qsats(:) ) -qsats(:) = zqsat(:,1) - -if (iflag_surface == 0 ) then - ! Latent and sensible heat fluxes imposed in 1D (exemple Arm cumulus) - ! radiative fluxes set to 0. - ! capcal and fluxgrd=0. replace generic_soil and allow to compute a surface temperature - ! A case close to Rico can be obtained by imposing fsens=-5E-3 et flat=-6E-5 - capcal(:)=1.e6 - fluxgrd(:)=0. - sensible(:)=-fsens - latent(:)=-flat - swdnsrf(:)=0. - lwdnsrf(:)=0. - fluxt(1:klon)=sensible/PHYEX%CST%XCPD - fluxq(1:klon)=latent(:)/PHYEX%CST%XLVTT - fluxu(1:klon)=-0.01*u(1:klon,1) - fluxv(1:klon)=-0.01*v(1:klon,1) -else - ! AI - if (iflag_radia.eq.0) then - ! Using the generic soil model - if ( nlon == 1 ) then - ! Imposed radiative fluxes to simulate a 1D case close to Arm Cumulus - swdnsrf(:)=fsw0*max(cos(2*rpi*(timesec-43200.)/(1.15*86400.)),0.) - lwdnsrf(:)=400. - else - ! Model forced in latitude by an idelaized radiative flux that depends on latitude - ! only. Used for 3D tests without radiation. - swdnsrf(:)=600*cos(latitude(:))+100 - lwdnsrf(:)=0. - endif - else - swdnsrf(:)=solsw(:) - lwdnsrf(:)=sollwdown(:) - endif - ! updating the soil temperature (below the surface) from the surface temperature - ! en returning the surface conductive flux and a heat capacity (accounting for the - ! sentivity of the conductive flux to the surface temperature. - call soil_conduct(klon,nsoil,debut,pdtphys,tsrf,tsoil,capcal,fluxgrd,zout) - - ! tsoil_out = tsoil with the vertical dimension of atmosphere outputs - - ! Computing the turbulent fluxes Fqx = rho w'x' = rho (kappa/ln(z/z0x))^2 [ x_surface - x ]. Upward - call soil_fluxes(klon,tsrf,qsats,pphi(:,1)/rg,ZRHOD(:,2),u(:,1),v(:,1),t(:,1),qx(:,1,1),fluxu,fluxv,fluxt,fluxq) - sensible(:)=fluxt(:)*PHYEX%CST%XCPD - latent(:)=fluxq(:)*PHYEX%CST%XLVTT -endif -tsoil_out=0. ; tsoil_out(:,1:min(nsoil,nlev))=tsoil(:,1:min(nsoil,nlev)) - -! Computing the surface temperature -lwupsrf=stefan*tsrf(:)**4 -tsrf(:)=tsrf(:)+pdtphys*(-latent-sensible+lwdnsrf+swdnsrf-lwupsrf+fluxgrd)/capcal(:) - -PSFTH(:) = fluxt(:) -PSFRV(:) = fluxq(:) -PSFU(:) = fluxu(:) -PSFV(:) = fluxv(:) -PSFSV(:,:) = 0. -! -!------------------------------------------------------------ -! Shallow convection -!------------------------------------------------------------ -! -CALL SHALLOW_MF(D, PHYEX%CST, PHYEX%NEBN, PHYEX%PARAM_MFSHALLN, PHYEX%TURBN, PHYEX%CSTURB, & - &KRR=KRR, KRRL=KRRL, KRRI=KRRI, KSV=KSV, & - &ONOMIXLG=PHYEX%MISC%ONOMIXLG,KSV_LGBEG=PHYEX%MISC%KSV_LGBEG,KSV_LGEND=PHYEX%MISC%KSV_LGEND, & - &PTSTEP=pdtphys, & - &PDZZ=zdzm(:,:),PZZ=zz_mass(:,:), & - &PRHODJ=PRHODJ(:,:),PRHODREF=ZRHOD(:,:), & - &PPABSM=ZPABST(:,:),PEXNM=ZEXN(:,:), & - &PSFTH=PSFTH(:),PSFRV=PSFRV(:), & - &PTHM=ZTHETA(:,:),PRM=ZRX(:,:,:),PUM=ZUT(:,:),PVM=ZVT(:,:), & - &PTKEM=PTKEM(:,:),PSVM=ZSVT(:,:,:), & - &PDUDT_MF=PDUDT_MF(:,:),PDVDT_MF=PDVDT_MF(:,:), & - &PDTHLDT_MF=PDTHLDT_MF(:,:),PDRTDT_MF=PDRTDT_MF(:,:),PDSVDT_MF=PDSVDT_MF(:,:,:), & - &PSIGMF=PSIGMF(:,:),PRC_MF=PRC_MF(:,:),PRI_MF=PRI_MF(:,:),PCF_MF=PCF_MF(:,:), & - &PFLXZTHVMF=ZFLXZTHVMF(:,:), & - &PFLXZTHMF=ZFLXZTHMF(:,:),PFLXZRMF=ZFLXZRMF(:,:),PFLXZUMF=ZFLXZUMF(:,:),PFLXZVMF=ZFLXZVMF(:,:), & - &PTHL_UP=PTHL_UP(:,:),PRT_UP=PRT_UP(:,:),PRV_UP=PRV_UP(:,:), & - &PRC_UP=PRC_UP(:,:),PRI_UP=PRI_UP(:,:), & - &PU_UP=PU_UP(:,:), PV_UP=PV_UP(:,:), PTHV_UP=PTHV_UP(:,:), PW_UP=PW_UP(:,:), & - &PFRAC_UP=PFRAC_UP(:,:),PEMF=PEMF(:,:),PDETR=PDETR(:,:),PENTR=PENTR(:,:), & - &KKLCL=IKLCL(:),KKETL=IKETL(:),KKCTL=IKCTL(:),PDX=1000.0,PDY=1000.0,KBUDGETS=PHYEX%MISC%NBUDGET ) - -! Add tendencies of shallow to total physics tendency -d_u(:,1:klev) = d_u(:,1:klev) + PDUDT_MF(:,2:klev+1) -d_v(:,1:klev) = d_v(:,1:klev) + PDVDT_MF(:,2:klev+1) -ZRXS(:,:,1)=ZRXS(:,:,1)+PDRTDT_MF(:,:) -ZTHETAS(:,:)=ZTHETAS(:,:)+PDTHLDT_MF(:,:) -! TODO add SV tendencies -! -!------------------------------------------------------------ -! Turbulence -!------------------------------------------------------------ -! out tendencies -ZRUS(:,:) = 0. -ZRVS(:,:) = 0. -ZRWS(:,:) = 0. -ZRSVS(:,:,:) = 0. -ZRTKES(:,:) = ZTKES(:,:) * PRHODJ(:,:) -DO JRR=1, KRR - ZRRS(:,:,JRR) = ZRXS(:,:,JRR) * PRHODJ(:,:) -ENDDO -ZRTHS(:,:) = ZTHETAS(:,:) * PRHODJ(:,:) -CALL TURB(PHYEX%CST, PHYEX%CSTURB, PHYEX%MISC%TBUCONF, PHYEX%TURBN, PHYEX%NEBN, D, PHYEX%MISC%TLES, & - & KRR, KRRL, KRRI, PHYEX%MISC%HLBCX, PHYEX%MISC%HLBCY, PHYEX%MISC%KGRADIENTS, PHYEX%MISC%KHALO, & - & PHYEX%TURBN%NTURBSPLIT, PHYEX%TURBN%LCLOUDMODIFLM, KSV, PHYEX%MISC%KSV_LGBEG, PHYEX%MISC%KSV_LGEND, & - & PHYEX%MISC%KSV_LIMA_NR, PHYEX%MISC%KSV_LIMA_NS, PHYEX%MISC%KSV_LIMA_NG, PHYEX%MISC%KSV_LIMA_NH, & - & PHYEX%MISC%O2D, PHYEX%MISC%ONOMIXLG, PHYEX%MISC%OFLAT, PHYEX%MISC%OCOUPLES, & - & PHYEX%MISC%OBLOWSNOW,PHYEX%MISC%OIBM, & - & PHYEX%MISC%OFLYER, PHYEX%MISC%COMPUTE_SRC, PHYEX%MISC%PRSNOW, & - & PHYEX%MISC%OOCEAN, PHYEX%MISC%ODEEPOC, PHYEX%MISC%ODIAG_IN_RUN, & - & PHYEX%TURBN%CTURBLEN_CLOUD, PHYEX%MISC%CMICRO, PHYEX%MISC%CELEC, & - & pdtphys,PHYEX%MISC%ZTFILE, & - & ZDXX(:,:),ZDYY(:,:),zdzm(:,:), & - & ZDZX(:,:),ZDZY(:,:),zz_flux(:,:), & - & ZDIRCOSXW(:),ZDIRCOSYW(:),ZDIRCOSZW(:),ZCOSSLOPE(:),ZSINSLOPE(:), & - & PRHODJ(:,:),PTHVREF(:,:), PHGRAD(:,:,:), zs(:), & - & PSFTH(:),PSFRV(:),PSFSV(:,:),PSFU(:),PSFV(:), & - & ZPABST(:,:),ZUT(:,:),ZVT(:,:),PWT(:,:),PTKEM(:,:),ZSVT(:,:,:),ZSRC(:,:), & - & PLENGTHM(:,:),PLENGTHH(:,:),MFMOIST(:,:), & - & ZBL_DEPTH(:),ZSBL_DEPTH(:), & - & ZCEI(:,:), PHYEX%TURBN%XCEI_MIN, PHYEX%TURBN%XCEI_MAX, PHYEX%TURBN%XCOEF_AMPL_SAT, & - & ZTHETA(:,:),ZRX(:,:,:), & - & ZRUS(:,:),ZRVS(:,:),ZRWS(:,:),ZRTHS(:,:),ZRRS(:,:,:),ZRSVS(:,:,:),ZRTKES(:,:), & - & PSIGS(:,:), & - & ZFLXZTHVMF(:,:),ZWTH(:,:),ZWRC(:,:),ZWSV(:,:,:),ZDP(:,:),ZTP(:,:),ZTDIFF(:,:),ZTDISS(:,:), & - & PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET ) -DO JRR=1, KRR - ZRXS(:,:,JRR) = ZRRS(:,:,JRR) / PRHODJ(:,:) -ENDDO -ZTHETAS(:,:) = ZRTHS(:,:) / PRHODJ(:,:) -ZTKES(:,:) = ZRTKES(:,:) / PRHODJ(:,:) -! Add tendencies of turb to total physics tendency -d_u(:,1:klev) = d_u(:,1:klev) + ZRUS(:,2:klev+1)/PRHODJ(:,2:klev+1) -d_v(:,1:klev) = d_v(:,1:klev) + ZRVS(:,2:klev+1)/PRHODJ(:,2:klev+1) -IF(PHYEX%PARAM_MFSHALLN%CMF_CLOUD=='STAT') THEN - PSIGS(:,:)=SQRT(PSIGS(:,:)**2 + PSIGMF(:,:)**2) -ENDIF -!------------------------------------------------------------ -! Microphysics -!------------------------------------------------------------ -ZSEA=1. -ZTOWN=0. -ZCIT=0. -ZTHVREFZIKB=0 -CALL RAIN_ICE (D, PHYEX%CST, PHYEX%PARAM_ICEN, PHYEX%RAIN_ICE_PARAMN, PHYEX%RAIN_ICE_DESCRN, & - PHYEX%ELEC_PARAM, PHYEX%ELEC_DESCR, PHYEX%MISC%TBUCONF, & - PHYEX%MISC%OELEC, PHYEX%MISC%OSEDIM_BEARD, & - ZTHVREFZIKB, & - pdtphys, KRR, ZEXN, & - zdzf, PRHODJ, ZRHOD, ZEXN, ZPABST, ZCIT, ZCLDFR, & - ZHLC_HRC, ZHLC_HCF, ZHLI_HRI, ZHLI_HCF, & - ztheta, ZRX(:,:,1), ZRX(:,:,2), ZRX(:,:,3), ZRX(:,:,4), ZRX(:,:,5), & - ZRX(:,:,6), zthetas, ZRXS(:,:,1), ZRXS(:,:,2), ZRXS(:,:,3), ZRXS(:,:,4), ZRXS(:,:,5), ZRXS(:,:,6), & - ZINPRC, ZINPRR, ZEVAP3D, & - ZINPRS, ZINPRG, ZINDEP, ZRAINFR, PSIGS, & - PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET, & - ZSEA, ZTOWN ) - -!------------------------------------------------------------ -! Tendencies and time evolution (values for next time step) -!------------------------------------------------------------ -! Tendencies, mixing ratio -> specific -d_qx(:,1:klev,1)=d_qx(:,1:klev,1) + (ZRXS(:,2:klev+1,1)-ZRXS0(:,2:klev+1,1))*ZQDM(:,2:klev+1) -d_qx(:,1:klev,2)=d_qx(:,1:klev,2) + (ZRXS(:,2:klev+1,2)-ZRXS0(:,2:klev+1,2))*ZQDM(:,2:klev+1) -d_qx(:,1:klev,3)=d_qx(:,1:klev,3) + (ZRXS(:,2:klev+1,4)-ZRXS0(:,2:klev+1,4))*ZQDM(:,2:klev+1) -d_qr(:,1:klev)=d_qr(:,1:klev) + (ZRXS(:,2:klev+1,3)-ZRXS0(:,2:klev+1,3))*ZQDM(:,2:klev+1) -d_qs(:,1:klev)=d_qs(:,1:klev) + (ZRXS(:,2:klev+1,5)-ZRXS0(:,2:klev+1,5))*ZQDM(:,2:klev+1) -d_qg(:,1:klev)=d_qg(:,1:klev) + (ZRXS(:,2:klev+1,6)-ZRXS0(:,2:klev+1,6))*ZQDM(:,2:klev+1) -! Tendency, theta -> T -d_t(:,1:klev)=d_t(:,1:klev) + (zthetas(:,2:klev+1)-zthetas0(:,2:klev+1))*ZEXN(:,2:klev+1) -! TKE -d_tke(:,1:klev)=d_tke(:,1:klev) + (ZTKES(:,2:klev+1) - ZTKES0(:,2:klev+1)) - -!Time evolution -ZQR(:,:)=ZQR(:,:)+d_qr(:,:)*pdtphys -ZQS(:,:)=ZQS(:,:)+d_qs(:,:)*pdtphys -ZQG(:,:)=ZQG(:,:)+d_qg(:,:)*pdtphys -PTKEM(:,2:klev+1)=PTKEM(:,2:klev+1)+d_tke(:,:)*pdtphys -! !------------------------------------------------------------ ! Entrees sorties !------------------------------------------------------------ -itap=itap+1 -!zjulian=zjulian+pdtphys/86400. -!print*,'avant mod(itap,iecri) , zjulian itap',zjulian,itap,iecri -if ( mod(itap,iecri) == 1 .or. iecri == 1 ) then -print*,'ECRITURE ITAP=',itap - call output_physiqex(debut,zjulian,pdtphys_,presnivs,paprs,u,v,t,qx,ZCLDFR,ZQR,ZQS,ZQG,PTKEM,ZTHETA) - call iophys_ecrit('tsrf', 1,'tsrf', ' ',tsrf) - call iophys_ecrit('capcal', 1,'capcal', ' ',capcal) - call iophys_ecrit('swdnsrf', 1,'swdnsrf', ' ',swdnsrf) - call iophys_ecrit('lwdnsrf', 1,'lwdnsrf', ' ',lwdnsrf) - call iophys_ecrit('lwupsrf', 1,'lwupsrf', ' ',lwupsrf) - call iophys_ecrit('sensible', 1,'sensible', ' ',sensible) - call iophys_ecrit('latent', 1,'latent', ' ',latent) - call iophys_ecrit('qsats', 1,'qsats', ' ',qsats) - call iophys_ecrit('tsoil', nlev,'tsoi', ' ',tsoil_out) -!AI -if (iflag_radia.ge.1) then - call iophys_ecrit('lwdn0', 1,'lwdn0', ' ',lwdn0) - call iophys_ecrit('lwdn', 1,'lwdn', ' ',lwdn) - call iophys_ecrit('lwup0', 1,'lwup0', ' ',lwup0) - call iophys_ecrit('lwup', 1,'lwup', ' ',lwup) - call iophys_ecrit('swdn0', 1,'swdn0', ' ',swdn0) - call iophys_ecrit('swdn', 1,'swdn', ' ',swdn) - call iophys_ecrit('swup0', 1,'swup0', ' ',swup0) - call iophys_ecrit('swup', 1,'swup', ' ',swup) - endif -endif + +call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t) + ! if lastcall, then it is time to write "restartphy.nc" file @@ -1227,14 +138,4 @@ end subroutine physiqex -! -SUBROUTINE VERTICAL_EXTEND(PX,KLEV) - - ! fill extra vetical levels to fit MNH interface - -REAL, DIMENSION(:,:), INTENT(INOUT) :: PX -INTEGER, INTENT(IN) :: KLEV -PX(:,1 )= PX(:,2) -PX(:,KLEV+2)= PX(:,KLEV+1) -END SUBROUTINE VERTICAL_EXTEND END MODULE physiqex_mod