Changeset 5197 for LMDZ6/trunk/libf
- Timestamp:
- Sep 18, 2024, 3:55:53 PM (2 months ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/physiqex_mod.F90
r4658 r5197 1 ! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $2 1 MODULE physiqex_mod 3 4 2 IMPLICIT NONE 5 6 3 CONTAINS 7 4 5 ! ------------------------------------------------------------------------------------------------ 6 ! 7 ! Main authors : F. Hourdin, S. Riette, Q. Libois, A. Idelkadi 8 ! -------------- 9 ! 10 ! Object : A version of the LMDZ physics that include PHYEX 11 ! ------- Tu be used in particular to be coupled to non hydrostatic version of the dynamical 12 ! core dynamico or its limited area version. 13 ! 14 ! 15 ! Development steps : 16 ! ------------------- 17 ! 18 ! physiqex_mod.F90 has been introducued in phylmd as an empty driver for the physics (with only 19 ! a newtonian relaxation for temperature and drag at the surface). 20 ! This driver has been used to introduce the PHYEX package, an externailzed physics package coming 21 ! from MesoNH, representing turbulence, convection, clouds and precipitations. 22 ! This work has been done for a large part in Dephy sessions (Oleron, 2022, Fréjus 2023, Lège 2024). 23 ! 24 ! Fréjus 2023 : first version runing in 1D for Arm cumulus 25 ! Lège 2024 : first versions with a generic soil model (developped at this occasion) and ECrad 26 ! 27 ! 28 ! To do list 29 ! --------- 30 ! To do list pour un branchement plus propre 31 ! * PHYEX considère toutes les espèces microphysiques et la tke comme pronostiques, des termes de tendances 32 ! sont calculés. L'avance temporelle est faite en fin de pas de temps mais pourrait être déplacée au dessus si 33 ! les tendances de ces variables passaient par l'interface 34 ! * PHYEX a besoin de variables avec mémoire d'un pas de temps à l'autre (en plus des variables pronostiques 35 ! décrites juste au dessus). Ce sont les tableaux ALLOCATABLE. Ces variables pourraient devenir 36 ! des traceurs. 37 ! * La variable ZDZMIN est ici calculée avec un MINVAL qui conduira à des résultats différents 38 ! lorsque la distribution (noeuds/procs) sera différente. Cette variable est utile pour déterminer 39 ! un critère CFL pour la sédimentation des précipitations. Il suffit donc d'avoir une valeur 40 ! approchante. 41 ! * Certains allocatable sont en klev, d'autres en klev+2. Il est possible de changer ceci mais il faut gérer 42 ! les recopies pour que les params voient des tableaux en klev+2 (en effectuant des recopies des 43 ! niveaux extrêmes comme fait pour le vent par exemple) 44 ! * L'eau tombée en surface (précipitations, sédimentation du nuage, terme de dépôt) se trouve dans 45 ! les variables ZINPRC, ZINPRR, ZINPRS, ZINPRG 46 ! 47 !------------------------------------------------------------------------------------------------- 48 49 8 50 SUBROUTINE physiqex (nlon,nlev, & 9 & debut,lafin,pdtphys , &51 & debut,lafin,pdtphys_, & 10 52 & paprs,pplay,pphi,pphis,presnivs, & 11 53 & u,v,rot,t,qx, & … … 15 57 USE dimphy, only : klon,klev 16 58 USE infotrac_phy, only : nqtot 17 USE geometry_mod, only : latitude 59 USE geometry_mod, only : latitude, cell_area, latitude_deg, longitude_deg 18 60 ! USE comcstphy, only : rg 19 61 USE ioipsl, only : ymds2ju 20 USE phys_state_var_mod, only : phys_state_var_init 62 !USE phys_state_var_mod, only : phys_state_var_init 63 USE phys_state_var_mod 21 64 USE phyetat0_mod, only: phyetat0 22 65 USE output_physiqex_mod, ONLY: output_physiqex 66 !!!!!! AI 04-2024 67 USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time 68 USE limit_read_mod, ONLY : init_limit_read 69 USE conf_phys_m, only: conf_phys 70 USE ozonecm_m, only: ozonecm 71 USE aero_mod 72 USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, & 73 year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour, calend 74 USE radlwsw_m, only: radlwsw 75 USE MODD_CST 76 USE phys_output_var_mod, ONLY : phys_output_var_init, cloud_cover_sw 77 USE write_field_phy 78 USE phyaqua_mod, only: zenang_an 79 !!!!!! 80 ! PHYEX internal modules 81 USE MODE_INIT_PHYEX, ONLY: INIT_PHYEX, FILL_DIMPHYEX 82 USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t 83 USE MODD_PHYEX, ONLY: PHYEX_t 84 USE MODI_INI_PHYEX, ONLY: INI_PHYEX 85 USE MODI_ICE_ADJUST, ONLY: ICE_ADJUST 86 USE MODI_RAIN_ICE, ONLY: RAIN_ICE 87 USE MODI_TURB 88 USE MODI_SHALLOW_MF 89 90 USE MODD_CST 91 USE ioipsl_getin_p_mod, ONLY : getin_p 92 USE generic_soil_ini, ONLY : soil_ini 93 USE generic_soil_conduct, ONLY : soil_conduct 94 USE generic_soil_fluxes, ONLY : soil_fluxes 23 95 24 96 IMPLICIT none … … 26 98 ! Routine argument: 27 99 ! 28 29 100 integer,intent(in) :: nlon ! number of atmospheric colums 30 101 integer,intent(in) :: nlev ! number of vertical levels (should be =klev) 31 102 logical,intent(in) :: debut ! signals first call to physics 32 103 logical,intent(in) :: lafin ! signals last call to physics 33 real,intent(in) :: pdtphys ! physics time step (s)104 real,intent(in) :: pdtphys_ ! physics time step (s) 34 105 real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) 35 106 real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) … … 48 119 real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers 49 120 real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure 50 51 ! include "clesphys.h" 121 real :: d_qr(klon, klev), d_qs(klon, klev), d_qg(klon, klev) ! tendency for rain, snow, graupel 122 real :: d_tke(klon, klev) 123 124 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 125 ! AI 04-2024 a nettoyer 126 ! 127 CHARACTER (LEN=20) :: modname='physiq_mod' 128 CHARACTER*80 abort_message 129 ! Declarations pas de temps 130 INTEGER, SAVE :: itap=0 ! compteur pour la physique 131 !$OMP THREADPRIVATE(itap) 132 INTEGER lmt_pas 133 SAVE lmt_pas ! frequence de mise a jour 134 !$OMP THREADPRIVATE(lmt_pas) 135 REAL zdtime, zdtime1, zdtime2, zlongi 136 ! Declarations pour appel de conf_phys 137 LOGICAL ok_hf 138 SAVE ok_hf 139 !$OMP THREADPRIVATE(ok_hf) 140 LOGICAL ok_journe ! sortir le fichier journalier 141 SAVE ok_journe 142 !$OMP THREADPRIVATE(ok_journe) 143 LOGICAL ok_mensuel ! sortir le fichier mensuel 144 SAVE ok_mensuel 145 !$OMP THREADPRIVATE(ok_mensuel) 146 LOGICAL ok_instan ! sortir le fichier instantane 147 SAVE ok_instan 148 !$OMP THREADPRIVATE(ok_instan) 149 LOGICAL ok_LES ! sortir le fichier LES 150 SAVE ok_LES 151 !$OMP THREADPRIVATE(ok_LES) 152 LOGICAL callstats ! sortir le fichier stats 153 SAVE callstats 154 !$OMP THREADPRIVATE(callstats) 155 LOGICAL ok_region ! sortir le fichier regional 156 PARAMETER (ok_region=.FALSE.) 157 REAL seuil_inversion 158 SAVE seuil_inversion 159 !$OMP THREADPRIVATE(seuil_inversion) 160 ! Parametres lies au nouveau schema de nuages (SB, PDF) 161 REAL, SAVE :: fact_cldcon 162 REAL, SAVE :: facttemps 163 !$OMP THREADPRIVATE(fact_cldcon,facttemps) 164 LOGICAL, SAVE :: ok_newmicro 165 !$OMP THREADPRIVATE(ok_newmicro) 166 INTEGER, SAVE :: iflag_cld_th 167 !$OMP THREADPRIVATE(iflag_cld_th) 168 REAL ratqsbas,ratqshaut,tau_ratqs 169 SAVE ratqsbas,ratqshaut,tau_ratqs 170 !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs) 171 ! ozone 172 INTEGER,SAVE :: read_climoz ! Read ozone climatology 173 ! (let it keep the default OpenMP shared attribute) 174 ! Allowed values are 0, 1 and 2 175 ! 0: do not read an ozone climatology 176 ! 1: read a single ozone climatology that will be used day and night 177 ! 2: read two ozone climatologies, the average day and night 178 ! climatology and the daylight climatology 179 INTEGER,SAVE :: ncid_climoz 180 REAL zzzo 181 REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for 182 ! the ozone fields, old method. 183 ! 184 ! aerosols 185 ! params et flags aerosols 186 LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not 187 LOGICAL ok_alw ! Apply aerosol LW effect or not 188 LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013) 189 REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995) 190 SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1 191 !$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1) 192 LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA 193 ! false : lecture des aerosol dans un fichier 194 INTEGER, SAVE :: flag_aerosol 195 !$OMP THREADPRIVATE(flag_aerosol) 196 LOGICAL, SAVE :: flag_bc_internal_mixture 197 !$OMP THREADPRIVATE(flag_bc_internal_mixture) 198 ! 199 !--STRAT AEROSOL 200 INTEGER, SAVE :: flag_aerosol_strat 201 !$OMP THREADPRIVATE(flag_aerosol_strat) 202 ! 203 !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION 204 LOGICAL, SAVE :: flag_aer_feedback 205 !$OMP THREADPRIVATE(flag_aer_feedback) 206 ! 207 LOGICAL, SAVE :: chemistry_couple ! true : use INCA chemistry O3 208 ! false : use offline chemistry O3 209 !$OMP THREADPRIVATE(chemistry_couple) 210 LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques 211 !$OMP THREADPRIVATE(ok_volcan) 212 INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato 213 !$OMP THREADPRIVATE(flag_volc_surfstrat) 214 ! 215 INTEGER, SAVE :: iflag_radia=0 ! active ou non le rayonnement (MPL) 216 !$OMP THREADPRIVATE(iflag_radia) 217 ! 218 REAL, SAVE :: alp_offset 219 !$OMP THREADPRIVATE(alp_offset) 220 ! Declarations pour le rayonnement 221 ! Le rayonnement n'est pas calcule tous les pas, il faut donc 222 ! sauvegarder les sorties du rayonnement 223 ! Surface 224 REAL zxtsol(KLON) 225 INTEGER itaprad 226 SAVE itaprad 227 !$OMP THREADPRIVATE(itaprad) 228 ! Date equinoxe de printemps 229 !LOGICAL, parameter :: ok_rad=.true. ! a nettoyer 230 REAL,SAVE :: solarlong0 231 !$OMP THREADPRIVATE(solarlong0) 232 ! 6 bandes SW : argument non utilise dans radlwsw a nettoyer ? 233 REAL,DIMENSION(6), parameter :: SFRWL = (/ 1.28432794E-03, 0.12304168, 0.33106142, & 234 & 0.32870591, 0.18568763, 3.02191470E-02 /) 235 LOGICAL, parameter :: new_orbit = .TRUE. 236 INTEGER, parameter :: mth_eq=3, day_eq=21 237 REAL :: day_since_equinox, jD_eq 238 !REAL zdtime, zdtime1, zdtime2, zlongi 239 REAL dist, rmu0(klon), fract(klon), zrmu0(klon), zfract(klon) 240 CHARACTER(len=512) :: namelist_ecrad_file 241 ! 242 REAL, dimension(klon, klev) :: zqsat 243 REAL, dimension(klon, klev) :: ref_liq, ref_ice, ref_liq_pi, ref_ice_pi ! rayons effectifs des gout 244 REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique 245 REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique 246 REAL, dimension(klon, klev) :: cldemirad ! emissivite pour 247 REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse 248 ! aerosols AI attention a les declarer ailleurs dans phylmd : phys_output ou phy_stats et non phys_local 249 REAL, dimension(klon,klev+1) :: ZLWFT0_i, ZSWFT0_i,ZFLDN0, ZFLUP0, ZFSDN0, ZFSUP0 250 REAL, dimension(klon) :: oplwad_aero, sollwad_aero,& 251 toplwai_aero, sollwai_aero, & 252 toplwad0_aero, sollwad0_aero 253 REAL, dimension(klon,9) :: solsw_aero, solsw0_aero,topsw_aero, topsw0_aero 254 255 REAL, dimension(klon) :: & 256 topswcf_aero, solswcf_aero, & 257 solswad_aero, & 258 solswad0_aero, & 259 topswai_aero, solswai_aero, & 260 toplwad_aero, topswad0_aero, & 261 topswad_aero 262 REAL, DIMENSION(klon,klev,naero_tot) :: m_allaer 263 264 REAL,dimension(klon,klev) :: d_t_swr, d_t_sw0, d_t_lwr, d_t_lw0 265 ! 266 EXTERNAL suphel ! AI attention 267 ! 268 include "YOMCST.h" 269 include "clesphys.h" 270 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Fin AI 271 272 include "flux_arp.h" 273 52 274 INTEGER length 53 275 PARAMETER ( length = 100 ) … … 57 279 !$OMP THREADPRIVATE(clesphy0) 58 280 59 281 ! Saved variables 282 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PSIGS !variance of s 283 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PCF_MF, PRC_MF, PRI_MF !shallow convection cloud 284 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ZQR, ZQS, ZQG !rain, snow, graupel specifiq contents 285 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: PTKEM ! TKE 286 TYPE(DIMPHYEX_t),SAVE :: D 287 TYPE(PHYEX_t), SAVE :: PHYEX 288 ! 289 INTEGER, PARAMETER :: KRR=6, KRRL=2, KRRI=3, KSV=0 290 INTEGER :: JRR 291 ! Time-State variables and Sources variables 292 REAL, DIMENSION(klon,klev+2) :: ZUT, ZVT ! wind component on klev+2 293 REAL, DIMENSION(klon,klev+2) :: ZPABST ! absolute pressure 294 REAL, DIMENSION(klon,klev+2,KRR) :: ZRX,ZRXS,ZRXS0 ! qx and source of q from LMDZ to rt 295 REAL, DIMENSION(klon,klev+2) :: ZRHOD ! rho_dry 296 REAL, DIMENSION(klon,klev+2,0) :: ZSVT ! passive scal. var. 297 REAL, DIMENSION(klon,klev+2) :: PWT ! vertical wind velocity (used only for diagnostic) 298 REAL, DIMENSION(klon,klev+2) :: ZRUS,ZRVS,ZRWS,ZRTHS,ZRTKES ! sources of momentum, conservative potential temperature, Turb. Kin. Energy 299 REAL, DIMENSION(KLON,KLEV+2) :: ZTHETAS !tendency 300 REAL, DIMENSION(KLON,KLEV+2) :: ZTHETAS0 301 REAL, DIMENSION(KLON,KLEV+2) :: ZTKES 302 REAL, DIMENSION(KLON,KLEV+2) :: ZTKES0 303 ! Source terms for all water kinds, PRRS(:,:,:,1) is used for the conservative ! mixing ratio 304 REAL, DIMENSION(klon,klev+2,KRR) :: ZRRS 305 REAL, DIMENSION(klon,klev+2) :: PTHVREF ! Virtual Potential Temperature of the reference state 306 ! Adjustment variables 307 REAL, DIMENSION(KLON,KLEV+2) :: zqdm !1-qt=1/(1+rt) 308 REAL, DIMENSION(KLON,KLEV+2) :: zqt !mixing ratio and total specifiq content 309 REAL, DIMENSION(KLON,KLEV+2) :: ZTHETA, ZEXN !theta and exner function 310 REAL, DIMENSION(KLON,KLEV+2) :: zz_mass !altitude above ground of mass points 311 REAL, DIMENSION(KLON,KLEV+2) :: zz_flux !altitude above ground of flux points 312 REAL, DIMENSION(KLON,KLEV+2) :: zdzm !distance between two consecutive mass points (expressed on flux points) 313 REAL, DIMENSION(KLON,KLEV+2) :: zdzf !distance between two consecutive flux points (expressed on mass points) 314 REAL, DIMENSION(KLON) :: zs !surface orography 315 REAL, DIMENSION(KLON) :: zsigqsat !coeff for the extra term for s variance 316 REAL, DIMENSION(0) :: ZMFCONV 317 REAL, DIMENSION(KLON,KLEV+2) :: ZICLDFR, ZWCLDFR, ZSSIO, ZSSIU, ZIFR, ZICE_CLD_WGT !used only in HIRLAM config 318 REAL, DIMENSION(KLON,KLEV+2) :: ZSRC ! bar(s'rc') 319 REAL, DIMENSION(KLON,KLEV+2) :: ZCLDFR !cloud fraction 320 REAL, DIMENSION(KLON,KLEV+2) :: ZHLC_HRC, ZHLC_HCF, ZHLI_HRI, ZHLI_HCF !subgrid autoconversion 321 ! Rain_ice variables 322 real :: ZDZMIN 323 real, dimension(klon, klev+2) :: ZCIT !ice concentration 324 real, dimension(klon) :: ZINPRC, ZINPRR, ZINPRS, ZINPRG !precipitation flux at ground 325 real, dimension(klon, klev+2) :: ZEVAP3D !evaporation (diag) 326 real, dimension(klon, klev+2) :: ZRAINFR 327 real, dimension(klon) :: ZINDEP !deposition flux, already contained in ZINPRC 328 real, dimension(klon) :: ZSEA, ZTOWN !sea and town fractions in the frid cell 329 ! Turbulence variables 330 REAL, DIMENSION(klon,klev+2) :: PDXX,PDYY,PDZX,PDZY ! metric coefficients 331 REAL, DIMENSION(klon,klev+2) :: PZZ ! physical distance between 2 succesive grid points along the K direction 332 REAL, DIMENSION(klon) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW ! Director Cosinus along x, y and z directions at surface w-point 333 REAL, DIMENSION(klon) :: PCOSSLOPE ! cosinus of the angle between i and the slope vector 334 REAL, DIMENSION(klon) :: PSINSLOPE ! sinus of the angle between i and the slope vector 335 REAL, DIMENSION(klon,klev+2) :: PRHODJ ! dry density * Grid size 336 REAL, DIMENSION(0,0) :: MFMOIST ! moist mass flux dual scheme 337 REAL, DIMENSION(0,0,0) :: PHGRAD ! horizontal gradients 338 REAL, DIMENSION(klon) :: PSFTH,PSFRV,PSFU,PSFV ! normal surface fluxes of theta, Rv, (u,v) parallel to the orography 339 REAL, DIMENSION(klon,0) :: PSFSV ! normal surface fluxes of Scalar var. KSV=0 340 REAL, DIMENSION(klon) :: ZBL_DEPTH ! BL height for TOMS 341 REAL, DIMENSION(klon) :: ZSBL_DEPTH ! SBL depth for RMC01 342 REAL, DIMENSION(klon,klev+2) :: ZCEI ! Cloud Entrainment instability index to emphasize localy turbulent fluxes 343 REAL, DIMENSION(klon,klev+2,0):: ZRSVS ! Source terms for all passive scalar variables 344 REAL, DIMENSION(klon,klev+2) :: ZFLXZTHVMF ! MF contribution for vert. turb. transport used in the buoy. prod. of TKE 345 REAL, DIMENSION(klon,klev+2) :: ZWTH ! heat flux 346 REAL, DIMENSION(klon,klev+2) :: ZWRC ! cloud water flux 347 REAL, DIMENSION(klon,klev+2,0):: ZWSV ! scalar flux 348 REAL, DIMENSION(klon,klev+2) :: ZTP ! Thermal TKE production (MassFlux + turb) 349 REAL, DIMENSION(klon,klev+2) :: ZDP ! Dynamic TKE production 350 REAL, DIMENSION(klon,klev+2) :: ZTDIFF ! Diffusion TKE term 351 REAL, DIMENSION(klon,klev+2) :: ZTDISS ! Dissipation TKE term 352 REAL, DIMENSION(0,0) :: PLENGTHM, PLENGTHH ! length scale from vdfexcu (HARMONIE-AROME) 353 REAL :: ZTHVREFZIKB ! for electricity scheme interface 354 REAL, DIMENSION(klon,klev+2) :: ZDXX,ZDYY,ZDZX,ZDZY,ZZZ 355 REAL, DIMENSION(klon) :: ZDIRCOSXW,ZDIRCOSYW,ZDIRCOSZW,ZCOSSLOPE,ZSINSLOPE 356 ! Shallow variables 357 REAL, DIMENSION(klon,klev+2) :: PDUDT_MF ! tendency of U by massflux scheme 358 REAL, DIMENSION(klon,klev+2) :: PDVDT_MF ! tendency of V by massflux scheme 359 REAL, DIMENSION(klon,klev+2) :: PDTHLDT_MF ! tendency of thl by massflux scheme 360 REAL, DIMENSION(klon,klev+2) :: PDRTDT_MF ! tendency of rt by massflux scheme 361 REAL, DIMENSION(klon,klev+2,KSV):: PDSVDT_MF ! tendency of Sv by massflux scheme 362 REAL, DIMENSION(klon,klev+2) :: PSIGMF 363 REAL, DIMENSION(klon,klev+2) :: ZFLXZTHMF 364 REAL, DIMENSION(klon,klev+2) :: ZFLXZRMF 365 REAL, DIMENSION(klon,klev+2) :: ZFLXZUMF 366 REAL, DIMENSION(klon,klev+2) :: ZFLXZVMF 367 REAL, DIMENSION(klon,klev+2) :: PTHL_UP ! Thl updraft characteristics 368 REAL, DIMENSION(klon,klev+2) :: PRT_UP ! Rt updraft characteristics 369 REAL, DIMENSION(klon,klev+2) :: PRV_UP ! Vapor updraft characteristics 370 REAL, DIMENSION(klon,klev+2) :: PU_UP ! U wind updraft characteristics 371 REAL, DIMENSION(klon,klev+2) :: PV_UP ! V wind updraft characteristics 372 REAL, DIMENSION(klon,klev+2) :: PRC_UP ! cloud content updraft characteristics 373 REAL, DIMENSION(klon,klev+2) :: PRI_UP ! ice content updraft characteristics 374 REAL, DIMENSION(klon,klev+2) :: PTHV_UP ! Thv updraft characteristics 375 REAL, DIMENSION(klon,klev+2) :: PW_UP ! vertical speed updraft characteristics 376 REAL, DIMENSION(klon,klev+2) :: PFRAC_UP ! updraft fraction 377 REAL, DIMENSION(klon,klev+2) :: PEMF ! updraft mass flux 378 REAL, DIMENSION(klon,klev+2) :: PDETR ! updraft detrainment 379 REAL, DIMENSION(klon,klev+2) :: PENTR ! updraft entrainment 380 INTEGER,DIMENSION(klon) ::IKLCL,IKETL,IKCTL ! level of LCL,ETL and CTL 381 ! Values after saturation adjustement 382 REAL, DIMENSION(klon,klev) :: t_adj ! Adjusted temperature 383 REAL, DIMENSION(klon,klev) :: qv_adj, ql_adj, qi_adj, qr_adj, qs_adj, qg_adj !specific contents after adjustement 384 ! 60 385 real :: temp_newton(klon,klev) 61 386 integer :: k 62 387 logical, save :: first=.true. 63 388 !$OMP THREADPRIVATE(first) 64 65 real,save :: rg=9.81 66 !$OMP THREADPRIVATE(rg) 389 !real,save :: rg=9.81 390 !!$OMP THREADPRIVATE(rg) 67 391 68 392 ! For I/Os 69 393 integer :: itau0 70 real :: zjulian 394 395 real, dimension(nlon) :: swdnsrf,lwdnsrf,lwupsrf,sensible,latent,qsats,capcal,fluxgrd, fluxu,fluxv,fluxt,fluxq 396 397 integer, parameter :: nsoil=11 398 real, dimension(nsoil) :: zout 399 real, dimension(nlon,nlev) :: tsoil_out 400 401 integer :: isoil 402 real,dimension(:), allocatable, save :: tsrf 403 real, dimension(:,:), allocatable, save :: tsoil 404 !$OMP THREADPRIVATE(tsrf,tsoil) 405 406 real :: thermal_inertia0,z0m0,z0h0,z0q0,betasoil 407 408 !AI 409 !real, save :: rpi, stefan, karman, fsw0, timesec 410 !!$OMP THREADPRIVATE(rpi, stefan, karman, fsw0, timesec) 411 real, save :: stefan, fsw0, timesec 412 !$OMP THREADPRIVATE(stefan, fsw0, timesec) 413 integer, save :: iecri=1,iflag_surface=0 414 !$OMP THREADPRIVATE(iecri,iflag_surface) 415 real,save :: zjulian 416 !$OMP THREADPRIVATE(zjulian) 417 ! AI 418 !integer, save :: itap=0 419 !!$OMP THREADPRIVATE(iecri,iflag_surface,rpi, stefan, karman, fsw0, timesec,tsrf, tsoil,itap) 71 420 72 421 … … 76 425 77 426 print*,'Debut physiqex',debut 427 ! AI 428 pdtphys=pdtphys_ 429 CALL update_time(pdtphys) 430 phys_tstep=NINT(pdtphys) 431 78 432 ! initializations 79 433 if (debut) then ! Things to do only for the first call to physics 80 print*,'Debut physiqex IN' 81 82 ! load initial conditions for physics (including the grid) 83 call phys_state_var_init(1) ! some initializations, required before calling phyetat0 84 call phyetat0("startphy.nc", clesphy0, tabcntr0) 85 86 ! Initialize outputs: 87 itau0=0 88 ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 89 !CALL ymds2ju(annee0, month, dayref, hour, zjulian) 90 call ymds2ju(1979, 1, 1, 0.0, zjulian) 91 434 print*,'Debut physiqex IN' 435 print*,'NOUVEAU PHYSIQEX 222222222222222222222' 436 CALL suphel ! initialiser constantes et parametres phys. 437 ! Read parameters and flags in files .def 438 CALL conf_phys(ok_journe, ok_mensuel, & 439 ok_instan, ok_hf, & 440 ok_LES, & 441 callstats, & 442 solarlong0,seuil_inversion, & 443 fact_cldcon, facttemps,ok_newmicro,iflag_radia, & 444 iflag_cld_th,ratqsbas,ratqshaut,tau_ratqs, & 445 ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, & 446 chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, & 447 flag_bc_internal_mixture, bl95_b0, bl95_b1, & 448 read_climoz, alp_offset) 449 450 ! IF (.NOT. create_etat0_limit) 451 CALL init_limit_read(days_elapsed) 452 453 ! load initial conditions for physics (including the grid) 454 ! AI 455 !call phys_state_var_init(1) ! some initializations, required before calling phyetat0 456 call phys_state_var_init(read_climoz) ! some initializations, required before calling phyetat0 457 call phyetat0("startphy.nc", clesphy0, tabcntr0) 458 CALL phys_output_var_init 459 460 ! AI 461 ! Init increments et pas 462 !itap = 0 463 itaprad = 0 464 lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours 465 WRITE(*,*)'La frequence de lecture surface est de ', & 466 lmt_pas 467 IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN 468 radpas = NINT( 86400./phys_tstep)/nbapp_rad 469 ELSE 470 WRITE(*,*) 'le nombre de pas de temps physique doit etre un ', & 471 'multiple de nbapp_rad' 472 WRITE(*,*) 'changer nbapp_rad ou alors commenter ce test ', & 473 'mais 1+1<>2' 474 abort_message='nbre de pas de temps physique n est pas multiple ' & 475 // 'de nbapp_rad' 476 CALL abort_physic(modname,abort_message,1) 477 ENDIF 478 479 ! Initialize outputs: 480 itau0=0 481 ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 482 !CALL ymds2ju(annee0, month, dayref, hour, zjulian) 483 !call ymds2ju(1979, 1, 1, 0.0, zjulian) 484 !print*,'zjulian=',zjulian 485 486 ZDZMIN=MINVAL((pphi(:,2:) - pphi(:,1:klev-1))/9.81) 487 CALL INIT_PHYEX(pdtphys, ZDZMIN, PHYEX) 488 CALL FILL_DIMPHYEX(KLON, KLEV, D) 489 490 ! 491 ! Variables saved 492 ALLOCATE(PTKEM(klon,klev+2)) 493 ALLOCATE(PSIGS(klon,klev+2)) 494 ALLOCATE(PCF_MF(klon,klev+2)) 495 ALLOCATE(PRC_MF(klon,klev+2)) 496 ALLOCATE(PRI_MF(klon,klev+2)) 497 ALLOCATE(ZQR(klon, klev)) 498 ALLOCATE(ZQS(klon, klev)) 499 ALLOCATE(ZQG(klon, klev)) 500 PSIGS=0. 501 PCF_MF=0. 502 PRC_MF=0. 503 PRI_MF=0. 504 ZQR=0. 505 ZQS=0. 506 ZQG=0. 507 PTKEM(:,:) = PHYEX%TURBN%XTKEMIN ! TODO: init from TKE at stationnary state 508 509 ! ---------------------------------------------------------------- 510 ! Initialisation du sol generique 511 ! ---------------------------------------------------------------- 512 allocate(tsrf(nlon)) 513 allocate(tsoil(nlon,nsoil)) 514 thermal_inertia0=2000. 515 z0m0=0.05 516 z0h0=0.05 517 z0q0=0.05 518 betasoil=0.2 519 call soil_ini(klon,nsoil,PHYEX%CST%XCPD,PHYEX%CST%XLVTT,thermal_inertia0,z0m0,z0h0,z0q0,betasoil) 520 ! 92 521 #ifndef CPP_IOIPSL_NO_OUTPUT 93 522 ! Initialize IOIPSL output file 94 523 #endif 524 ! AI Partie Surface depplacee 525 ! compute tendencies to return to the dynamics: 526 ! "friction" on the first layer 527 ! tsrf(1:nlon)=t(1:nlon,1) 528 ! do isoil=1,nsoil 529 ! tsoil(1:nlon,isoil)=tsrf(1:nlon) 530 ! !tsoil(1:nlon,isoil)=302. 531 ! enddo 532 ! rpi=2.*asin(1.) 533 ! timesec=5*3600 534 ! stefan=5.67e-8 535 ! fsw0=900. 536 ! call getin_p('iecri',iecri) 537 ! call getin_p('iflag_surface',iflag_surface) 95 538 96 539 endif ! of if (debut) 97 540 541 ! AI 542 ! Increment 543 !itap = itap + 1 98 544 !------------------------------------------------------------ 99 545 ! Initialisations a chaque pas de temps 100 546 !------------------------------------------------------------ 101 102 547 103 548 ! set all tendencies to zero … … 106 551 d_t(1:klon,1:klev)=0. 107 552 d_qx(1:klon,1:klev,1:nqtot)=0. 553 d_qr(1:klon,1:klev)=0. 554 d_qs(1:klon,1:klev)=0. 555 d_qg(1:klon,1:klev)=0. 108 556 d_ps(1:klon)=0. 109 110 !------------------------------------------------------------ 111 ! Calculs 112 !------------------------------------------------------------ 113 557 d_tke(1:klon,1:klev)=0. 558 ! 559 ! AI attention 560 d_t_swr = 0. 561 d_t_sw0 = 0. 562 d_t_lwr = 0. 563 d_t_lw0 = 0. 564 ! 565 ZDXX(:,:) = 0. 566 ZDYY(:,:) = 0. 567 ZDZX(:,:) = 0. 568 ZDZY(:,:) = 0. 569 ZDIRCOSXW(:) = 1. 570 ZDIRCOSYW(:) = 1. 571 ZDIRCOSZW(:) = 1. 572 ZCOSSLOPE(:) = 0. 573 ZSINSLOPE(:) = 1. 574 PHGRAD(:,:,:) = 0. 575 ZBL_DEPTH(:) = 0. ! needed only with LRMC01 key (correction in the surface boundary layer) 576 ZSBL_DEPTH(:) = 0. 577 ZCEI(:,:) = 0. ! needed only if HTURBLEN_CL /= 'NONE' modification of mixing lengh inside clouds 578 ZSVT(:,:,:) = 0. 579 PWT(:,:) = 0. 580 ZUT(:,2:klev+1) = u(:,:) 581 ZVT(:,2:klev+1) = v(:,:) 582 ! 583 !------------------------------------------------------------ 584 ! Conversions and extra levels 585 !------------------------------------------------------------ 586 !TODO check in Meso-NH how values are extrapolated outside of the physical domain 587 zqt(:,2:klev+1) = qx(:,:,1) + qx(:,:,2) + qx(:,:,3) 588 zqt(:,1)=0. 589 zqt(:,klev+2)=0. 590 zqdm(:,:)=1.-zqt(:,:) !equal to 1/(1+rt) 591 ZRX(:,2:klev+1,1) = qx(:,:,1) / zqdm(:,2:klev+1) 592 ZRX(:,2:klev+1,2) = qx(:,:,2) / zqdm(:,2:klev+1) 593 ZRX(:,2:klev+1,4) = qx(:,:,3) / zqdm(:,2:klev+1) 594 ZRX(:,2:klev+1,3) = ZQR(:,:) 595 ZRX(:,2:klev+1,5) = ZQS(:,:) 596 ZRX(:,2:klev+1,6) = ZQG(:,:) 597 DO JRR=1,KRR 598 CALL VERTICAL_EXTEND(ZRX(:,:,JRR),klev) 599 END DO 600 ! 601 ZEXN(:,2:klev+1) = (pplay(:,:) / PHYEX%CST%XP00) ** (PHYEX%CST%XRD/PHYEX%CST%XCPD) 602 ZTHETA(:,2:klev+1) = t(:,:) / ZEXN(:,2:klev+1) 603 CALL VERTICAL_EXTEND(ZEXN,klev) 604 CALL VERTICAL_EXTEND(ZTHETA,klev) 605 606 !TODO check in Meso-NH how zz_mass and zz_flux are initialized outside of the physical domain 607 zs(:) = pphis(:)/PHYEX%CST%XG 608 zz_mass(:,2:klev+1) = pphi(:,:) / PHYEX%CST%XG 609 zz_mass(:,1) = 2*zs-zz_mass(:,2) 610 zz_mass(:,klev+2)=2.*zz_mass(:,klev+1)-zz_mass(:,klev) 611 612 do k=2, klev+2 613 zz_flux(:,k)=(zz_mass(:,k-1)+zz_mass(:,k))/2. 614 enddo 615 zz_flux(:,1)=2*zz_mass(:,1)-zz_flux(:,2) 616 617 !zdzf is the distance between two consecutive flux points (expressed on mass points) 618 do k=1,klev+1 619 zdzf(:,k)=zz_flux(:,k+1)-zz_flux(:,k) 620 enddo 621 zdzf(:,klev+2)=(zz_mass(:,klev+2)-zz_flux(:,klev+2))*2. 622 623 !zdzm distance between two consecutive mass points (expressed on flux points) 624 do k=2,klev+2 625 zdzm(:,k)=zz_mass(:,k)-zz_mass(:,k-1) 626 enddo 627 zdzm(:,1)=(zz_mass(:,1)-zz_flux(:,1))*2. 628 629 ZPABST(:,2:klev+1) = pplay(:,:) 630 ZRHOD(:,2:klev+1)=ZPABST(:,2:klev+1)/(t*(PHYEX%CST%XRD+ZRX(:,2:klev+1,1)*PHYEX%CST%XRV)) 631 DO k=2,klev+1 632 PRHODJ(:,k) = ZRHOD(:,k) * (zdzf(:,k)*cell_area(:)) 633 END DO 634 PTHVREF(:,:) = ZTHETA(:,:) * (1. + PHYEX%CST%XRV/PHYEX%CST%XRD * ZRX(:,:,1)) * ZQDM(:,:) 635 636 CALL VERTICAL_EXTEND(ZPABST,klev) 637 CALL VERTICAL_EXTEND(PRHODJ,klev) 638 CALL VERTICAL_EXTEND(ZRHOD,klev) 639 CALL VERTICAL_EXTEND(ZUT,klev) 640 CALL VERTICAL_EXTEND(ZVT,klev) 641 642 643 !------------------------------------------------------------ 644 ! Tendencies 645 !------------------------------------------------------------ 646 !For Meso-NH, initialia values for the tendencies are filled with 647 !a pseudo-tendecy computed by dividing the state variable by the time step 648 !This mechanism enables the possibility for the different parametrisations 649 !to guess the value at the end of the time step. 650 !For the wind components, we could do the same way but it is not needed 651 !as the parametrisations don't use the S varaible to guess the futur value of the wind. 652 ZRXS(:,:,:) = ZRX(:,:,:)/pdtphys 653 ZTHETAS(:,:)=ZTHETA(:,:)/pdtphys 654 ZTKES(:,:)=PTKEM(:,:)/pdtphys 655 !To compute the actual tendency, we save the initial values of these variables 656 ZRXS0(:,:,:) = ZRXS(:,:,:) 657 ZTHETAS0(:,:)=ZTHETAS(:,:) 658 ZTKES0(:,:)=ZTKES(:,:) 659 !------------------------------------------------------------ 660 ! Adjustment 661 !------------------------------------------------------------ 662 ! 663 ZSRC(:,:) = 0. 664 ZSIGQSAT=PHYEX%NEBN%VSIGQSAT 665 CALL ICE_ADJUST (D, PHYEX%CST, PHYEX%RAIN_ICE_PARAMN, PHYEX%NEBN, PHYEX%TURBN, PHYEX%PARAM_ICEN, & 666 &PHYEX%MISC%TBUCONF, KRR, & 667 &'ADJU', & 668 &pdtphys, ZSIGQSAT, & 669 &PRHODJ, ZEXN, ZRHOD, PSIGS, .FALSE., zmfconv, & 670 &ZPABST, ZZ_MASS, & 671 &ZEXN, PCF_MF, PRC_MF, PRI_MF, & 672 &ZICLDFR, ZWCLDFR, ZSSIO, ZSSIU, ZIFR, & 673 &ZRX(:,:,1), ZRX(:,:,2), ZRXS(:,:,1), ZRXS(:,:,2), ZTHETA, ZTHETAS, & 674 &PHYEX%MISC%COMPUTE_SRC, ZSRC, ZCLDFR, & 675 &ZRX(:,:,3), ZRX(:,:,4), ZRXS(:,:,4), ZRX(:,:,5), ZRX(:,:,6), & 676 &PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET, & 677 &ZICE_CLD_WGT, & 678 &PHLC_HRC=ZHLC_HRC, PHLC_HCF=ZHLC_HCF, PHLI_HRI=ZHLI_HRI, PHLI_HCF=ZHLI_HCF ) 679 ! 680 !Variables are updated with their adjusted values (to be used by the other parametrisations) 681 ZTHETA(:,:)=ZTHETAS(:,:)*pdtphys 682 ZRX(:,:,:)=ZRXS(:,:,:)*pdtphys 683 t_adj=ZTHETA(:,2:klev+1)*ZEXN(:,2:klev+1) 684 qv_adj=ZRX(:,2:klev+1,1)*zqdm(:,2:klev+1) 685 ql_adj=ZRX(:,2:klev+1,2)*zqdm(:,2:klev+1) 686 qi_adj=ZRX(:,2:klev+1,4)*zqdm(:,2:klev+1) 687 qr_adj=ZRX(:,2:klev+1,3)*zqdm(:,2:klev+1) 688 qs_adj=ZRX(:,2:klev+1,5)*zqdm(:,2:klev+1) 689 qg_adj=ZRX(:,2:klev+1,6)*zqdm(:,2:klev+1) 690 ZRHOD(:,2:klev+1)=ZPABST(:,2:klev+1)/(t*(PHYEX%CST%XRD+ZRX(:,2:klev+1,1)*PHYEX%CST%XRV)) 691 CALL VERTICAL_EXTEND(ZRHOD,klev) 692 ! 693 !------------------------------------------------------------ 694 ! Radiation 695 !------------------------------------------------------------ 696 ! 697 !AI a nettoyer apres 698 ! calcul de qsat Fred mai 2024 699 zqsat(:,:) = EXP( XALPW - XBETAW / t(:,:) - XGAMW * ALOG( t(:,:) ) ) 700 zqsat(:,:) = XRD / XRV * zqsat(:,:) / ( pplay(:,:) - zqsat(:,:)) 701 ! 702 if (iflag_radia.ge.1) then 703 ! 704 !ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors 705 !Cloud fraction is available in the ZCLDFR variable 706 ! 707 ! AI attention : dans un 1er temps read_climoz=0 708 ! Update ozone if day change 709 !print*,'solarlong0, itap, lmt_pas, days_elapsed : ', solarlong0, itap, lmt_pas, days_elapsed 710 !stop 711 IF (MOD(itap-1,lmt_pas) == 0) THEN 712 ! IF (read_climoz <= 0) THEN 713 ! Once per day, update ozone from Royer: 714 IF (solarlong0<-999.) then 715 ! Generic case with evolvoing season 716 zzzo=real(days_elapsed+1) 717 ELSE IF (abs(solarlong0-1000.)<1.e-4) then 718 ! Particular case with annual mean insolation 719 zzzo=real(90) ! could be revisited 720 IF (read_climoz/=-1) THEN 721 abort_message ='read_climoz=-1 is recommended when ' & 722 // 'solarlong0=1000.' 723 CALL abort_physic (modname,abort_message,1) 724 ENDIF 725 ELSE 726 ! Case where the season is imposed with solarlong0 727 zzzo=real(90) ! could be revisited 728 ENDIF 729 zzzo=real(days_elapsed+1) 730 wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzzo) 731 !print*,'wo = ', wo 732 ! ELSE 733 ! !--- ro3i = elapsed days number since current year 1st january, 0h 734 ! ro3i=days_elapsed+jh_cur-jh_1jan 735 ! !--- scaling for old style files (360 records) 736 ! IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len 737 ! IF(adjust_tropopause) THEN 738 ! CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), & 739 ! ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), & 740 ! time_climoz , longitude_deg, latitude_deg, & 741 ! dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)) 742 ! ELSE 743 ! CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), & 744 ! ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), & 745 ! time_climoz ) 746 ! ENDIF 747 ! ! Convert from mole fraction of ozone to column density of ozone in a 748 ! ! cell, in kDU: 749 ! FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd & 750 ! * zmasse / dobson_u / 1e3 751 ! ! (By regridding ozone values for LMDZ only once a day, we 752 ! ! have already neglected the variation of pressure in one 753 ! ! day. So do not recompute "wo" at each time step even if 754 ! ! "zmasse" changes a little.) 755 ! ENDIF 756 ENDIF 757 !========================================================================= 758 ! Calculs de l'orbite. 759 ! Necessaires pour le rayonnement et la surface (calcul de l'albedo). 760 ! doit donc etre plac\'e avant radlwsw et pbl_surface 761 !CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq) 762 CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq) 763 print*,'mth_eq, day_eq, jD_eq, jH_cur =', mth_eq, day_eq, jD_eq, jH_cur 764 day_since_equinox = (jD_cur + jH_cur) - jD_eq 765 !choix entre calcul de la longitude solaire vraie ou valeur fixee = solarlong0 766 IF (solarlong0<-999.) THEN 767 IF (new_orbit) THEN 768 ! calcul selon la routine utilisee pour les planetes 769 CALL solarlong(day_since_equinox, zlongi, dist) 770 print*,'day_since_equinox, zlongi, dist : ', day_since_equinox, zlongi, dist 771 ELSE 772 ! calcul selon la routine utilisee pour l'AR4 773 CALL orbite(REAL(days_elapsed+1),zlongi,dist) 774 ENDIF 775 ELSE 776 zlongi=solarlong0 ! longitude solaire vraie 777 dist=1. ! distance au soleil / moyenne 778 ENDIF 779 !IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist 780 781 ! AI a refaire en definissant les itap, itaprad 782 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 783 ! ============================ 784 ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur 785 ! l'annee a partir d'une formule analytique. 786 ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et 787 ! non nul aux poles. 788 IF (abs(solarlong0-1000.)<1.e-4) THEN 789 CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, & 790 latitude_deg,longitude_deg,rmu0,fract) 791 swradcorr(:) = 1.0 792 ! JrNt(:) = 1.0 793 zrmu0(:) = rmu0(:) 794 ELSE 795 ! recode par Olivier Boucher en sept 2015 796 SELECT CASE (iflag_cycle_diurne) 797 CASE(0) 798 ! Sans cycle diurne 799 CALL angle(zlongi, latitude_deg, fract, rmu0) 800 swradcorr = 1.0 801 ! JrNt = 1.0 802 zrmu0 = rmu0 803 CASE(1) 804 ! Avec cycle diurne sans application des poids 805 ! bit comparable a l ancienne formulation cycle_diurne=true 806 ! on integre entre gmtime et gmtime+radpas 807 zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s) 808 CALL zenang(zlongi,jH_cur,0.0,zdtime, & 809 latitude_deg,longitude_deg,rmu0,fract) 810 zrmu0 = rmu0 811 swradcorr = 1.0 812 ! Calcul du flag jour-nuit 813 ! JrNt = 0.0 814 !WHERE (fract.GT.0.0) JrNt = 1.0 815 CASE(2) 816 ! Avec cycle diurne sans application des poids 817 ! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1) 818 ! Comme cette routine est appele a tous les pas de temps de 819 ! la physique meme si le rayonnement n'est pas appele je 820 ! remonte en arriere les radpas-1 pas de temps 821 ! suivant. Petite ruse avec MOD pour prendre en compte le 822 ! premier pas de temps de la physique pendant lequel 823 ! itaprad=0 824 zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1) 825 zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1) 826 CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, & 827 latitude_deg,longitude_deg,rmu0,fract) 828 ! 829 ! Calcul des poids 830 ! 831 zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le 832 zdtime2=0.0 !--pas de temps de la physique qui se termine 833 CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, & 834 latitude_deg,longitude_deg,zrmu0,zfract) 835 swradcorr = 0.0 836 WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) & 837 swradcorr=zfract/fract*zrmu0/rmu0 838 ! Calcul du flag jour-nuit 839 !JrNt = 0.0 840 !WHERE (zfract.GT.0.0) JrNt = 1.0 841 END SELECT 842 ENDIF 843 !print*,'fract = ',fract 844 !sza_o = ACOS (rmu0) *180./pi 845 846 ! Calcul de l'ensoleillement : 847 IF (MOD(itaprad,radpas).EQ.0) THEN 848 namelist_ecrad_file='namelist_ecrad' 849 ! 850 ! Q0 gestion des pas de temps : 851 ! itap, itaprad ????? 852 ! 853 ! Q1 champs lies a la surface ? sorties appel au modele de surface 854 ! zxtsol temp au sol : dvrait etre recuperee de la partie surface 855 ! albsol_dir, albsol_dif aussi 856 albsol_dir = 1. ! AI attention 857 albsol_dif = 1. ! AI attention 858 zxtsol(:) = t(:,1) ! AI attention 859 ! 860 ! Q2 ozone ? lmdz ok a completer 861 ! wo : pour l'instant appel ozonecm 862 ! 863 ! Q3 orbite ? lmdz ok a comleter 864 ! dist, rmu0 et fract 865 ! pour l'instant cas solarlong0<-999. 866 ! comment introduire les pas de temps rayonnement 867 ! 868 ! Q3 proprietes optiques des nuages ? call cloud_optics_prop.F90 ? non necessaires pour ecrad 869 ! epaisseur optique, emissivite : cldfrarad, cldemirad, cldtaurad , non necessaire pour ecrad 870 ! rayons effectifs : ref_liq, ref_ice ????? 871 ! cloud liq/Ice water mixing contenent (kg/kg) ???? 872 ! les rayons effectifs seront calcules en utilisant la routine de l'IFS 873 ! en fonct des contenus en eau, ..... 874 cldemirad = 1. 875 cldtaurad = 5. 876 cldtaupirad = 5. 877 ref_liq = 1. 878 ref_ice = 1. 879 ref_liq_pi = 1. 880 ref_ice_pi = 1. 881 ! 882 ! Q4 aerosols ? 883 ! simplifier l'ecriture : mettre tous les flags et parametres dans aero_mod.F90 ? 884 ! ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, & 885 ! flag_aerosol, flag_aerosol_strat, flag_aer_feedback, & 886 ! tau_aero, piz_aero, cg_aero, & 887 ! tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & 888 ! tau_aero_lw_rrtm, & 889 ! 890 ! Q5 thermo : zqsat 891 ! zqsat specific humidity at liq saturation => 892 ! calcul de Fred ou a calculer dans ecrad enf fct de t et p 893 ! flwc, fiwc : qx( 2) et qx( 3) 894 895 print*,'Entree radlwsw a itap : ', itap 896 CALL radlwsw & 897 (debut, dist, rmu0, fract, & 898 & paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, & 899 & t,qx(:,:,1),wo, & 900 & ZCLDFR, cldemirad, cldtaurad, & 901 & ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, & 902 & flag_aerosol, flag_aerosol_strat, flag_aer_feedback, & 903 & tau_aero, piz_aero, cg_aero, & 904 & tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & 905 ! Rajoute par OB pour RRTM 906 & tau_aero_lw_rrtm, & 907 & cldtaupirad, m_allaer, & 908 ! 909 ! & zqsat, flwc, fiwc, & 910 & zqsat, qx(:,:,2), qx(:,:,2), & 911 & ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, & 912 & namelist_ecrad_file, & 913 & heat,heat0,cool,cool0,albpla, & 914 & heat_volc,cool_volc, & 915 & topsw,toplw,solsw,solswfdiff,sollw, & 916 & sollwdown, & 917 & topsw0,toplw0,solsw0,sollw0, & 918 & lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, & 919 & swdnc0, swdn0, swdn, swupc0, swup0, swup, & 920 & topswad_aero, solswad_aero, & 921 & topswai_aero, solswai_aero, & 922 & topswad0_aero, solswad0_aero, & 923 & topsw_aero, topsw0_aero, & 924 & solsw_aero, solsw0_aero, & 925 & topswcf_aero, solswcf_aero, & 926 !-C. Kleinschmitt for LW diagnostics 927 & toplwad_aero, sollwad_aero,& 928 & toplwai_aero, sollwai_aero, & 929 & toplwad0_aero, sollwad0_aero,& 930 !-end 931 & ZLWFT0_i, ZFLDN0, ZFLUP0, & 932 & ZSWFT0_i, ZFSDN0, ZFSUP0, & 933 & cloud_cover_sw) 934 itaprad = 0 935 ENDIF ! IF (MOD(itaprad,radpas).EQ.0) 936 itaprad = itaprad + 1 937 938 IF (iflag_radia.eq.0) THEN 939 heat=0. 940 cool=0. 941 sollw=0. ! MPL 01032011 942 solsw=0. 943 radsol=0. 944 swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars 945 swup0=0. 946 lwup=0. 947 lwup0=0. 948 lwdn=0. 949 lwdn0=0. 950 ENDIF 951 ! Ajouter la tendance des rayonnements (tous les pas) 952 ! avec une correction pour le cycle diurne dans le SW 953 ! 954 DO k=1, klev 955 ! AI attention 956 !d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY 957 !d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY 958 !d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY 959 !d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY 960 d_t_swr(:,k)=swradcorr(:)*heat(:,k)/RDAY 961 d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)/RDAY 962 d_t_lwr(:,k)=-cool(:,k)/RDAY 963 d_t_lw0(:,k)=-cool0(:,k)/RDAY 964 ENDDO 965 966 !t_seri(:,:) = t_seri(:,:) + d_t_swr(:,:) + d_t_lwr(:,:) 967 d_t = d_t + d_t_swr + d_t_lwr 968 969 !CALL writefield_phy('d_t_swr',d_t_swr,klev) 970 !CALL writefield_phy('d_t_lwr',d_t_lwr,klev) 971 !CALL writefield_phy('temp',t,klev) 972 973 endif ! iflag_radia a nettoyer 974 ! 975 ! AI passer peut-etre la surface avant radiation 976 ! albedo, temp surf, ... pour le rayonnement 977 ! 978 !ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors 979 !Cloud fraction is available in the ZCLDFR variable 980 981 ! 982 !------------------------------------------------------------ 983 ! Surface 984 !------------------------------------------------------------ 985 ! 986 ! AI 114 987 ! compute tendencies to return to the dynamics: 115 988 ! "friction" on the first layer 116 d_u(1:klon,1)=-u(1:klon,1)/86400. 117 d_v(1:klon,1)=-v(1:klon,1)/86400. 118 ! newtonian relaxation towards temp_newton() 119 do k=1,klev 120 temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 121 d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5 122 enddo 123 124 989 if (debut) then 990 tsrf(1:nlon)=t(1:nlon,1) 991 do isoil=1,nsoil 992 tsoil(1:nlon,isoil)=tsrf(1:nlon) 993 !tsoil(1:nlon,isoil)=302. 994 enddo 995 rpi=2.*asin(1.) 996 timesec=5*3600 997 stefan=5.67e-8 998 fsw0=900. 999 call getin_p('iecri',iecri) 1000 call getin_p('iflag_surface',iflag_surface) 1001 endif 1002 1003 timesec=timesec+pdtphys 1004 1005 ! Computation of qsats. Done even it is not used for ouptuts 1006 ! AI 1007 !qsats(:) = EXP( XALPW - XBETAW / t(:,1) - XGAMW * ALOG( t(:,1) ) ) 1008 !qsats(:) = XRD / XRV * qsats(:) / ( pplay(:,1) - qsats(:) ) 1009 qsats(:) = zqsat(:,1) 1010 1011 if (iflag_surface == 0 ) then 1012 ! Latent and sensible heat fluxes imposed in 1D (exemple Arm cumulus) 1013 ! radiative fluxes set to 0. 1014 ! capcal and fluxgrd=0. replace generic_soil and allow to compute a surface temperature 1015 ! A case close to Rico can be obtained by imposing fsens=-5E-3 et flat=-6E-5 1016 capcal(:)=1.e6 1017 fluxgrd(:)=0. 1018 sensible(:)=-fsens 1019 latent(:)=-flat 1020 swdnsrf(:)=0. 1021 lwdnsrf(:)=0. 1022 fluxt(1:klon)=sensible/PHYEX%CST%XCPD 1023 fluxq(1:klon)=latent(:)/PHYEX%CST%XLVTT 1024 fluxu(1:klon)=-0.01*u(1:klon,1) 1025 fluxv(1:klon)=-0.01*v(1:klon,1) 1026 else 1027 ! AI 1028 if (iflag_radia.eq.0) then 1029 ! Using the generic soil model 1030 if ( nlon == 1 ) then 1031 ! Imposed radiative fluxes to simulate a 1D case close to Arm Cumulus 1032 swdnsrf(:)=fsw0*max(cos(2*rpi*(timesec-43200.)/(1.15*86400.)),0.) 1033 lwdnsrf(:)=400. 1034 else 1035 ! Model forced in latitude by an idelaized radiative flux that depends on latitude 1036 ! only. Used for 3D tests without radiation. 1037 swdnsrf(:)=600*cos(latitude(:))+100 1038 lwdnsrf(:)=0. 1039 endif 1040 else 1041 swdnsrf(:)=solsw(:) 1042 lwdnsrf(:)=sollwdown(:) 1043 endif 1044 ! updating the soil temperature (below the surface) from the surface temperature 1045 ! en returning the surface conductive flux and a heat capacity (accounting for the 1046 ! sentivity of the conductive flux to the surface temperature. 1047 call soil_conduct(klon,nsoil,debut,pdtphys,tsrf,tsoil,capcal,fluxgrd,zout) 1048 1049 ! tsoil_out = tsoil with the vertical dimension of atmosphere outputs 1050 1051 ! Computing the turbulent fluxes Fqx = rho w'x' = rho (kappa/ln(z/z0x))^2 [ x_surface - x ]. Upward 1052 call soil_fluxes(klon,tsrf,qsats,pphi(:,1)/rg,ZRHOD(:,2),u(:,1),v(:,1),t(:,1),qx(:,1,1),fluxu,fluxv,fluxt,fluxq) 1053 sensible(:)=fluxt(:)*PHYEX%CST%XCPD 1054 latent(:)=fluxq(:)*PHYEX%CST%XLVTT 1055 endif 1056 tsoil_out=0. ; tsoil_out(:,1:min(nsoil,nlev))=tsoil(:,1:min(nsoil,nlev)) 1057 1058 ! Computing the surface temperature 1059 lwupsrf=stefan*tsrf(:)**4 1060 tsrf(:)=tsrf(:)+pdtphys*(-latent-sensible+lwdnsrf+swdnsrf-lwupsrf+fluxgrd)/capcal(:) 1061 1062 PSFTH(:) = fluxt(:) 1063 PSFRV(:) = fluxq(:) 1064 PSFU(:) = fluxu(:) 1065 PSFV(:) = fluxv(:) 1066 PSFSV(:,:) = 0. 1067 ! 1068 !------------------------------------------------------------ 1069 ! Shallow convection 1070 !------------------------------------------------------------ 1071 ! 1072 CALL SHALLOW_MF(D, PHYEX%CST, PHYEX%NEBN, PHYEX%PARAM_MFSHALLN, PHYEX%TURBN, PHYEX%CSTURB, & 1073 &KRR=KRR, KRRL=KRRL, KRRI=KRRI, KSV=KSV, & 1074 &ONOMIXLG=PHYEX%MISC%ONOMIXLG,KSV_LGBEG=PHYEX%MISC%KSV_LGBEG,KSV_LGEND=PHYEX%MISC%KSV_LGEND, & 1075 &PTSTEP=pdtphys, & 1076 &PDZZ=zdzm(:,:),PZZ=zz_mass(:,:), & 1077 &PRHODJ=PRHODJ(:,:),PRHODREF=ZRHOD(:,:), & 1078 &PPABSM=ZPABST(:,:),PEXNM=ZEXN(:,:), & 1079 &PSFTH=PSFTH(:),PSFRV=PSFRV(:), & 1080 &PTHM=ZTHETA(:,:),PRM=ZRX(:,:,:),PUM=ZUT(:,:),PVM=ZVT(:,:), & 1081 &PTKEM=PTKEM(:,:),PSVM=ZSVT(:,:,:), & 1082 &PDUDT_MF=PDUDT_MF(:,:),PDVDT_MF=PDVDT_MF(:,:), & 1083 &PDTHLDT_MF=PDTHLDT_MF(:,:),PDRTDT_MF=PDRTDT_MF(:,:),PDSVDT_MF=PDSVDT_MF(:,:,:), & 1084 &PSIGMF=PSIGMF(:,:),PRC_MF=PRC_MF(:,:),PRI_MF=PRI_MF(:,:),PCF_MF=PCF_MF(:,:), & 1085 &PFLXZTHVMF=ZFLXZTHVMF(:,:), & 1086 &PFLXZTHMF=ZFLXZTHMF(:,:),PFLXZRMF=ZFLXZRMF(:,:),PFLXZUMF=ZFLXZUMF(:,:),PFLXZVMF=ZFLXZVMF(:,:), & 1087 &PTHL_UP=PTHL_UP(:,:),PRT_UP=PRT_UP(:,:),PRV_UP=PRV_UP(:,:), & 1088 &PRC_UP=PRC_UP(:,:),PRI_UP=PRI_UP(:,:), & 1089 &PU_UP=PU_UP(:,:), PV_UP=PV_UP(:,:), PTHV_UP=PTHV_UP(:,:), PW_UP=PW_UP(:,:), & 1090 &PFRAC_UP=PFRAC_UP(:,:),PEMF=PEMF(:,:),PDETR=PDETR(:,:),PENTR=PENTR(:,:), & 1091 &KKLCL=IKLCL(:),KKETL=IKETL(:),KKCTL=IKCTL(:),PDX=1000.0,PDY=1000.0,KBUDGETS=PHYEX%MISC%NBUDGET ) 1092 1093 ! Add tendencies of shallow to total physics tendency 1094 d_u(:,1:klev) = d_u(:,1:klev) + PDUDT_MF(:,2:klev+1) 1095 d_v(:,1:klev) = d_v(:,1:klev) + PDVDT_MF(:,2:klev+1) 1096 ZRXS(:,:,1)=ZRXS(:,:,1)+PDRTDT_MF(:,:) 1097 ZTHETAS(:,:)=ZTHETAS(:,:)+PDTHLDT_MF(:,:) 1098 ! TODO add SV tendencies 1099 ! 1100 !------------------------------------------------------------ 1101 ! Turbulence 1102 !------------------------------------------------------------ 1103 ! out tendencies 1104 ZRUS(:,:) = 0. 1105 ZRVS(:,:) = 0. 1106 ZRWS(:,:) = 0. 1107 ZRSVS(:,:,:) = 0. 1108 ZRTKES(:,:) = ZTKES(:,:) * PRHODJ(:,:) 1109 DO JRR=1, KRR 1110 ZRRS(:,:,JRR) = ZRXS(:,:,JRR) * PRHODJ(:,:) 1111 ENDDO 1112 ZRTHS(:,:) = ZTHETAS(:,:) * PRHODJ(:,:) 1113 CALL TURB(PHYEX%CST, PHYEX%CSTURB, PHYEX%MISC%TBUCONF, PHYEX%TURBN, PHYEX%NEBN, D, PHYEX%MISC%TLES, & 1114 & KRR, KRRL, KRRI, PHYEX%MISC%HLBCX, PHYEX%MISC%HLBCY, PHYEX%MISC%KGRADIENTS, PHYEX%MISC%KHALO, & 1115 & PHYEX%TURBN%NTURBSPLIT, PHYEX%TURBN%LCLOUDMODIFLM, KSV, PHYEX%MISC%KSV_LGBEG, PHYEX%MISC%KSV_LGEND, & 1116 & PHYEX%MISC%KSV_LIMA_NR, PHYEX%MISC%KSV_LIMA_NS, PHYEX%MISC%KSV_LIMA_NG, PHYEX%MISC%KSV_LIMA_NH, & 1117 & PHYEX%MISC%O2D, PHYEX%MISC%ONOMIXLG, PHYEX%MISC%OFLAT, PHYEX%MISC%OCOUPLES, & 1118 & PHYEX%MISC%OBLOWSNOW,PHYEX%MISC%OIBM, & 1119 & PHYEX%MISC%OFLYER, PHYEX%MISC%COMPUTE_SRC, PHYEX%MISC%PRSNOW, & 1120 & PHYEX%MISC%OOCEAN, PHYEX%MISC%ODEEPOC, PHYEX%MISC%ODIAG_IN_RUN, & 1121 & PHYEX%TURBN%CTURBLEN_CLOUD, PHYEX%MISC%CMICRO, PHYEX%MISC%CELEC, & 1122 & pdtphys,PHYEX%MISC%ZTFILE, & 1123 & ZDXX(:,:),ZDYY(:,:),zdzm(:,:), & 1124 & ZDZX(:,:),ZDZY(:,:),zz_flux(:,:), & 1125 & ZDIRCOSXW(:),ZDIRCOSYW(:),ZDIRCOSZW(:),ZCOSSLOPE(:),ZSINSLOPE(:), & 1126 & PRHODJ(:,:),PTHVREF(:,:), PHGRAD(:,:,:), zs(:), & 1127 & PSFTH(:),PSFRV(:),PSFSV(:,:),PSFU(:),PSFV(:), & 1128 & ZPABST(:,:),ZUT(:,:),ZVT(:,:),PWT(:,:),PTKEM(:,:),ZSVT(:,:,:),ZSRC(:,:), & 1129 & PLENGTHM(:,:),PLENGTHH(:,:),MFMOIST(:,:), & 1130 & ZBL_DEPTH(:),ZSBL_DEPTH(:), & 1131 & ZCEI(:,:), PHYEX%TURBN%XCEI_MIN, PHYEX%TURBN%XCEI_MAX, PHYEX%TURBN%XCOEF_AMPL_SAT, & 1132 & ZTHETA(:,:),ZRX(:,:,:), & 1133 & ZRUS(:,:),ZRVS(:,:),ZRWS(:,:),ZRTHS(:,:),ZRRS(:,:,:),ZRSVS(:,:,:),ZRTKES(:,:), & 1134 & PSIGS(:,:), & 1135 & ZFLXZTHVMF(:,:),ZWTH(:,:),ZWRC(:,:),ZWSV(:,:,:),ZDP(:,:),ZTP(:,:),ZTDIFF(:,:),ZTDISS(:,:), & 1136 & PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET ) 1137 DO JRR=1, KRR 1138 ZRXS(:,:,JRR) = ZRRS(:,:,JRR) / PRHODJ(:,:) 1139 ENDDO 1140 ZTHETAS(:,:) = ZRTHS(:,:) / PRHODJ(:,:) 1141 ZTKES(:,:) = ZRTKES(:,:) / PRHODJ(:,:) 1142 ! Add tendencies of turb to total physics tendency 1143 d_u(:,1:klev) = d_u(:,1:klev) + ZRUS(:,2:klev+1)/PRHODJ(:,2:klev+1) 1144 d_v(:,1:klev) = d_v(:,1:klev) + ZRVS(:,2:klev+1)/PRHODJ(:,2:klev+1) 1145 IF(PHYEX%PARAM_MFSHALLN%CMF_CLOUD=='STAT') THEN 1146 PSIGS(:,:)=SQRT(PSIGS(:,:)**2 + PSIGMF(:,:)**2) 1147 ENDIF 1148 !------------------------------------------------------------ 1149 ! Microphysics 1150 !------------------------------------------------------------ 1151 ZSEA=1. 1152 ZTOWN=0. 1153 ZCIT=0. 1154 ZTHVREFZIKB=0 1155 CALL RAIN_ICE (D, PHYEX%CST, PHYEX%PARAM_ICEN, PHYEX%RAIN_ICE_PARAMN, PHYEX%RAIN_ICE_DESCRN, & 1156 PHYEX%ELEC_PARAM, PHYEX%ELEC_DESCR, PHYEX%MISC%TBUCONF, & 1157 PHYEX%MISC%OELEC, PHYEX%MISC%OSEDIM_BEARD, & 1158 ZTHVREFZIKB, & 1159 pdtphys, KRR, ZEXN, & 1160 zdzf, PRHODJ, ZRHOD, ZEXN, ZPABST, ZCIT, ZCLDFR, & 1161 ZHLC_HRC, ZHLC_HCF, ZHLI_HRI, ZHLI_HCF, & 1162 ztheta, ZRX(:,:,1), ZRX(:,:,2), ZRX(:,:,3), ZRX(:,:,4), ZRX(:,:,5), & 1163 ZRX(:,:,6), zthetas, ZRXS(:,:,1), ZRXS(:,:,2), ZRXS(:,:,3), ZRXS(:,:,4), ZRXS(:,:,5), ZRXS(:,:,6), & 1164 ZINPRC, ZINPRR, ZEVAP3D, & 1165 ZINPRS, ZINPRG, ZINDEP, ZRAINFR, PSIGS, & 1166 PHYEX%MISC%YLBUDGET, PHYEX%MISC%NBUDGET, & 1167 ZSEA, ZTOWN ) 1168 1169 !------------------------------------------------------------ 1170 ! Tendencies and time evolution (values for next time step) 1171 !------------------------------------------------------------ 1172 ! Tendencies, mixing ratio -> specific 1173 d_qx(:,1:klev,1)=d_qx(:,1:klev,1) + (ZRXS(:,2:klev+1,1)-ZRXS0(:,2:klev+1,1))*ZQDM(:,2:klev+1) 1174 d_qx(:,1:klev,2)=d_qx(:,1:klev,2) + (ZRXS(:,2:klev+1,2)-ZRXS0(:,2:klev+1,2))*ZQDM(:,2:klev+1) 1175 d_qx(:,1:klev,3)=d_qx(:,1:klev,3) + (ZRXS(:,2:klev+1,4)-ZRXS0(:,2:klev+1,4))*ZQDM(:,2:klev+1) 1176 d_qr(:,1:klev)=d_qr(:,1:klev) + (ZRXS(:,2:klev+1,3)-ZRXS0(:,2:klev+1,3))*ZQDM(:,2:klev+1) 1177 d_qs(:,1:klev)=d_qs(:,1:klev) + (ZRXS(:,2:klev+1,5)-ZRXS0(:,2:klev+1,5))*ZQDM(:,2:klev+1) 1178 d_qg(:,1:klev)=d_qg(:,1:klev) + (ZRXS(:,2:klev+1,6)-ZRXS0(:,2:klev+1,6))*ZQDM(:,2:klev+1) 1179 ! Tendency, theta -> T 1180 d_t(:,1:klev)=d_t(:,1:klev) + (zthetas(:,2:klev+1)-zthetas0(:,2:klev+1))*ZEXN(:,2:klev+1) 1181 ! TKE 1182 d_tke(:,1:klev)=d_tke(:,1:klev) + (ZTKES(:,2:klev+1) - ZTKES0(:,2:klev+1)) 1183 1184 !Time evolution 1185 ZQR(:,:)=ZQR(:,:)+d_qr(:,:)*pdtphys 1186 ZQS(:,:)=ZQS(:,:)+d_qs(:,:)*pdtphys 1187 ZQG(:,:)=ZQG(:,:)+d_qg(:,:)*pdtphys 1188 PTKEM(:,2:klev+1)=PTKEM(:,2:klev+1)+d_tke(:,:)*pdtphys 1189 ! 125 1190 !------------------------------------------------------------ 126 1191 ! Entrees sorties 127 1192 !------------------------------------------------------------ 128 1193 129 130 call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t) 131 1194 itap=itap+1 1195 !zjulian=zjulian+pdtphys/86400. 1196 !print*,'avant mod(itap,iecri) , zjulian itap',zjulian,itap,iecri 1197 if ( mod(itap,iecri) == 1 .or. iecri == 1 ) then 1198 print*,'ECRITURE ITAP=',itap 1199 call output_physiqex(debut,zjulian,pdtphys_,presnivs,paprs,u,v,t,qx,ZCLDFR,ZQR,ZQS,ZQG,PTKEM,ZTHETA) 1200 call iophys_ecrit('tsrf', 1,'tsrf', ' ',tsrf) 1201 call iophys_ecrit('capcal', 1,'capcal', ' ',capcal) 1202 call iophys_ecrit('swdnsrf', 1,'swdnsrf', ' ',swdnsrf) 1203 call iophys_ecrit('lwdnsrf', 1,'lwdnsrf', ' ',lwdnsrf) 1204 call iophys_ecrit('lwupsrf', 1,'lwupsrf', ' ',lwupsrf) 1205 call iophys_ecrit('sensible', 1,'sensible', ' ',sensible) 1206 call iophys_ecrit('latent', 1,'latent', ' ',latent) 1207 call iophys_ecrit('qsats', 1,'qsats', ' ',qsats) 1208 call iophys_ecrit('tsoil', nlev,'tsoi', ' ',tsoil_out) 1209 !AI 1210 if (iflag_radia.ge.1) then 1211 call iophys_ecrit('lwdn0', 1,'lwdn0', ' ',lwdn0) 1212 call iophys_ecrit('lwdn', 1,'lwdn', ' ',lwdn) 1213 call iophys_ecrit('lwup0', 1,'lwup0', ' ',lwup0) 1214 call iophys_ecrit('lwup', 1,'lwup', ' ',lwup) 1215 call iophys_ecrit('swdn0', 1,'swdn0', ' ',swdn0) 1216 call iophys_ecrit('swdn', 1,'swdn', ' ',swdn) 1217 call iophys_ecrit('swup0', 1,'swup0', ' ',swup0) 1218 call iophys_ecrit('swup', 1,'swup', ' ',swup) 1219 endif 1220 endif 132 1221 133 1222 ! if lastcall, then it is time to write "restartphy.nc" file … … 138 1227 139 1228 end subroutine physiqex 1229 ! 1230 SUBROUTINE VERTICAL_EXTEND(PX,KLEV) 1231 1232 ! fill extra vetical levels to fit MNH interface 1233 1234 REAL, DIMENSION(:,:), INTENT(INOUT) :: PX 1235 INTEGER, INTENT(IN) :: KLEV 1236 PX(:,1 )= PX(:,2) 1237 PX(:,KLEV+2)= PX(:,KLEV+1) 1238 END SUBROUTINE VERTICAL_EXTEND 140 1239 141 1240 END MODULE physiqex_mod
Note: See TracChangeset
for help on using the changeset viewer.
