- Timestamp:
- Sep 18, 2024, 4:08:09 PM (7 weeks 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
r5197 r5198 1 ! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $ 1 2 MODULE physiqex_mod 3 2 4 IMPLICIT NONE 5 3 6 CONTAINS 4 7 5 ! ------------------------------------------------------------------------------------------------6 !7 ! Main authors : F. Hourdin, S. Riette, Q. Libois, A. Idelkadi8 ! --------------9 !10 ! Object : A version of the LMDZ physics that include PHYEX11 ! ------- Tu be used in particular to be coupled to non hydrostatic version of the dynamical12 ! 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 only19 ! 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 coming21 ! 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 cumulus25 ! Lège 2024 : first versions with a generic soil model (developped at this occasion) and ECrad26 !27 !28 ! To do list29 ! ---------30 ! To do list pour un branchement plus propre31 ! * PHYEX considère toutes les espèces microphysiques et la tke comme pronostiques, des termes de tendances32 ! sont calculés. L'avance temporelle est faite en fin de pas de temps mais pourrait être déplacée au dessus si33 ! les tendances de ces variables passaient par l'interface34 ! * PHYEX a besoin de variables avec mémoire d'un pas de temps à l'autre (en plus des variables pronostiques35 ! décrites juste au dessus). Ce sont les tableaux ALLOCATABLE. Ces variables pourraient devenir36 ! des traceurs.37 ! * La variable ZDZMIN est ici calculée avec un MINVAL qui conduira à des résultats différents38 ! lorsque la distribution (noeuds/procs) sera différente. Cette variable est utile pour déterminer39 ! un critère CFL pour la sédimentation des précipitations. Il suffit donc d'avoir une valeur40 ! approchante.41 ! * Certains allocatable sont en klev, d'autres en klev+2. Il est possible de changer ceci mais il faut gérer42 ! les recopies pour que les params voient des tableaux en klev+2 (en effectuant des recopies des43 ! 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 dans45 ! les variables ZINPRC, ZINPRR, ZINPRS, ZINPRG46 !47 !-------------------------------------------------------------------------------------------------48 49 50 8 SUBROUTINE physiqex (nlon,nlev, & 51 & debut,lafin,pdtphys _, &9 & debut,lafin,pdtphys, & 52 10 & paprs,pplay,pphi,pphis,presnivs, & 53 11 & u,v,rot,t,qx, & … … 57 15 USE dimphy, only : klon,klev 58 16 USE infotrac_phy, only : nqtot 59 USE geometry_mod, only : latitude , cell_area, latitude_deg, longitude_deg17 USE geometry_mod, only : latitude 60 18 ! USE comcstphy, only : rg 61 19 USE ioipsl, only : ymds2ju 62 !USE phys_state_var_mod, only : phys_state_var_init 63 USE phys_state_var_mod 20 USE phys_state_var_mod, only : phys_state_var_init 64 21 USE phyetat0_mod, only: phyetat0 65 22 USE output_physiqex_mod, ONLY: output_physiqex 66 !!!!!! AI 04-202467 USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time68 USE limit_read_mod, ONLY : init_limit_read69 USE conf_phys_m, only: conf_phys70 USE ozonecm_m, only: ozonecm71 USE aero_mod72 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, calend74 USE radlwsw_m, only: radlwsw75 USE MODD_CST76 USE phys_output_var_mod, ONLY : phys_output_var_init, cloud_cover_sw77 USE write_field_phy78 USE phyaqua_mod, only: zenang_an79 !!!!!!80 ! PHYEX internal modules81 USE MODE_INIT_PHYEX, ONLY: INIT_PHYEX, FILL_DIMPHYEX82 USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t83 USE MODD_PHYEX, ONLY: PHYEX_t84 USE MODI_INI_PHYEX, ONLY: INI_PHYEX85 USE MODI_ICE_ADJUST, ONLY: ICE_ADJUST86 USE MODI_RAIN_ICE, ONLY: RAIN_ICE87 USE MODI_TURB88 USE MODI_SHALLOW_MF89 90 USE MODD_CST91 USE ioipsl_getin_p_mod, ONLY : getin_p92 USE generic_soil_ini, ONLY : soil_ini93 USE generic_soil_conduct, ONLY : soil_conduct94 USE generic_soil_fluxes, ONLY : soil_fluxes95 23 96 24 IMPLICIT none … … 98 26 ! Routine argument: 99 27 ! 28 100 29 integer,intent(in) :: nlon ! number of atmospheric colums 101 30 integer,intent(in) :: nlev ! number of vertical levels (should be =klev) 102 31 logical,intent(in) :: debut ! signals first call to physics 103 32 logical,intent(in) :: lafin ! signals last call to physics 104 real,intent(in) :: pdtphys _! physics time step (s)33 real,intent(in) :: pdtphys ! physics time step (s) 105 34 real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) 106 35 real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) … … 119 48 real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers 120 49 real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure 121 real :: d_qr(klon, klev), d_qs(klon, klev), d_qg(klon, klev) ! tendency for rain, snow, graupel122 real :: d_tke(klon, klev)123 50 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 51 ! include "clesphys.h" 274 52 INTEGER length 275 53 PARAMETER ( length = 100 ) … … 279 57 !$OMP THREADPRIVATE(clesphy0) 280 58 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 ! 59 385 60 real :: temp_newton(klon,klev) 386 61 integer :: k 387 62 logical, save :: first=.true. 388 63 !$OMP THREADPRIVATE(first) 389 !real,save :: rg=9.81 390 !!$OMP THREADPRIVATE(rg) 64 65 real,save :: rg=9.81 66 !$OMP THREADPRIVATE(rg) 391 67 392 68 ! For I/Os 393 69 integer :: itau0 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) 70 real :: zjulian 420 71 421 72 … … 425 76 426 77 print*,'Debut physiqex',debut 427 ! AI428 pdtphys=pdtphys_429 CALL update_time(pdtphys)430 phys_tstep=NINT(pdtphys)431 432 78 ! initializations 433 79 if (debut) then ! Things to do only for the first call to physics 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) 80 print*,'Debut physiqex IN' 449 81 450 ! IF (.NOT. create_etat0_limit) 451 CALL init_limit_read(days_elapsed) 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) 452 85 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 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) 478 91 479 ! Initialize outputs:480 itau0=0481 ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0482 !CALL ymds2ju(annee0, month, dayref, hour, zjulian)483 !call ymds2ju(1979, 1, 1, 0.0, zjulian)484 !print*,'zjulian=',zjulian485 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 saved492 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 state508 509 ! ----------------------------------------------------------------510 ! Initialisation du sol generique511 ! ----------------------------------------------------------------512 allocate(tsrf(nlon))513 allocate(tsoil(nlon,nsoil))514 thermal_inertia0=2000.515 z0m0=0.05516 z0h0=0.05517 z0q0=0.05518 betasoil=0.2519 call soil_ini(klon,nsoil,PHYEX%CST%XCPD,PHYEX%CST%XLVTT,thermal_inertia0,z0m0,z0h0,z0q0,betasoil)520 !521 92 #ifndef CPP_IOIPSL_NO_OUTPUT 522 93 ! Initialize IOIPSL output file 523 94 #endif 524 ! AI Partie Surface depplacee525 ! compute tendencies to return to the dynamics:526 ! "friction" on the first layer527 ! tsrf(1:nlon)=t(1:nlon,1)528 ! do isoil=1,nsoil529 ! tsoil(1:nlon,isoil)=tsrf(1:nlon)530 ! !tsoil(1:nlon,isoil)=302.531 ! enddo532 ! rpi=2.*asin(1.)533 ! timesec=5*3600534 ! stefan=5.67e-8535 ! fsw0=900.536 ! call getin_p('iecri',iecri)537 ! call getin_p('iflag_surface',iflag_surface)538 95 539 96 endif ! of if (debut) 540 97 541 ! AI542 ! Increment543 !itap = itap + 1544 98 !------------------------------------------------------------ 545 99 ! Initialisations a chaque pas de temps 546 100 !------------------------------------------------------------ 101 547 102 548 103 ! set all tendencies to zero … … 551 106 d_t(1:klon,1:klev)=0. 552 107 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.556 108 d_ps(1:klon)=0. 557 d_tke(1:klon,1:klev)=0.558 !559 ! AI attention560 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 clouds578 ZSVT(:,:,:) = 0.579 PWT(:,:) = 0.580 ZUT(:,2:klev+1) = u(:,:)581 ZVT(:,2:klev+1) = v(:,:)582 !583 !------------------------------------------------------------584 ! Conversions and extra levels585 !------------------------------------------------------------586 !TODO check in Meso-NH how values are extrapolated outside of the physical domain587 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,KRR598 CALL VERTICAL_EXTEND(ZRX(:,:,JRR),klev)599 END DO600 !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 domain607 zs(:) = pphis(:)/PHYEX%CST%XG608 zz_mass(:,2:klev+1) = pphi(:,:) / PHYEX%CST%XG609 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+2613 zz_flux(:,k)=(zz_mass(:,k-1)+zz_mass(:,k))/2.614 enddo615 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+1619 zdzf(:,k)=zz_flux(:,k+1)-zz_flux(:,k)620 enddo621 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+2625 zdzm(:,k)=zz_mass(:,k)-zz_mass(:,k-1)626 enddo627 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+1632 PRHODJ(:,k) = ZRHOD(:,k) * (zdzf(:,k)*cell_area(:))633 END DO634 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 109 643 110 !------------------------------------------------------------ 644 ! Tendencies111 ! Calculs 645 112 !------------------------------------------------------------ 646 !For Meso-NH, initialia values for the tendencies are filled with647 !a pseudo-tendecy computed by dividing the state variable by the time step648 !This mechanism enables the possibility for the different parametrisations649 !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 needed651 !as the parametrisations don't use the S varaible to guess the futur value of the wind.652 ZRXS(:,:,:) = ZRX(:,:,:)/pdtphys653 ZTHETAS(:,:)=ZTHETA(:,:)/pdtphys654 ZTKES(:,:)=PTKEM(:,:)/pdtphys655 !To compute the actual tendency, we save the initial values of these variables656 ZRXS0(:,:,:) = ZRXS(:,:,:)657 ZTHETAS0(:,:)=ZTHETAS(:,:)658 ZTKES0(:,:)=ZTKES(:,:)659 !------------------------------------------------------------660 ! Adjustment661 !------------------------------------------------------------662 !663 ZSRC(:,:) = 0.664 ZSIGQSAT=PHYEX%NEBN%VSIGQSAT665 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(:,:)*pdtphys682 ZRX(:,:,:)=ZRXS(:,:,:)*pdtphys683 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 ! Radiation695 !------------------------------------------------------------696 !697 !AI a nettoyer apres698 ! calcul de qsat Fred mai 2024699 zqsat(:,:) = EXP( XALPW - XBETAW / t(:,:) - XGAMW * ALOG( t(:,:) ) )700 zqsat(:,:) = XRD / XRV * zqsat(:,:) / ( pplay(:,:) - zqsat(:,:))701 !702 if (iflag_radia.ge.1) then703 !704 !ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors705 !Cloud fraction is available in the ZCLDFR variable706 !707 ! AI attention : dans un 1er temps read_climoz=0708 ! Update ozone if day change709 !print*,'solarlong0, itap, lmt_pas, days_elapsed : ', solarlong0, itap, lmt_pas, days_elapsed710 !stop711 IF (MOD(itap-1,lmt_pas) == 0) THEN712 ! IF (read_climoz <= 0) THEN713 ! Once per day, update ozone from Royer:714 IF (solarlong0<-999.) then715 ! Generic case with evolvoing season716 zzzo=real(days_elapsed+1)717 ELSE IF (abs(solarlong0-1000.)<1.e-4) then718 ! Particular case with annual mean insolation719 zzzo=real(90) ! could be revisited720 IF (read_climoz/=-1) THEN721 abort_message ='read_climoz=-1 is recommended when ' &722 // 'solarlong0=1000.'723 CALL abort_physic (modname,abort_message,1)724 ENDIF725 ELSE726 ! Case where the season is imposed with solarlong0727 zzzo=real(90) ! could be revisited728 ENDIF729 zzzo=real(days_elapsed+1)730 wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzzo)731 !print*,'wo = ', wo732 ! ELSE733 ! !--- ro3i = elapsed days number since current year 1st january, 0h734 ! ro3i=days_elapsed+jh_cur-jh_1jan735 ! !--- scaling for old style files (360 records)736 ! IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len737 ! IF(adjust_tropopause) THEN738 ! 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 ! ELSE743 ! 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 ! ENDIF747 ! ! Convert from mole fraction of ozone to column density of ozone in a748 ! ! cell, in kDU:749 ! FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &750 ! * zmasse / dobson_u / 1e3751 ! ! (By regridding ozone values for LMDZ only once a day, we752 ! ! have already neglected the variation of pressure in one753 ! ! day. So do not recompute "wo" at each time step even if754 ! ! "zmasse" changes a little.)755 ! ENDIF756 ENDIF757 !=========================================================================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_surface761 !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_cur764 day_since_equinox = (jD_cur + jH_cur) - jD_eq765 !choix entre calcul de la longitude solaire vraie ou valeur fixee = solarlong0766 IF (solarlong0<-999.) THEN767 IF (new_orbit) THEN768 ! calcul selon la routine utilisee pour les planetes769 CALL solarlong(day_since_equinox, zlongi, dist)770 print*,'day_since_equinox, zlongi, dist : ', day_since_equinox, zlongi, dist771 ELSE772 ! calcul selon la routine utilisee pour l'AR4773 CALL orbite(REAL(days_elapsed+1),zlongi,dist)774 ENDIF775 ELSE776 zlongi=solarlong0 ! longitude solaire vraie777 dist=1. ! distance au soleil / moyenne778 ENDIF779 !IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist780 113 781 ! AI a refaire en definissant les itap, itaprad782 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!783 ! ============================784 ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur785 ! l'annee a partir d'une formule analytique.786 ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et787 ! non nul aux poles.788 IF (abs(solarlong0-1000.)<1.e-4) THEN789 CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &790 latitude_deg,longitude_deg,rmu0,fract)791 swradcorr(:) = 1.0792 ! JrNt(:) = 1.0793 zrmu0(:) = rmu0(:)794 ELSE795 ! recode par Olivier Boucher en sept 2015796 SELECT CASE (iflag_cycle_diurne)797 CASE(0)798 ! Sans cycle diurne799 CALL angle(zlongi, latitude_deg, fract, rmu0)800 swradcorr = 1.0801 ! JrNt = 1.0802 zrmu0 = rmu0803 CASE(1)804 ! Avec cycle diurne sans application des poids805 ! bit comparable a l ancienne formulation cycle_diurne=true806 ! on integre entre gmtime et gmtime+radpas807 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 = rmu0811 swradcorr = 1.0812 ! Calcul du flag jour-nuit813 ! JrNt = 0.0814 !WHERE (fract.GT.0.0) JrNt = 1.0815 CASE(2)816 ! Avec cycle diurne sans application des poids817 ! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)818 ! Comme cette routine est appele a tous les pas de temps de819 ! la physique meme si le rayonnement n'est pas appele je820 ! remonte en arriere les radpas-1 pas de temps821 ! suivant. Petite ruse avec MOD pour prendre en compte le822 ! premier pas de temps de la physique pendant lequel823 ! itaprad=0824 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 poids830 !831 zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le832 zdtime2=0.0 !--pas de temps de la physique qui se termine833 CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &834 latitude_deg,longitude_deg,zrmu0,zfract)835 swradcorr = 0.0836 WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &837 swradcorr=zfract/fract*zrmu0/rmu0838 ! Calcul du flag jour-nuit839 !JrNt = 0.0840 !WHERE (zfract.GT.0.0) JrNt = 1.0841 END SELECT842 ENDIF843 !print*,'fract = ',fract844 !sza_o = ACOS (rmu0) *180./pi845 846 ! Calcul de l'ensoleillement :847 IF (MOD(itaprad,radpas).EQ.0) THEN848 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 surface854 ! zxtsol temp au sol : dvrait etre recuperee de la partie surface855 ! albsol_dir, albsol_dif aussi856 albsol_dir = 1. ! AI attention857 albsol_dif = 1. ! AI attention858 zxtsol(:) = t(:,1) ! AI attention859 !860 ! Q2 ozone ? lmdz ok a completer861 ! wo : pour l'instant appel ozonecm862 !863 ! Q3 orbite ? lmdz ok a comleter864 ! dist, rmu0 et fract865 ! pour l'instant cas solarlong0<-999.866 ! comment introduire les pas de temps rayonnement867 !868 ! Q3 proprietes optiques des nuages ? call cloud_optics_prop.F90 ? non necessaires pour ecrad869 ! epaisseur optique, emissivite : cldfrarad, cldemirad, cldtaurad , non necessaire pour ecrad870 ! 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'IFS873 ! 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 : zqsat891 ! zqsat specific humidity at liq saturation =>892 ! calcul de Fred ou a calculer dans ecrad enf fct de t et p893 ! flwc, fiwc : qx( 2) et qx( 3)894 895 print*,'Entree radlwsw a itap : ', itap896 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 RRTM906 & 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 diagnostics927 & toplwad_aero, sollwad_aero,&928 & toplwai_aero, sollwai_aero, &929 & toplwad0_aero, sollwad0_aero,&930 !-end931 & ZLWFT0_i, ZFLDN0, ZFLUP0, &932 & ZSWFT0_i, ZFSDN0, ZFSUP0, &933 & cloud_cover_sw)934 itaprad = 0935 ENDIF ! IF (MOD(itaprad,radpas).EQ.0)936 itaprad = itaprad + 1937 938 IF (iflag_radia.eq.0) THEN939 heat=0.940 cool=0.941 sollw=0. ! MPL 01032011942 solsw=0.943 radsol=0.944 swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars945 swup0=0.946 lwup=0.947 lwup0=0.948 lwdn=0.949 lwdn0=0.950 ENDIF951 ! Ajouter la tendance des rayonnements (tous les pas)952 ! avec une correction pour le cycle diurne dans le SW953 !954 DO k=1, klev955 ! AI attention956 !d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY957 !d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY958 !d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY959 !d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY960 d_t_swr(:,k)=swradcorr(:)*heat(:,k)/RDAY961 d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)/RDAY962 d_t_lwr(:,k)=-cool(:,k)/RDAY963 d_t_lw0(:,k)=-cool0(:,k)/RDAY964 ENDDO965 966 !t_seri(:,:) = t_seri(:,:) + d_t_swr(:,:) + d_t_lwr(:,:)967 d_t = d_t + d_t_swr + d_t_lwr968 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 nettoyer974 !975 ! AI passer peut-etre la surface avant radiation976 ! albedo, temp surf, ... pour le rayonnement977 !978 !ECRAD can be plugged here and can use the adjusted values for temperature and hydrometeors979 !Cloud fraction is available in the ZCLDFR variable980 981 !982 !------------------------------------------------------------983 ! Surface984 !------------------------------------------------------------985 !986 ! AI987 114 ! compute tendencies to return to the dynamics: 988 115 ! "friction" on the first layer 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 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 1002 123 1003 timesec=timesec+pdtphys1004 124 1005 ! Computation of qsats. Done even it is not used for ouptuts1006 ! AI1007 !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 ) then1012 ! 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 temperature1015 ! A case close to Rico can be obtained by imposing fsens=-5E-3 et flat=-6E-51016 capcal(:)=1.e61017 fluxgrd(:)=0.1018 sensible(:)=-fsens1019 latent(:)=-flat1020 swdnsrf(:)=0.1021 lwdnsrf(:)=0.1022 fluxt(1:klon)=sensible/PHYEX%CST%XCPD1023 fluxq(1:klon)=latent(:)/PHYEX%CST%XLVTT1024 fluxu(1:klon)=-0.01*u(1:klon,1)1025 fluxv(1:klon)=-0.01*v(1:klon,1)1026 else1027 ! AI1028 if (iflag_radia.eq.0) then1029 ! Using the generic soil model1030 if ( nlon == 1 ) then1031 ! Imposed radiative fluxes to simulate a 1D case close to Arm Cumulus1032 swdnsrf(:)=fsw0*max(cos(2*rpi*(timesec-43200.)/(1.15*86400.)),0.)1033 lwdnsrf(:)=400.1034 else1035 ! Model forced in latitude by an idelaized radiative flux that depends on latitude1036 ! only. Used for 3D tests without radiation.1037 swdnsrf(:)=600*cos(latitude(:))+1001038 lwdnsrf(:)=0.1039 endif1040 else1041 swdnsrf(:)=solsw(:)1042 lwdnsrf(:)=sollwdown(:)1043 endif1044 ! updating the soil temperature (below the surface) from the surface temperature1045 ! en returning the surface conductive flux and a heat capacity (accounting for the1046 ! 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 outputs1050 1051 ! Computing the turbulent fluxes Fqx = rho w'x' = rho (kappa/ln(z/z0x))^2 [ x_surface - x ]. Upward1052 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%XCPD1054 latent(:)=fluxq(:)*PHYEX%CST%XLVTT1055 endif1056 tsoil_out=0. ; tsoil_out(:,1:min(nsoil,nlev))=tsoil(:,1:min(nsoil,nlev))1057 1058 ! Computing the surface temperature1059 lwupsrf=stefan*tsrf(:)**41060 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 convection1070 !------------------------------------------------------------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 tendency1094 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 tendencies1099 !1100 !------------------------------------------------------------1101 ! Turbulence1102 !------------------------------------------------------------1103 ! out tendencies1104 ZRUS(:,:) = 0.1105 ZRVS(:,:) = 0.1106 ZRWS(:,:) = 0.1107 ZRSVS(:,:,:) = 0.1108 ZRTKES(:,:) = ZTKES(:,:) * PRHODJ(:,:)1109 DO JRR=1, KRR1110 ZRRS(:,:,JRR) = ZRXS(:,:,JRR) * PRHODJ(:,:)1111 ENDDO1112 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, KRR1138 ZRXS(:,:,JRR) = ZRRS(:,:,JRR) / PRHODJ(:,:)1139 ENDDO1140 ZTHETAS(:,:) = ZRTHS(:,:) / PRHODJ(:,:)1141 ZTKES(:,:) = ZRTKES(:,:) / PRHODJ(:,:)1142 ! Add tendencies of turb to total physics tendency1143 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') THEN1146 PSIGS(:,:)=SQRT(PSIGS(:,:)**2 + PSIGMF(:,:)**2)1147 ENDIF1148 !------------------------------------------------------------1149 ! Microphysics1150 !------------------------------------------------------------1151 ZSEA=1.1152 ZTOWN=0.1153 ZCIT=0.1154 ZTHVREFZIKB=01155 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 -> specific1173 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 -> T1180 d_t(:,1:klev)=d_t(:,1:klev) + (zthetas(:,2:klev+1)-zthetas0(:,2:klev+1))*ZEXN(:,2:klev+1)1181 ! TKE1182 d_tke(:,1:klev)=d_tke(:,1:klev) + (ZTKES(:,2:klev+1) - ZTKES0(:,2:klev+1))1183 1184 !Time evolution1185 ZQR(:,:)=ZQR(:,:)+d_qr(:,:)*pdtphys1186 ZQS(:,:)=ZQS(:,:)+d_qs(:,:)*pdtphys1187 ZQG(:,:)=ZQG(:,:)+d_qg(:,:)*pdtphys1188 PTKEM(:,2:klev+1)=PTKEM(:,2:klev+1)+d_tke(:,:)*pdtphys1189 !1190 125 !------------------------------------------------------------ 1191 126 ! Entrees sorties 1192 127 !------------------------------------------------------------ 1193 128 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 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 1221 132 1222 133 ! if lastcall, then it is time to write "restartphy.nc" file … … 1227 138 1228 139 end subroutine physiqex 1229 !1230 SUBROUTINE VERTICAL_EXTEND(PX,KLEV)1231 1232 ! fill extra vetical levels to fit MNH interface1233 1234 REAL, DIMENSION(:,:), INTENT(INOUT) :: PX1235 INTEGER, INTENT(IN) :: KLEV1236 PX(:,1 )= PX(:,2)1237 PX(:,KLEV+2)= PX(:,KLEV+1)1238 END SUBROUTINE VERTICAL_EXTEND1239 140 1240 141 END MODULE physiqex_mod
Note: See TracChangeset
for help on using the changeset viewer.
