| 1 | ! |
|---|
| 2 | ! $Id: $ |
|---|
| 3 | ! |
|---|
| 4 | MODULE physiq_mod |
|---|
| 5 | |
|---|
| 6 | IMPLICIT NONE |
|---|
| 7 | |
|---|
| 8 | CONTAINS |
|---|
| 9 | |
|---|
| 10 | SUBROUTINE physiq (nlon,nlev,nqmax, |
|---|
| 11 | . debut,lafin,rjourvrai,gmtime,pdtphys, |
|---|
| 12 | . paprs,pplay,ppk,pphi,pphis,presnivs, |
|---|
| 13 | . u,v,t,qx, |
|---|
| 14 | . flxmw, |
|---|
| 15 | . d_u, d_v, d_t, d_qx, d_ps) |
|---|
| 16 | |
|---|
| 17 | c====================================================================== |
|---|
| 18 | c |
|---|
| 19 | c Modifications pour la physique de Venus |
|---|
| 20 | c S. Lebonnois (LMD/CNRS) Septembre 2005 |
|---|
| 21 | c |
|---|
| 22 | c --------------------------------------------------------------------- |
|---|
| 23 | c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 |
|---|
| 24 | c |
|---|
| 25 | c Objet: Moniteur general de la physique du modele |
|---|
| 26 | cAA Modifications quant aux traceurs : |
|---|
| 27 | cAA - uniformisation des parametrisations ds phytrac |
|---|
| 28 | cAA - stockage des moyennes des champs necessaires |
|---|
| 29 | cAA en mode traceur off-line |
|---|
| 30 | c modif ( P. Le Van , 12/10/98 ) |
|---|
| 31 | c |
|---|
| 32 | c Arguments: |
|---|
| 33 | c |
|---|
| 34 | c nlon----input-I-nombre de points horizontaux |
|---|
| 35 | c nlev----input-I-nombre de couches verticales |
|---|
| 36 | c nqmax---input-I-nombre de traceurs |
|---|
| 37 | c debut---input-L-variable logique indiquant le premier passage |
|---|
| 38 | c lafin---input-L-variable logique indiquant le dernier passage |
|---|
| 39 | c rjour---input-R-numero du jour de l'experience |
|---|
| 40 | c gmtime--input-R-fraction de la journee (0 a 1) |
|---|
| 41 | c pdtphys-input-R-pas d'integration pour la physique (seconde) |
|---|
| 42 | c paprs---input-R-pression pour chaque inter-couche (en Pa) |
|---|
| 43 | c pplay---input-R-pression pour le mileu de chaque couche (en Pa) |
|---|
| 44 | c ppk ---input-R-fonction d'Exner au milieu de couche |
|---|
| 45 | c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol) |
|---|
| 46 | c pphis---input-R-geopotentiel du sol |
|---|
| 47 | c presnivs-input_R_pressions approximat. des milieux couches ( en PA) |
|---|
| 48 | c u-------input-R-vitesse dans la direction X (de O a E) en m/s |
|---|
| 49 | c v-------input-R-vitesse Y (de S a N) en m/s |
|---|
| 50 | c t-------input-R-temperature (K) |
|---|
| 51 | c qx------input-R-mass mixing ratio traceurs (kg/kg) |
|---|
| 52 | c d_t_dyn-input-R-tendance dynamique pour "t" (K/s) |
|---|
| 53 | c flxmw---input-R-flux de masse vertical en kg/s |
|---|
| 54 | c |
|---|
| 55 | c d_u-----output-R-tendance physique de "u" (m/s/s) |
|---|
| 56 | c d_v-----output-R-tendance physique de "v" (m/s/s) |
|---|
| 57 | c d_t-----output-R-tendance physique de "t" (K/s) |
|---|
| 58 | c d_qx----output-R-tendance physique de "qx" (kg/kg/s) |
|---|
| 59 | c d_ps----output-R-tendance physique de la pression au sol |
|---|
| 60 | c====================================================================== |
|---|
| 61 | USE ioipsl |
|---|
| 62 | ! USE histcom ! not needed; histcom is included in ioipsl |
|---|
| 63 | use dimphy |
|---|
| 64 | USE geometry_mod,only: longitude, latitude, ! in radians |
|---|
| 65 | & longitude_deg,latitude_deg, ! in degrees |
|---|
| 66 | & cell_area,dx,dy |
|---|
| 67 | USE phys_state_var_mod ! Variables sauvegardees de la physique |
|---|
| 68 | USE cpdet_phy_mod, only: cpdet, t2tpot |
|---|
| 69 | USE chemparam_mod |
|---|
| 70 | USE age_of_air_mod, ONLY: ok_aoa, reinit_aoa, i_aoa, init_age |
|---|
| 71 | USE age_of_air_mod, ONLY: aoa_ini, age_of_air |
|---|
| 72 | USE conc |
|---|
| 73 | USE param_v4_h |
|---|
| 74 | USE compo_hedin83_mod2 |
|---|
| 75 | use radlwsw_newtoncool_mod, only: radlwsw_newtoncool |
|---|
| 76 | ! use ieee_arithmetic |
|---|
| 77 | use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy |
|---|
| 78 | use mod_grid_phy_lmdz, only: nbp_lon |
|---|
| 79 | use infotrac_phy, only: iflag_trac, tname, ttext |
|---|
| 80 | use vertical_layers_mod, only: pseudoalt |
|---|
| 81 | use turb_mod, only : sens, turb_resolved |
|---|
| 82 | use nonoro_gwd_ran_mod, only: nonoro_gwd_ran |
|---|
| 83 | use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation |
|---|
| 84 | use iono_h, only: temp_elect, temp_ion |
|---|
| 85 | #ifdef CPP_XIOS |
|---|
| 86 | use xios_output_mod, only: initialize_xios_output, |
|---|
| 87 | & update_xios_timestep, |
|---|
| 88 | & send_xios_field |
|---|
| 89 | use wxios, only: wxios_context_init, xios_context_finalize |
|---|
| 90 | #endif |
|---|
| 91 | #ifdef MESOSCALE |
|---|
| 92 | use comm_wrf |
|---|
| 93 | #else |
|---|
| 94 | use iophy |
|---|
| 95 | use write_field_phy |
|---|
| 96 | use mod_phys_lmdz_omp_data, ONLY: is_omp_master |
|---|
| 97 | USE mod_phys_lmdz_para, only : is_parallel,jj_nb, |
|---|
| 98 | & is_north_pole_phy, |
|---|
| 99 | & is_south_pole_phy, |
|---|
| 100 | & is_master |
|---|
| 101 | #endif |
|---|
| 102 | IMPLICIT none |
|---|
| 103 | c====================================================================== |
|---|
| 104 | c CLEFS CPP POUR LES IO |
|---|
| 105 | c ===================== |
|---|
| 106 | #ifndef MESOSCALE |
|---|
| 107 | c#define histhf |
|---|
| 108 | #define histday |
|---|
| 109 | #define histmth |
|---|
| 110 | #define histins |
|---|
| 111 | #endif |
|---|
| 112 | c====================================================================== |
|---|
| 113 | #include "dimsoil.h" |
|---|
| 114 | #include "clesphys.h" |
|---|
| 115 | #include "iniprint.h" |
|---|
| 116 | #include "timerad.h" |
|---|
| 117 | #include "tabcontrol.h" |
|---|
| 118 | #include "nirdata.h" |
|---|
| 119 | #include "hedin.h" |
|---|
| 120 | c====================================================================== |
|---|
| 121 | LOGICAL ok_journe ! sortir le fichier journalier |
|---|
| 122 | save ok_journe |
|---|
| 123 | c PARAMETER (ok_journe=.true.) |
|---|
| 124 | c |
|---|
| 125 | LOGICAL ok_mensuel ! sortir le fichier mensuel |
|---|
| 126 | save ok_mensuel |
|---|
| 127 | c PARAMETER (ok_mensuel=.true.) |
|---|
| 128 | c |
|---|
| 129 | LOGICAL ok_instan ! sortir le fichier instantane |
|---|
| 130 | save ok_instan |
|---|
| 131 | c PARAMETER (ok_instan=.true.) |
|---|
| 132 | c |
|---|
| 133 | c====================================================================== |
|---|
| 134 | c |
|---|
| 135 | c Variables argument: |
|---|
| 136 | c |
|---|
| 137 | INTEGER,INTENT(IN) :: nlon |
|---|
| 138 | INTEGER,INTENT(IN) :: nlev |
|---|
| 139 | INTEGER,INTENT(IN) :: nqmax |
|---|
| 140 | REAL,INTENT(IN) :: rjourvrai |
|---|
| 141 | REAL,INTENT(IN) :: gmtime |
|---|
| 142 | REAL,INTENT(IN) :: pdtphys |
|---|
| 143 | LOGICAL,INTENT(IN) :: debut, lafin |
|---|
| 144 | REAL,INTENT(IN) :: paprs(klon,klev+1) |
|---|
| 145 | REAL,INTENT(IN) :: pplay(klon,klev) |
|---|
| 146 | REAL,INTENT(IN) :: pphi(klon,klev) |
|---|
| 147 | REAL,INTENT(IN) :: pphis(klon) |
|---|
| 148 | REAL,INTENT(IN) :: presnivs(klev) |
|---|
| 149 | |
|---|
| 150 | ! ADAPTATION GCM POUR CP(T) |
|---|
| 151 | REAL ppk(klon,klev) |
|---|
| 152 | |
|---|
| 153 | REAL,INTENT(IN) :: u(klon,klev) |
|---|
| 154 | REAL,INTENT(IN) :: v(klon,klev) |
|---|
| 155 | REAL,INTENT(IN) :: t(klon,klev) |
|---|
| 156 | REAL,INTENT(IN) :: qx(klon,klev,nqmax) |
|---|
| 157 | |
|---|
| 158 | REAL d_u_dyn(klon,klev) |
|---|
| 159 | REAL d_t_dyn(klon,klev) |
|---|
| 160 | |
|---|
| 161 | REAL,INTENT(IN) :: flxmw(klon,klev) |
|---|
| 162 | |
|---|
| 163 | REAL,INTENT(OUT) :: d_u(klon,klev) |
|---|
| 164 | REAL,INTENT(OUT) :: d_v(klon,klev) |
|---|
| 165 | REAL,INTENT(OUT) :: d_t(klon,klev) |
|---|
| 166 | REAL,INTENT(OUT) :: d_qx(klon,klev,nqmax) |
|---|
| 167 | REAL,INTENT(OUT) :: d_ps(klon) |
|---|
| 168 | |
|---|
| 169 | logical ok_hf |
|---|
| 170 | real ecrit_hf |
|---|
| 171 | integer nid_hf |
|---|
| 172 | save ok_hf, ecrit_hf, nid_hf |
|---|
| 173 | |
|---|
| 174 | #ifdef histhf |
|---|
| 175 | data ok_hf,ecrit_hf/.true.,0.25/ |
|---|
| 176 | #else |
|---|
| 177 | data ok_hf/.false./ |
|---|
| 178 | #endif |
|---|
| 179 | |
|---|
| 180 | c Variables propres a la physique |
|---|
| 181 | |
|---|
| 182 | integer,save :: itap ! physics counter |
|---|
| 183 | REAL delp(klon,klev) ! epaisseur d'une couche |
|---|
| 184 | REAL omega(klon,klev) ! vitesse verticale en Pa/s (+ downward) |
|---|
| 185 | REAL vertwind(klon,klev) ! vitesse verticale en m/s (+ upward) |
|---|
| 186 | |
|---|
| 187 | INTEGER igwd,idx(klon),itest(klon) |
|---|
| 188 | |
|---|
| 189 | c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro |
|---|
| 190 | |
|---|
| 191 | REAL zulow(klon),zvlow(klon) |
|---|
| 192 | REAL zustrdr(klon), zvstrdr(klon) |
|---|
| 193 | REAL zustrli(klon), zvstrli(klon) |
|---|
| 194 | REAL zustrhi(klon), zvstrhi(klon) |
|---|
| 195 | REAL zublstrdr(klon), zvblstrdr(klon) |
|---|
| 196 | REAL znlow(klon), zeff(klon) |
|---|
| 197 | REAL zbl(klon), knu2(klon),kbreak(nlon) |
|---|
| 198 | REAL tau0(klon), ztau(klon,klev) |
|---|
| 199 | |
|---|
| 200 | c Pour calcul GW drag oro et nonoro: CALCUL de N2: |
|---|
| 201 | real zdtlev(klon,klev),zdzlev(klon,klev) |
|---|
| 202 | real ztlev(klon,klev) |
|---|
| 203 | real zn2(klon,klev) ! BV^2 at plev |
|---|
| 204 | |
|---|
| 205 | c Pour les bilans de moment angulaire, |
|---|
| 206 | integer bilansmc |
|---|
| 207 | c Pour le transport de ballons |
|---|
| 208 | integer ballons |
|---|
| 209 | c j'ai aussi besoin |
|---|
| 210 | c du stress de couche limite a la surface: |
|---|
| 211 | |
|---|
| 212 | REAL zustrcl(klon),zvstrcl(klon) |
|---|
| 213 | |
|---|
| 214 | c et du stress total c de la physique: |
|---|
| 215 | |
|---|
| 216 | REAL zustrph(klon),zvstrph(klon) |
|---|
| 217 | |
|---|
| 218 | c Variables locales: |
|---|
| 219 | c |
|---|
| 220 | REAL cdragh(klon) ! drag coefficient pour T and Q |
|---|
| 221 | REAL cdragm(klon) ! drag coefficient pour vent |
|---|
| 222 | c |
|---|
| 223 | cAA Pour TRACEURS |
|---|
| 224 | cAA |
|---|
| 225 | REAL,save,allocatable :: source(:,:) |
|---|
| 226 | REAL ycoefh(klon,klev) ! coef d'echange pour phytrac |
|---|
| 227 | REAL :: kzz_p(klev) ! coef d'echange pour phytrac pour le 1D |
|---|
| 228 | REAL yu1(klon) ! vents dans la premiere couche U |
|---|
| 229 | REAL yv1(klon) ! vents dans la premiere couche V |
|---|
| 230 | |
|---|
| 231 | REAL dsens(klon) ! derivee chaleur sensible |
|---|
| 232 | REAL ve(klon) ! integr. verticale du transport meri. de l'energie |
|---|
| 233 | REAL vq(klon) ! integr. verticale du transport meri. de l'eau |
|---|
| 234 | REAL ue(klon) ! integr. verticale du transport zonal de l'energie |
|---|
| 235 | REAL uq(klon) ! integr. verticale du transport zonal de l'eau |
|---|
| 236 | c |
|---|
| 237 | REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2) |
|---|
| 238 | |
|---|
| 239 | c====================================================================== |
|---|
| 240 | c |
|---|
| 241 | c Declaration des procedures appelees |
|---|
| 242 | c |
|---|
| 243 | EXTERNAL ajsec ! ajustement sec |
|---|
| 244 | EXTERNAL clmain ! couche limite |
|---|
| 245 | EXTERNAL clmain_ideal ! couche limite simple |
|---|
| 246 | EXTERNAL hgardfou ! verifier les temperatures |
|---|
| 247 | c EXTERNAL orbite ! calculer l'orbite |
|---|
| 248 | EXTERNAL phyetat0 ! lire l'etat initial de la physique |
|---|
| 249 | EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique |
|---|
| 250 | EXTERNAL radlwsw ! rayonnements solaire et infrarouge |
|---|
| 251 | ! EXTERNAL suphec ! initialiser certaines constantes |
|---|
| 252 | EXTERNAL transp ! transport total de l'eau et de l'energie |
|---|
| 253 | EXTERNAL printflag |
|---|
| 254 | EXTERNAL zenang |
|---|
| 255 | EXTERNAL diagetpq |
|---|
| 256 | EXTERNAL conf_phys |
|---|
| 257 | EXTERNAL diagphy |
|---|
| 258 | EXTERNAL mucorr |
|---|
| 259 | EXTERNAL nirco2abs |
|---|
| 260 | EXTERNAL nir_leedat |
|---|
| 261 | EXTERNAL nltecool |
|---|
| 262 | EXTERNAL nlte_tcool |
|---|
| 263 | EXTERNAL nlte_setup |
|---|
| 264 | EXTERNAL blendrad |
|---|
| 265 | EXTERNAL nlthermeq |
|---|
| 266 | EXTERNAL euvheat |
|---|
| 267 | EXTERNAL param_read_e107 |
|---|
| 268 | EXTERNAL conduction |
|---|
| 269 | EXTERNAL molvis |
|---|
| 270 | EXTERNAL moldiff_red |
|---|
| 271 | |
|---|
| 272 | c |
|---|
| 273 | c Variables locales |
|---|
| 274 | c |
|---|
| 275 | CXXX PB |
|---|
| 276 | REAL fluxt(klon,klev) ! flux turbulent de chaleur |
|---|
| 277 | REAL fluxu(klon,klev) ! flux turbulent de vitesse u |
|---|
| 278 | REAL fluxv(klon,klev) ! flux turbulent de vitesse v |
|---|
| 279 | c |
|---|
| 280 | REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique |
|---|
| 281 | REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec |
|---|
| 282 | REAL flux_ec(klon,klev) ! flux de chaleur Ec |
|---|
| 283 | c |
|---|
| 284 | REAL tmpout(klon,klev) ! [K/s] |
|---|
| 285 | c |
|---|
| 286 | REAL dist, rmu0(klon), fract(klon) |
|---|
| 287 | REAL zdtime, zlongi |
|---|
| 288 | c |
|---|
| 289 | INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil |
|---|
| 290 | c |
|---|
| 291 | REAL zphi(klon,klev) |
|---|
| 292 | REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2 |
|---|
| 293 | REAL tlaymean ! valeur temporaire pour calculer zzlay |
|---|
| 294 | real tsurf(klon) |
|---|
| 295 | |
|---|
| 296 | c va avec nlte_tcool |
|---|
| 297 | INTEGER ierr_nlte |
|---|
| 298 | REAL varerr |
|---|
| 299 | |
|---|
| 300 | ! photochemistry |
|---|
| 301 | |
|---|
| 302 | integer :: chempas |
|---|
| 303 | real :: zctime |
|---|
| 304 | |
|---|
| 305 | ! sedimentation |
|---|
| 306 | |
|---|
| 307 | REAL :: m0_mode1(klev,2),m0_mode2(klev,2) |
|---|
| 308 | REAL :: m3_mode1(klev,3),m3_mode2 (klev,3) |
|---|
| 309 | REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2) |
|---|
| 310 | REAL :: aer_flux(klev) |
|---|
| 311 | c |
|---|
| 312 | c Variables du changement |
|---|
| 313 | c |
|---|
| 314 | c ajs: ajustement sec |
|---|
| 315 | c vdf: couche limite (Vertical DiFfusion) |
|---|
| 316 | REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax) |
|---|
| 317 | REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev) |
|---|
| 318 | c |
|---|
| 319 | REAL d_ts(klon) |
|---|
| 320 | c |
|---|
| 321 | REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev) |
|---|
| 322 | REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax) |
|---|
| 323 | c |
|---|
| 324 | CMOD LOTT: Tendances Orography Sous-maille |
|---|
| 325 | REAL d_u_oro(klon,klev), d_v_oro(klon,klev) |
|---|
| 326 | REAL d_t_oro(klon,klev) |
|---|
| 327 | REAL d_u_lif(klon,klev), d_v_lif(klon,klev) |
|---|
| 328 | REAL d_t_lif(klon,klev) |
|---|
| 329 | C Tendances Ondes de G non oro (runs strato). |
|---|
| 330 | REAL d_u_hin(klon,klev), d_v_hin(klon,klev) |
|---|
| 331 | REAL d_t_hin(klon,klev) |
|---|
| 332 | |
|---|
| 333 | c Tendencies due to radiative scheme [K/s] |
|---|
| 334 | c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv |
|---|
| 335 | c are not computed at each physical timestep |
|---|
| 336 | c therefore, they are defined and saved in phys_state_var_mod |
|---|
| 337 | |
|---|
| 338 | c Tendencies due to molecular viscosity and conduction |
|---|
| 339 | real d_t_conduc(klon,klev) ! [K/s] |
|---|
| 340 | real d_u_molvis(klon,klev) ! [m/s] /s |
|---|
| 341 | real d_v_molvis(klon,klev) ! [m/s] /s |
|---|
| 342 | |
|---|
| 343 | c Tendencies due to molecular diffusion |
|---|
| 344 | real d_q_moldif(klon,klev,nqmax) |
|---|
| 345 | |
|---|
| 346 | c Tendencies due to ambipolar ion diffusion |
|---|
| 347 | real d_q_iondif(klon,klev,nqmax) |
|---|
| 348 | |
|---|
| 349 | c |
|---|
| 350 | c Variables liees a l'ecriture de la bande histoire physique |
|---|
| 351 | c |
|---|
| 352 | INTEGER ecrit_mth |
|---|
| 353 | SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel) |
|---|
| 354 | c |
|---|
| 355 | INTEGER ecrit_day |
|---|
| 356 | SAVE ecrit_day ! frequence d'ecriture (fichier journalier) |
|---|
| 357 | c |
|---|
| 358 | INTEGER ecrit_ins |
|---|
| 359 | SAVE ecrit_ins ! frequence d'ecriture (fichier instantane) |
|---|
| 360 | c |
|---|
| 361 | integer itau_w ! pas de temps ecriture = itap + itau_phy |
|---|
| 362 | |
|---|
| 363 | c Variables locales pour effectuer les appels en serie |
|---|
| 364 | c |
|---|
| 365 | REAL t_seri(klon,klev) |
|---|
| 366 | REAL u_seri(klon,klev), v_seri(klon,klev) |
|---|
| 367 | c |
|---|
| 368 | REAL :: tr_seri(klon,klev,nqmax) |
|---|
| 369 | REAL :: tr_hedin(klon,klev,nqmax) |
|---|
| 370 | REAL :: d_tr(klon,klev,nqmax) |
|---|
| 371 | c pour sorties |
|---|
| 372 | REAL :: col_dens_tr(klon,nqmax) |
|---|
| 373 | REAL,allocatable,save :: prod_tr(:,:,:) |
|---|
| 374 | REAL,allocatable,save :: loss_tr(:,:,:) |
|---|
| 375 | |
|---|
| 376 | c pour ioipsl |
|---|
| 377 | INTEGER nid_day, nid_mth, nid_ins |
|---|
| 378 | SAVE nid_day, nid_mth, nid_ins |
|---|
| 379 | INTEGER nhori, nvert, idayref |
|---|
| 380 | REAL zsto, zout, zsto1, zsto2, zero |
|---|
| 381 | parameter (zero=0.0e0) |
|---|
| 382 | real zjulian |
|---|
| 383 | save zjulian |
|---|
| 384 | |
|---|
| 385 | CHARACTER*2 str2 |
|---|
| 386 | character*20 modname |
|---|
| 387 | character*80 abort_message |
|---|
| 388 | logical ok_sync |
|---|
| 389 | |
|---|
| 390 | character*30 nom_fichier |
|---|
| 391 | character*10 varname |
|---|
| 392 | character*40 vartitle |
|---|
| 393 | character*20 varunits |
|---|
| 394 | C Variables liees au bilan d'energie et d'enthalpi |
|---|
| 395 | REAL ztsol(klon) |
|---|
| 396 | REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot |
|---|
| 397 | $ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot |
|---|
| 398 | SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot |
|---|
| 399 | $ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot |
|---|
| 400 | REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec |
|---|
| 401 | REAL d_h_vcol_phy |
|---|
| 402 | REAL fs_bound, fq_bound |
|---|
| 403 | SAVE d_h_vcol_phy |
|---|
| 404 | REAL zero_v(klon),zero_v2(klon,klev) |
|---|
| 405 | CHARACTER*15 ztit |
|---|
| 406 | INTEGER ip_ebil ! PRINT level for energy conserv. diag. |
|---|
| 407 | SAVE ip_ebil |
|---|
| 408 | DATA ip_ebil/2/ |
|---|
| 409 | INTEGER if_ebil ! level for energy conserv. dignostics |
|---|
| 410 | SAVE if_ebil |
|---|
| 411 | c+jld ec_conser |
|---|
| 412 | REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique |
|---|
| 413 | c-jld ec_conser |
|---|
| 414 | |
|---|
| 415 | c ALBEDO VARIATIONS (VCD) |
|---|
| 416 | REAL factAlb |
|---|
| 417 | c TEST VENUS... |
|---|
| 418 | REAL mang(klon,klev) ! moment cinetique |
|---|
| 419 | REAL mangtot ! moment cinetique total |
|---|
| 420 | |
|---|
| 421 | c cell_area for outputs in hist* |
|---|
| 422 | REAL cell_area_out(klon) |
|---|
| 423 | #ifdef MESOSCALE |
|---|
| 424 | REAL :: dt_dyn(klev) |
|---|
| 425 | #endif |
|---|
| 426 | c Declaration des constantes et des fonctions thermodynamiques |
|---|
| 427 | c |
|---|
| 428 | #include "YOMCST.h" |
|---|
| 429 | |
|---|
| 430 | c====================================================================== |
|---|
| 431 | c====================================================================== |
|---|
| 432 | c INITIALISATIONS |
|---|
| 433 | c====================================================================== |
|---|
| 434 | |
|---|
| 435 | modname = 'physiq' |
|---|
| 436 | ok_sync=.TRUE. |
|---|
| 437 | |
|---|
| 438 | bilansmc = 0 |
|---|
| 439 | ballons = 0 |
|---|
| 440 | ! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!! |
|---|
| 441 | #ifndef MESOSCALE |
|---|
| 442 | if (is_parallel) then |
|---|
| 443 | bilansmc = 0 |
|---|
| 444 | ballons = 0 |
|---|
| 445 | endif |
|---|
| 446 | #endif |
|---|
| 447 | IF (if_ebil.ge.1) THEN |
|---|
| 448 | DO i=1,klon |
|---|
| 449 | zero_v(i)=0. |
|---|
| 450 | END DO |
|---|
| 451 | DO i=1,klon |
|---|
| 452 | DO j=1,klev |
|---|
| 453 | zero_v2(i,j)=0. |
|---|
| 454 | END DO |
|---|
| 455 | END DO |
|---|
| 456 | END IF |
|---|
| 457 | |
|---|
| 458 | c====================================================================== |
|---|
| 459 | c DONE ONLY AT FIRST CALL |
|---|
| 460 | c======================== |
|---|
| 461 | IF (debut) THEN |
|---|
| 462 | allocate(source(klon,nqmax)) |
|---|
| 463 | allocate(prod_tr(klon,klev,nqmax)) |
|---|
| 464 | allocate(loss_tr(klon,klev,nqmax)) |
|---|
| 465 | allocate(no_emission(klon,klev)) |
|---|
| 466 | allocate(o2_emission(klon,klev)) |
|---|
| 467 | |
|---|
| 468 | #ifdef CPP_XIOS |
|---|
| 469 | ! Initialize XIOS context |
|---|
| 470 | write(*,*) "physiq: call wxios_context_init" |
|---|
| 471 | CALL wxios_context_init |
|---|
| 472 | #endif |
|---|
| 473 | |
|---|
| 474 | ! The call to suphec is now done in iniphysiq_mod (interface) |
|---|
| 475 | ! CALL suphec ! initialiser constantes et parametres phys. |
|---|
| 476 | |
|---|
| 477 | IF (if_ebil.ge.1) d_h_vcol_phy=0. |
|---|
| 478 | ! |
|---|
| 479 | ! load flag and parameter values from physiq.def |
|---|
| 480 | ! |
|---|
| 481 | call conf_phys(ok_journe, ok_mensuel, |
|---|
| 482 | . ok_instan, |
|---|
| 483 | . if_ebil) |
|---|
| 484 | |
|---|
| 485 | call phys_state_var_init(nqmax) |
|---|
| 486 | c |
|---|
| 487 | c Initialising Hedin model for upper atm |
|---|
| 488 | c (to be revised when coupled to chemistry) : |
|---|
| 489 | call conc_init |
|---|
| 490 | |
|---|
| 491 | ! initialise physics counter |
|---|
| 492 | |
|---|
| 493 | itap = 0 |
|---|
| 494 | |
|---|
| 495 | #ifdef MESOSCALE |
|---|
| 496 | print*,'check pdtphys',pdtphys |
|---|
| 497 | PRINT*,'check phisfi ',pphis(1),pphis(klon) |
|---|
| 498 | PRINT*,'check geop',pphi(1,1),pphi(klon,klev) |
|---|
| 499 | PRINT*,'check radsol',radsol(1),radsol(klon) |
|---|
| 500 | print*,'check ppk',ppk(1,1),ppk(klon,klev) |
|---|
| 501 | print*,'check ftsoil',ftsoil(1,1),ftsoil(klon,nsoilmx) |
|---|
| 502 | print*,'check ftsol',ftsol(1),ftsol(klon) |
|---|
| 503 | print*, "check temp", t(1,1),t(klon,klev) |
|---|
| 504 | print*, "check pres",paprs(1,1),paprs(klon,klev),pplay(1,1), |
|---|
| 505 | . pplay(klon,klev) |
|---|
| 506 | print*, "check u", u(1,1),u(klon,klev) |
|---|
| 507 | print*, "check v", v(1,1),v(klon,klev) |
|---|
| 508 | print*,'check falbe',falbe(1),falbe(klon) |
|---|
| 509 | !nqtot=nqmax |
|---|
| 510 | !ALLOCATE(tname(nqtot)) |
|---|
| 511 | !tname=noms |
|---|
| 512 | zmea=0. |
|---|
| 513 | zstd=0. |
|---|
| 514 | zsig=0. |
|---|
| 515 | zgam=0. |
|---|
| 516 | zthe=0. |
|---|
| 517 | dtime=pdtphys |
|---|
| 518 | #else |
|---|
| 519 | c |
|---|
| 520 | c Lecture startphy.nc : |
|---|
| 521 | c |
|---|
| 522 | CALL phyetat0 ("startphy.nc") |
|---|
| 523 | IF (.not.startphy_file) THEN |
|---|
| 524 | ! Additionnal academic initializations |
|---|
| 525 | ftsol(:)=t(:,1) ! surface temperature as in first atm. layer |
|---|
| 526 | DO isoil=1, nsoilmx |
|---|
| 527 | ! subsurface temperatures equal to surface temperature |
|---|
| 528 | ftsoil(:,isoil)=ftsol(:) |
|---|
| 529 | ENDDO |
|---|
| 530 | ENDIF |
|---|
| 531 | #endif |
|---|
| 532 | |
|---|
| 533 | c dtime est defini dans tabcontrol.h et lu dans startphy |
|---|
| 534 | c pdtphys est calcule a partir des nouvelles conditions: |
|---|
| 535 | c Reinitialisation du pas de temps physique quand changement |
|---|
| 536 | IF (ABS(dtime-pdtphys).GT.0.001) THEN |
|---|
| 537 | WRITE(lunout,*) 'Pas physique a change',dtime, |
|---|
| 538 | . pdtphys |
|---|
| 539 | c abort_message='Pas physique n est pas correct ' |
|---|
| 540 | c call abort_physic(modname,abort_message,1) |
|---|
| 541 | c---------------- |
|---|
| 542 | c pour initialiser convenablement le time_counter, il faut tenir compte |
|---|
| 543 | c du changement de dtime en changeant itau_phy (point de depart) |
|---|
| 544 | itau_phy = NINT(itau_phy*dtime/pdtphys) |
|---|
| 545 | c---------------- |
|---|
| 546 | dtime=pdtphys |
|---|
| 547 | ENDIF |
|---|
| 548 | |
|---|
| 549 | radpas = NINT(RDAY/pdtphys/nbapp_rad) |
|---|
| 550 | |
|---|
| 551 | CALL printflag( ok_journe,ok_instan ) |
|---|
| 552 | |
|---|
| 553 | #ifdef CPP_XIOS |
|---|
| 554 | write(*,*) "physiq: call initialize_xios_output" |
|---|
| 555 | call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY, |
|---|
| 556 | & presnivs,pseudoalt) |
|---|
| 557 | #endif |
|---|
| 558 | |
|---|
| 559 | c |
|---|
| 560 | c--------- |
|---|
| 561 | c FLOTT |
|---|
| 562 | IF (ok_orodr) THEN |
|---|
| 563 | DO i=1,klon |
|---|
| 564 | rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0) |
|---|
| 565 | ENDDO |
|---|
| 566 | CALL SUGWD(klon,klev,paprs,pplay) |
|---|
| 567 | DO i=1,klon |
|---|
| 568 | zuthe(i)=0. |
|---|
| 569 | zvthe(i)=0. |
|---|
| 570 | if(zstd(i).gt.10.)then |
|---|
| 571 | zuthe(i)=(1.-zgam(i))*cos(zthe(i)) |
|---|
| 572 | zvthe(i)=(1.-zgam(i))*sin(zthe(i)) |
|---|
| 573 | endif |
|---|
| 574 | ENDDO |
|---|
| 575 | ENDIF |
|---|
| 576 | |
|---|
| 577 | if (bilansmc.eq.1) then |
|---|
| 578 | C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES |
|---|
| 579 | C DU BILAN DE MOMENT ANGULAIRE. |
|---|
| 580 | open(27,file='aaam_bud.out',form='formatted') |
|---|
| 581 | open(28,file='fields_2d.out',form='formatted') |
|---|
| 582 | write(*,*)'Ouverture de aaam_bud.out (FL Vous parle)' |
|---|
| 583 | write(*,*)'Ouverture de fields_2d.out (FL Vous parle)' |
|---|
| 584 | endif !bilansmc |
|---|
| 585 | |
|---|
| 586 | c--------------SLEBONNOIS |
|---|
| 587 | C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES |
|---|
| 588 | C DES BALLONS |
|---|
| 589 | if (ballons.eq.1) then |
|---|
| 590 | open(30,file='ballons-lat.out',form='formatted') |
|---|
| 591 | open(31,file='ballons-lon.out',form='formatted') |
|---|
| 592 | open(32,file='ballons-u.out',form='formatted') |
|---|
| 593 | open(33,file='ballons-v.out',form='formatted') |
|---|
| 594 | open(34,file='ballons-alt.out',form='formatted') |
|---|
| 595 | write(*,*)'Ouverture des ballons*.out' |
|---|
| 596 | endif !ballons |
|---|
| 597 | c------------- |
|---|
| 598 | |
|---|
| 599 | c--------- |
|---|
| 600 | C TRACEURS |
|---|
| 601 | C source dans couche limite |
|---|
| 602 | source(:,:) = 0.0 ! pas de source, pour l'instant |
|---|
| 603 | c--------- |
|---|
| 604 | |
|---|
| 605 | c--------- |
|---|
| 606 | c INITIALIZE THERMOSPHERIC PARAMETERS |
|---|
| 607 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 608 | |
|---|
| 609 | if (callthermos) then |
|---|
| 610 | call fill_data_thermos |
|---|
| 611 | call allocate_param_thermos(klev) |
|---|
| 612 | call param_read_e107 |
|---|
| 613 | endif |
|---|
| 614 | |
|---|
| 615 | c Initialisation (recomputed in concentration2) |
|---|
| 616 | do ig=1,klon |
|---|
| 617 | do j=1,klev |
|---|
| 618 | rnew(ig,j)=RD |
|---|
| 619 | cpnew(ig,j)=cpdet(t(ig,j)) |
|---|
| 620 | mmean(ig,j)=RMD |
|---|
| 621 | akknew(ig,j)=1.e-4 |
|---|
| 622 | rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j)) |
|---|
| 623 | enddo |
|---|
| 624 | |
|---|
| 625 | enddo |
|---|
| 626 | |
|---|
| 627 | IF(callthermos.or.callnlte.or.callnirco2) THEN |
|---|
| 628 | call compo_hedin83_init2 |
|---|
| 629 | ENDIF |
|---|
| 630 | if (callnlte.and.nltemodel.eq.2) call nlte_setup |
|---|
| 631 | if (callnirco2.and.nircorr.eq.1) call nir_leedat |
|---|
| 632 | c--------- |
|---|
| 633 | |
|---|
| 634 | c |
|---|
| 635 | c Verifications: |
|---|
| 636 | c |
|---|
| 637 | IF (nlon .NE. klon) THEN |
|---|
| 638 | WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, |
|---|
| 639 | . klon |
|---|
| 640 | abort_message='nlon et klon ne sont pas coherents' |
|---|
| 641 | call abort_physic(modname,abort_message,1) |
|---|
| 642 | ENDIF |
|---|
| 643 | IF (nlev .NE. klev) THEN |
|---|
| 644 | WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, |
|---|
| 645 | . klev |
|---|
| 646 | abort_message='nlev et klev ne sont pas coherents' |
|---|
| 647 | call abort_physic(modname,abort_message,1) |
|---|
| 648 | ENDIF |
|---|
| 649 | c |
|---|
| 650 | IF (dtime*REAL(radpas).GT.(RDAY*0.25).AND.cycle_diurne) |
|---|
| 651 | $ THEN |
|---|
| 652 | WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant' |
|---|
| 653 | WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne" |
|---|
| 654 | abort_message='Nbre d appels au rayonnement insuffisant' |
|---|
| 655 | call abort_physic(modname,abort_message,1) |
|---|
| 656 | ENDIF |
|---|
| 657 | c |
|---|
| 658 | WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=", |
|---|
| 659 | . iflag_ajs |
|---|
| 660 | c |
|---|
| 661 | ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours |
|---|
| 662 | IF (ok_mensuel) THEN |
|---|
| 663 | WRITE(lunout,*)'La frequence de sortie mensuelle est de ', |
|---|
| 664 | . ecrit_mth |
|---|
| 665 | ENDIF |
|---|
| 666 | |
|---|
| 667 | ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours |
|---|
| 668 | IF (ok_journe) THEN |
|---|
| 669 | WRITE(lunout,*)'La frequence de sortie journaliere est de ', |
|---|
| 670 | . ecrit_day |
|---|
| 671 | ENDIF |
|---|
| 672 | |
|---|
| 673 | ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable |
|---|
| 674 | IF (ok_instan) THEN |
|---|
| 675 | WRITE(lunout,*)'La frequence de sortie instant. est de ', |
|---|
| 676 | . ecrit_ins |
|---|
| 677 | ENDIF |
|---|
| 678 | |
|---|
| 679 | c Verification synchronize AMBIPOLAR DIFFUSION & CHEMISTRY |
|---|
| 680 | |
|---|
| 681 | IF ((ok_iondiff) .and. (NINT(RDAY/dtime).ne.nbapp_chem)) THEN |
|---|
| 682 | WRITE(lunout,*)'nbapp_chem .NE. day_step' |
|---|
| 683 | WRITE(lunout,*)'nbapp_chem ', nbapp_chem |
|---|
| 684 | WRITE(lunout,*)'day_step ', NINT(RDAY/dtime) |
|---|
| 685 | WRITE(lunout,*)'nbapp_chem must be equal to day_step' |
|---|
| 686 | STOP |
|---|
| 687 | ENDIF |
|---|
| 688 | |
|---|
| 689 | c Initialisation des sorties |
|---|
| 690 | c=========================== |
|---|
| 691 | |
|---|
| 692 | #ifdef CPP_IOIPSL |
|---|
| 693 | |
|---|
| 694 | #ifdef histhf |
|---|
| 695 | #include "ini_histhf.h" |
|---|
| 696 | #endif |
|---|
| 697 | |
|---|
| 698 | #ifdef histday |
|---|
| 699 | #include "ini_histday.h" |
|---|
| 700 | #endif |
|---|
| 701 | |
|---|
| 702 | #ifdef histmth |
|---|
| 703 | #include "ini_histmth.h" |
|---|
| 704 | #endif |
|---|
| 705 | |
|---|
| 706 | #ifdef histins |
|---|
| 707 | #include "ini_histins.h" |
|---|
| 708 | #endif |
|---|
| 709 | |
|---|
| 710 | #endif |
|---|
| 711 | |
|---|
| 712 | c |
|---|
| 713 | c Initialiser les valeurs de u pour calculs tendances |
|---|
| 714 | c (pour T, c'est fait dans phyetat0) |
|---|
| 715 | c |
|---|
| 716 | DO k = 1, klev |
|---|
| 717 | DO i = 1, klon |
|---|
| 718 | u_ancien(i,k) = u(i,k) |
|---|
| 719 | ENDDO |
|---|
| 720 | ENDDO |
|---|
| 721 | |
|---|
| 722 | ! initialisation of microphysical and chemical parameters |
|---|
| 723 | |
|---|
| 724 | if (ok_chem .and. .not. ok_cloud) then |
|---|
| 725 | print*,"chemistry requires clouds" |
|---|
| 726 | print*,"ok_cloud must be .true." |
|---|
| 727 | stop |
|---|
| 728 | end if |
|---|
| 729 | |
|---|
| 730 | if (.not. ok_chem .and. ok_cloud .and. (cl_scheme == 1)) then |
|---|
| 731 | print*,"cl_scheme = 1 requires chemistry" |
|---|
| 732 | print*,"ok_chem must be .true." |
|---|
| 733 | stop |
|---|
| 734 | end if |
|---|
| 735 | |
|---|
| 736 | ! number of microphysical tracers |
|---|
| 737 | |
|---|
| 738 | nmicro = 0 |
|---|
| 739 | if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2 |
|---|
| 740 | if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12 |
|---|
| 741 | |
|---|
| 742 | ! initialise chemical parameters. includes the indexation of microphys tracers |
|---|
| 743 | |
|---|
| 744 | if (ok_chem .or. cl_scheme == 2) then |
|---|
| 745 | call chemparam_ini() |
|---|
| 746 | end if |
|---|
| 747 | |
|---|
| 748 | if ((iflag_trac.eq.1) .and. ok_aoa) then |
|---|
| 749 | write(lunout,*) 'Initialising age of air subroutine' |
|---|
| 750 | allocate(init_age(klon,klev)) |
|---|
| 751 | call aoa_ini(ok_chem, tr_scheme) ! get index of age of air tracer in the tracer array |
|---|
| 752 | if (reinit_aoa) then |
|---|
| 753 | age(:,:) = 0. ! Set initial value of age of air to 0 everywhere |
|---|
| 754 | tr_seri(:,:,i_aoa) = 1e-30 ! Set initial value of tracer to tiny everywhere |
|---|
| 755 | endif ! reinitialisation loop |
|---|
| 756 | init_age(:,:) = age(:,:) ! save initial value of age, either read in from restartphy or set to 0 |
|---|
| 757 | endif ! age of air initialisation |
|---|
| 758 | |
|---|
| 759 | ! initialise cloud parameters (for cl_scheme = 1) |
|---|
| 760 | |
|---|
| 761 | if (ok_cloud .and. (cl_scheme == 1)) then |
|---|
| 762 | call cloud_ini(nlon,nlev,nb_mode) |
|---|
| 763 | end if |
|---|
| 764 | |
|---|
| 765 | ! initialise mmean |
|---|
| 766 | |
|---|
| 767 | if (callthermos) then |
|---|
| 768 | call concentrations2(pplay,t,qx,nqmax) |
|---|
| 769 | end if |
|---|
| 770 | |
|---|
| 771 | c======================== |
|---|
| 772 | ENDIF ! debut |
|---|
| 773 | c======================== |
|---|
| 774 | |
|---|
| 775 | c====================================================================== |
|---|
| 776 | ! ------------------------------------------------------ |
|---|
| 777 | ! Initializations done at every physical timestep: |
|---|
| 778 | ! ------------------------------------------------------ |
|---|
| 779 | |
|---|
| 780 | c Mettre a zero des variables de sortie (pour securite) |
|---|
| 781 | c |
|---|
| 782 | DO i = 1, klon |
|---|
| 783 | d_ps(i) = 0.0 |
|---|
| 784 | ENDDO |
|---|
| 785 | DO k = 1, klev |
|---|
| 786 | DO i = 1, klon |
|---|
| 787 | d_t(i,k) = 0.0 |
|---|
| 788 | d_u(i,k) = 0.0 |
|---|
| 789 | d_v(i,k) = 0.0 |
|---|
| 790 | ENDDO |
|---|
| 791 | ENDDO |
|---|
| 792 | DO iq = 1, nqmax |
|---|
| 793 | DO k = 1, klev |
|---|
| 794 | DO i = 1, klon |
|---|
| 795 | d_qx(i,k,iq) = 0.0 |
|---|
| 796 | ENDDO |
|---|
| 797 | ENDDO |
|---|
| 798 | ENDDO |
|---|
| 799 | c |
|---|
| 800 | c Ne pas affecter les valeurs entrees de u, v, h, et q |
|---|
| 801 | c |
|---|
| 802 | DO k = 1, klev |
|---|
| 803 | DO i = 1, klon |
|---|
| 804 | t_seri(i,k) = t(i,k) |
|---|
| 805 | u_seri(i,k) = u(i,k) |
|---|
| 806 | v_seri(i,k) = v(i,k) |
|---|
| 807 | ENDDO |
|---|
| 808 | ENDDO |
|---|
| 809 | DO iq = 1, nqmax |
|---|
| 810 | DO k = 1, klev |
|---|
| 811 | DO i = 1, klon |
|---|
| 812 | tr_seri(i,k,iq) = qx(i,k,iq) |
|---|
| 813 | ENDDO |
|---|
| 814 | ENDDO |
|---|
| 815 | ENDDO |
|---|
| 816 | C |
|---|
| 817 | DO i = 1, klon |
|---|
| 818 | ztsol(i) = ftsol(i) |
|---|
| 819 | ENDDO |
|---|
| 820 | C |
|---|
| 821 | IF (if_ebil.ge.1) THEN |
|---|
| 822 | ztit='after dynamic' |
|---|
| 823 | CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime |
|---|
| 824 | e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay |
|---|
| 825 | s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) |
|---|
| 826 | C Comme les tendances de la physique sont ajoute dans la dynamique, |
|---|
| 827 | C on devrait avoir que la variation d'entalpie par la dynamique |
|---|
| 828 | C est egale a la variation de la physique au pas de temps precedent. |
|---|
| 829 | C Donc la somme de ces 2 variations devrait etre nulle. |
|---|
| 830 | call diagphy(cell_area,ztit,ip_ebil |
|---|
| 831 | e , zero_v, zero_v, zero_v, zero_v, zero_v |
|---|
| 832 | e , zero_v, zero_v, zero_v, ztsol |
|---|
| 833 | e , d_h_vcol+d_h_vcol_phy, d_qt, 0. |
|---|
| 834 | s , fs_bound, fq_bound ) |
|---|
| 835 | END IF |
|---|
| 836 | |
|---|
| 837 | c==================================================================== |
|---|
| 838 | c XIOS outputs |
|---|
| 839 | |
|---|
| 840 | #ifdef CPP_XIOS |
|---|
| 841 | ! update XIOS time/calendar |
|---|
| 842 | call update_xios_timestep |
|---|
| 843 | #endif |
|---|
| 844 | |
|---|
| 845 | c==================================================================== |
|---|
| 846 | c Diagnostiquer la tendance dynamique |
|---|
| 847 | c |
|---|
| 848 | IF (ancien_ok) THEN |
|---|
| 849 | DO k = 1, klev |
|---|
| 850 | DO i = 1, klon |
|---|
| 851 | d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime |
|---|
| 852 | d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime |
|---|
| 853 | ENDDO |
|---|
| 854 | ENDDO |
|---|
| 855 | |
|---|
| 856 | ! ADAPTATION GCM POUR CP(T) |
|---|
| 857 | do i=1,klon |
|---|
| 858 | flux_dyn(i,1) = 0.0 |
|---|
| 859 | do j=2,klev |
|---|
| 860 | flux_dyn(i,j) = flux_dyn(i,j-1) |
|---|
| 861 | . +cpnew(i,j-1)/RG*d_t_dyn(i,j-1)*(paprs(i,j-1)-paprs(i,j)) |
|---|
| 862 | enddo |
|---|
| 863 | enddo |
|---|
| 864 | |
|---|
| 865 | ELSE |
|---|
| 866 | DO k = 1, klev |
|---|
| 867 | DO i = 1, klon |
|---|
| 868 | d_u_dyn(i,k) = 0.0 |
|---|
| 869 | d_t_dyn(i,k) = 0.0 |
|---|
| 870 | ENDDO |
|---|
| 871 | ENDDO |
|---|
| 872 | ancien_ok = .TRUE. |
|---|
| 873 | ENDIF |
|---|
| 874 | c==================================================================== |
|---|
| 875 | |
|---|
| 876 | ! Compute vertical velocity (Pa/s) from vertical mass flux |
|---|
| 877 | ! Need to linearly interpolate mass flux to mid-layers |
|---|
| 878 | do k=1,klev-1 |
|---|
| 879 | omega(1:klon,k) = 0.5*RG*(flxmw(1:klon,k)+flxmw(1:klon,k+1)) |
|---|
| 880 | . / cell_area(1:klon) |
|---|
| 881 | enddo |
|---|
| 882 | omega(1:klon,klev) = 0.5*RG*flxmw(1:klon,klev) / cell_area(1:klon) |
|---|
| 883 | |
|---|
| 884 | c====== |
|---|
| 885 | c GEOP CORRECTION |
|---|
| 886 | c |
|---|
| 887 | c Ajouter le geopotentiel du sol: |
|---|
| 888 | c |
|---|
| 889 | DO k = 1, klev |
|---|
| 890 | DO i = 1, klon |
|---|
| 891 | zphi(i,k) = pphi(i,k) + pphis(i) |
|---|
| 892 | ENDDO |
|---|
| 893 | ENDDO |
|---|
| 894 | |
|---|
| 895 | c............................ |
|---|
| 896 | c CETTE CORRECTION VA DE PAIR AVEC DES MODIFS DE LEAPFROG(_p) |
|---|
| 897 | c ELLE MARCHE A 50 NIVEAUX (si mmean constante...) |
|---|
| 898 | c MAIS PAS A 78 NIVEAUX (quand mmean varie...) |
|---|
| 899 | c A ANALYSER PLUS EN DETAIL AVANT D'UTILISER |
|---|
| 900 | c............................ |
|---|
| 901 | c zphi is recomputed (pphi is not ok if mean molecular mass varies) |
|---|
| 902 | c with dphi = RT/mmean d(ln p) [evaluated at interface] |
|---|
| 903 | |
|---|
| 904 | c DO i = 1, klon |
|---|
| 905 | c zphi(i,1) = pphis(i) + R*t_seri(i,1)/mmean(i,1)*1000. |
|---|
| 906 | c * *( log(paprs(i,1)) - log(pplay(i,1)) ) |
|---|
| 907 | c DO k = 2, klev |
|---|
| 908 | c zphi(i,k) = zphi(i,k-1) |
|---|
| 909 | c * + R*500.*(t_seri(i,k)/mmean(i,k)+t_seri(i,k-1)/mmean(i,k-1)) |
|---|
| 910 | c * * (log(pplay(i,k-1)) - log(pplay(i,k))) |
|---|
| 911 | c ENDDO |
|---|
| 912 | c ENDDO |
|---|
| 913 | c............................ |
|---|
| 914 | c===== |
|---|
| 915 | |
|---|
| 916 | c calcul de l'altitude aux niveaux intercouches |
|---|
| 917 | c ponderation des altitudes au niveau des couches en dp/p |
|---|
| 918 | |
|---|
| 919 | DO i=1,klon |
|---|
| 920 | zzlay(i,1)=zphi(i,1)/RG ! [m] |
|---|
| 921 | zzlev(i,1)=pphis(i)/RG ! [m] |
|---|
| 922 | ENDDO |
|---|
| 923 | DO k=2,klev |
|---|
| 924 | DO i=1,klon |
|---|
| 925 | tlaymean=t_seri(i,k) |
|---|
| 926 | IF (t_seri(i,k).ne.t_seri(i,k-1)) |
|---|
| 927 | & tlaymean=(t_seri(i,k)-t_seri(i,k-1)) |
|---|
| 928 | & /log(t_seri(i,k)/t_seri(i,k-1)) |
|---|
| 929 | |
|---|
| 930 | zzlay(i,k)=zzlay(i,k-1) |
|---|
| 931 | & -(log(pplay(i,k)/pplay(i,k-1))*rnew(i,k-1)*tlaymean |
|---|
| 932 | & /(RG*(RA/(RA+zzlay(i,k-1)))**2)) |
|---|
| 933 | ENDDO |
|---|
| 934 | ENDDO |
|---|
| 935 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 936 | ! Old : from geopotential. Problem when mu varies in the upper atm... |
|---|
| 937 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 938 | ! DO k=1,klev |
|---|
| 939 | ! DO i=1,klon |
|---|
| 940 | ! zzlay(i,k)=zphi(i,k)/RG ! [m] |
|---|
| 941 | ! ENDDO |
|---|
| 942 | ! ENDDO |
|---|
| 943 | ! DO i=1,klon |
|---|
| 944 | ! zzlev(i,1)=pphis(i)/RG ! [m] |
|---|
| 945 | ! ENDDO |
|---|
| 946 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 947 | DO k=2,klev |
|---|
| 948 | DO i=1,klon |
|---|
| 949 | z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k)) |
|---|
| 950 | z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k)) |
|---|
| 951 | zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2) |
|---|
| 952 | ENDDO |
|---|
| 953 | ENDDO |
|---|
| 954 | DO i=1,klon |
|---|
| 955 | zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev)) |
|---|
| 956 | ENDDO |
|---|
| 957 | |
|---|
| 958 | c==================================================================== |
|---|
| 959 | c |
|---|
| 960 | c Verifier les temperatures |
|---|
| 961 | c |
|---|
| 962 | CALL hgardfou(t_seri,ftsol,'debutphy') |
|---|
| 963 | |
|---|
| 964 | c==================================================================== |
|---|
| 965 | c Orbite et eclairement |
|---|
| 966 | c======================= |
|---|
| 967 | |
|---|
| 968 | c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0. |
|---|
| 969 | c donc pas de variations de Ls, ni de dist. |
|---|
| 970 | c La seule chose qui compte, c'est la rotation de la planete devant |
|---|
| 971 | c le Soleil... |
|---|
| 972 | |
|---|
| 973 | zlongi = 0.0 |
|---|
| 974 | dist = 0.72333 ! en UA |
|---|
| 975 | |
|---|
| 976 | c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite |
|---|
| 977 | c a sa valeur, et prendre en compte leur evolution, |
|---|
| 978 | c il faudra refaire un orbite.F... |
|---|
| 979 | c CALL orbite(zlongi,dist) |
|---|
| 980 | |
|---|
| 981 | IF (cycle_diurne) THEN |
|---|
| 982 | zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s) |
|---|
| 983 | CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg, |
|---|
| 984 | & rmu0,fract) |
|---|
| 985 | ELSE |
|---|
| 986 | call mucorr(klon,zlongi,latitude_deg,rmu0,fract) |
|---|
| 987 | ENDIF |
|---|
| 988 | |
|---|
| 989 | ! print fraction of venus day |
|---|
| 990 | |
|---|
| 991 | if (is_master) then |
|---|
| 992 | print*, 'gmtime = ', gmtime |
|---|
| 993 | end if |
|---|
| 994 | |
|---|
| 995 | c====================================================================== |
|---|
| 996 | c FIN INITIALISATIONS |
|---|
| 997 | c====================================================================== |
|---|
| 998 | c====================================================================== |
|---|
| 999 | |
|---|
| 1000 | c======================================================= |
|---|
| 1001 | ! CONDUCTION and MOLECULAR VISCOSITY |
|---|
| 1002 | c======================================================= |
|---|
| 1003 | |
|---|
| 1004 | d_t_conduc(:,:)=0. |
|---|
| 1005 | d_u_molvis(:,:)=0. |
|---|
| 1006 | d_v_molvis(:,:)=0. |
|---|
| 1007 | |
|---|
| 1008 | IF (callthermos) THEN |
|---|
| 1009 | |
|---|
| 1010 | tsurf(:)=t_seri(:,1) |
|---|
| 1011 | call conduction(klon, klev,pdtphys, |
|---|
| 1012 | $ pplay,paprs,t_seri, |
|---|
| 1013 | $ tsurf,zzlev,zzlay,d_t_conduc) |
|---|
| 1014 | |
|---|
| 1015 | call molvis(klon, klev,pdtphys, |
|---|
| 1016 | $ pplay,paprs,t_seri, |
|---|
| 1017 | $ u,tsurf,zzlev,zzlay,d_u_molvis) |
|---|
| 1018 | |
|---|
| 1019 | call molvis(klon, klev, pdtphys, |
|---|
| 1020 | $ pplay,paprs,t_seri, |
|---|
| 1021 | $ v,tsurf,zzlev,zzlay,d_v_molvis) |
|---|
| 1022 | |
|---|
| 1023 | DO k=1,klev |
|---|
| 1024 | DO ig=1,klon |
|---|
| 1025 | t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K] |
|---|
| 1026 | u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! [m/s] |
|---|
| 1027 | v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! [m/s] |
|---|
| 1028 | ENDDO |
|---|
| 1029 | ENDDO |
|---|
| 1030 | ENDIF |
|---|
| 1031 | |
|---|
| 1032 | c==================================================================== |
|---|
| 1033 | c Compute mean mass, cp and R : |
|---|
| 1034 | c------------------------------------ |
|---|
| 1035 | |
|---|
| 1036 | if(callthermos) then |
|---|
| 1037 | call concentrations2(pplay,t_seri,tr_seri, nqmax) |
|---|
| 1038 | else |
|---|
| 1039 | ! one still needs to recompute cp and rho which depend on temperature |
|---|
| 1040 | ! (rnew, mmean and akknew remain already initialized constants) |
|---|
| 1041 | do ig=1,klon |
|---|
| 1042 | do j=1,klev |
|---|
| 1043 | cpnew(ig,j)=cpdet(t(ig,j)) |
|---|
| 1044 | rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j)) |
|---|
| 1045 | enddo |
|---|
| 1046 | enddo |
|---|
| 1047 | endif |
|---|
| 1048 | |
|---|
| 1049 | |
|---|
| 1050 | c======================================================= |
|---|
| 1051 | ! CHEMISTRY AND MICROPHYSICS |
|---|
| 1052 | c======================================================= |
|---|
| 1053 | |
|---|
| 1054 | if (iflag_trac.eq.1) then |
|---|
| 1055 | |
|---|
| 1056 | if ( ok_aoa ) then |
|---|
| 1057 | call age_of_air(tr_seri(:,:,i_aoa),klon,klev, |
|---|
| 1058 | I itap,pdtphys,init_age, |
|---|
| 1059 | O age) |
|---|
| 1060 | end if |
|---|
| 1061 | |
|---|
| 1062 | !==================================================================== |
|---|
| 1063 | ! Case 1: pseudo-chemistry with relaxation toward fixed profile |
|---|
| 1064 | !========= |
|---|
| 1065 | if (tr_scheme.eq.1) then |
|---|
| 1066 | |
|---|
| 1067 | call phytrac_relax (debut,lafin,nqmax, |
|---|
| 1068 | I nlon,nlev,dtime,pplay, |
|---|
| 1069 | O tr_seri) |
|---|
| 1070 | |
|---|
| 1071 | elseif (tr_scheme.eq.2) then |
|---|
| 1072 | !==================================================================== |
|---|
| 1073 | ! Case 2: surface emission |
|---|
| 1074 | ! For the moment, inspired from Mars version |
|---|
| 1075 | ! However, the variable 'source' could be used in physiq |
|---|
| 1076 | ! so the call to phytrac_emiss could be to initialise it. |
|---|
| 1077 | !========= |
|---|
| 1078 | |
|---|
| 1079 | call phytrac_emiss (debut,lafin,nqmax, |
|---|
| 1080 | I nlon,nlev,dtime,paprs, |
|---|
| 1081 | I latitude_deg,longitude_deg, |
|---|
| 1082 | O tr_seri) |
|---|
| 1083 | |
|---|
| 1084 | else if (tr_scheme.eq.3) then |
|---|
| 1085 | !==================================================================== |
|---|
| 1086 | ! Case 3: Full chemistry and/or clouds. |
|---|
| 1087 | ! routines are called every "chempas" physical timestep. |
|---|
| 1088 | ! |
|---|
| 1089 | ! if the physics is called 96000 times per venus day: |
|---|
| 1090 | ! |
|---|
| 1091 | ! nbapp_chem = 24000 => chempas = 4 => zctime = 420 s |
|---|
| 1092 | ! nbapp_chem = 12000 => chempas = 8 => zctime = 840 s |
|---|
| 1093 | !========= |
|---|
| 1094 | |
|---|
| 1095 | |
|---|
| 1096 | chempas = nint(rday/pdtphys/nbapp_chem) |
|---|
| 1097 | zctime = dtime*real(chempas) ! chemical timestep |
|---|
| 1098 | |
|---|
| 1099 | if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling |
|---|
| 1100 | |
|---|
| 1101 | ! photochemistry and microphysics |
|---|
| 1102 | |
|---|
| 1103 | call phytrac_chimie(debut, |
|---|
| 1104 | $ gmtime, |
|---|
| 1105 | $ nqmax, |
|---|
| 1106 | $ klon, |
|---|
| 1107 | $ latitude_deg, |
|---|
| 1108 | $ longitude_deg, |
|---|
| 1109 | $ zzlay, |
|---|
| 1110 | $ nlev, |
|---|
| 1111 | $ zctime, |
|---|
| 1112 | $ t_seri, |
|---|
| 1113 | $ pplay, |
|---|
| 1114 | $ tr_seri, |
|---|
| 1115 | $ d_tr_chem, |
|---|
| 1116 | $ iter, |
|---|
| 1117 | $ prod_tr, |
|---|
| 1118 | $ loss_tr, |
|---|
| 1119 | $ no_emission, |
|---|
| 1120 | $ o2_emission) |
|---|
| 1121 | |
|---|
| 1122 | if (ok_sedim) then |
|---|
| 1123 | if (cl_scheme == 1) then |
|---|
| 1124 | |
|---|
| 1125 | ! sedimentation for simplified microphysics |
|---|
| 1126 | |
|---|
| 1127 | #ifndef MESOSCALE |
|---|
| 1128 | call new_cloud_sedim(nb_mode, |
|---|
| 1129 | $ klon, |
|---|
| 1130 | $ nlev, |
|---|
| 1131 | $ zctime, |
|---|
| 1132 | $ pplay, |
|---|
| 1133 | $ paprs, |
|---|
| 1134 | $ t_seri, |
|---|
| 1135 | $ tr_seri, |
|---|
| 1136 | $ d_tr_chem, |
|---|
| 1137 | $ d_tr_sed(:,:,1:2), |
|---|
| 1138 | $ nqmax, |
|---|
| 1139 | $ Fsedim) |
|---|
| 1140 | |
|---|
| 1141 | ! test to avoid nans |
|---|
| 1142 | |
|---|
| 1143 | do k = 1, klev |
|---|
| 1144 | do i = 1, klon |
|---|
| 1145 | if ((d_tr_sed(i,k,1) /= d_tr_sed(i,k,1)) .or. |
|---|
| 1146 | $ (d_tr_sed(i,k,2) /= d_tr_sed(i,k,2))) then |
|---|
| 1147 | print*,'sedim NaN PROBLEM' |
|---|
| 1148 | print*,'d_tr_sed Nan?',d_tr_sed(i,k,:) |
|---|
| 1149 | print*,'Temp',t_seri(i,k) |
|---|
| 1150 | print*,'lat-lon',i,'level',k,'zctime',zctime |
|---|
| 1151 | print*,'F_sed',Fsedim(i,k) |
|---|
| 1152 | d_tr_sed(i,k,:) = 0. |
|---|
| 1153 | end if |
|---|
| 1154 | end do |
|---|
| 1155 | end do |
|---|
| 1156 | |
|---|
| 1157 | ! tendency due to condensation and sedimentation |
|---|
| 1158 | |
|---|
| 1159 | d_tr_sed(:,:,1:2) = d_tr_sed(:,:,1:2)/zctime |
|---|
| 1160 | Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime |
|---|
| 1161 | Fsedim(:,klev+1) = 0. |
|---|
| 1162 | |
|---|
| 1163 | else if (cl_scheme == 2) then |
|---|
| 1164 | |
|---|
| 1165 | ! sedimentation for detailed microphysics |
|---|
| 1166 | |
|---|
| 1167 | d_tr_sed(:,:,:) = 0. |
|---|
| 1168 | |
|---|
| 1169 | do i = 1, klon |
|---|
| 1170 | |
|---|
| 1171 | ! mode 1 |
|---|
| 1172 | |
|---|
| 1173 | m0_mode1(:,1) = tr_seri(i,:,i_m0_mode1drop) |
|---|
| 1174 | m0_mode1(:,2) = tr_seri(i,:,i_m0_mode1ccn) |
|---|
| 1175 | m3_mode1(:,1) = tr_seri(i,:,i_m3_mode1sa) |
|---|
| 1176 | m3_mode1(:,2) = tr_seri(i,:,i_m3_mode1w) |
|---|
| 1177 | m3_mode1(:,3) = tr_seri(i,:,i_m3_mode1ccn) |
|---|
| 1178 | |
|---|
| 1179 | call drop_sedimentation(zctime,klev,m0_mode1,m3_mode1, |
|---|
| 1180 | $ paprs(i,:),zzlev(i,:), |
|---|
| 1181 | $ zzlay(i,:),t_seri(i,:),1, |
|---|
| 1182 | $ d_ccn_sed(:,1),d_drop_sed, |
|---|
| 1183 | $ d_ccn_sed(:,2),d_liq_sed) |
|---|
| 1184 | |
|---|
| 1185 | d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop) |
|---|
| 1186 | $ + d_drop_sed |
|---|
| 1187 | d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn) |
|---|
| 1188 | $ + d_ccn_sed(:,1) |
|---|
| 1189 | d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn) |
|---|
| 1190 | $ + d_ccn_sed(:,2) |
|---|
| 1191 | d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa) |
|---|
| 1192 | $ + d_liq_sed(:,1) |
|---|
| 1193 | d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w) |
|---|
| 1194 | $ + d_liq_sed(:,2) |
|---|
| 1195 | |
|---|
| 1196 | ! mode 2 |
|---|
| 1197 | |
|---|
| 1198 | m0_mode2(:,1) = tr_seri(i,:,i_m0_mode2drop) |
|---|
| 1199 | m0_mode2(:,2) = tr_seri(i,:,i_m0_mode2ccn) |
|---|
| 1200 | m3_mode2(:,1) = tr_seri(i,:,i_m3_mode2sa) |
|---|
| 1201 | m3_mode2(:,2) = tr_seri(i,:,i_m3_mode2w) |
|---|
| 1202 | m3_mode2(:,3) = tr_seri(i,:,i_m3_mode2ccn) |
|---|
| 1203 | |
|---|
| 1204 | call drop_sedimentation(zctime,klev,m0_mode2,m3_mode2, |
|---|
| 1205 | $ paprs(i,:),zzlev(i,:), |
|---|
| 1206 | $ zzlay(i,:),t_seri(i,:),2, |
|---|
| 1207 | $ d_ccn_sed(:,1),d_drop_sed, |
|---|
| 1208 | $ d_ccn_sed(:,2),d_liq_sed) |
|---|
| 1209 | |
|---|
| 1210 | d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop) |
|---|
| 1211 | $ + d_drop_sed |
|---|
| 1212 | d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn) |
|---|
| 1213 | $ + d_ccn_sed(:,1) |
|---|
| 1214 | d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn) |
|---|
| 1215 | $ + d_ccn_sed(:,2) |
|---|
| 1216 | d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa) |
|---|
| 1217 | $ + d_liq_sed(:,1) |
|---|
| 1218 | d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w) |
|---|
| 1219 | $ + d_liq_sed(:,2) |
|---|
| 1220 | |
|---|
| 1221 | ! aer |
|---|
| 1222 | |
|---|
| 1223 | call aer_sedimentation(zctime,klev, |
|---|
| 1224 | $ tr_seri(i,:,i_m0_aer), |
|---|
| 1225 | $ tr_seri(i,:,i_m3_aer), |
|---|
| 1226 | $ paprs(i,:),zzlev(i,:), |
|---|
| 1227 | $ zzlay(i,:),t_seri(i,:), |
|---|
| 1228 | $ d_tr_sed(i,:,i_m0_aer), |
|---|
| 1229 | $ d_tr_sed(i,:,i_m3_aer), |
|---|
| 1230 | $ aer_flux) |
|---|
| 1231 | |
|---|
| 1232 | end do |
|---|
| 1233 | |
|---|
| 1234 | ! tendency due to sedimentation |
|---|
| 1235 | |
|---|
| 1236 | do iq = nqmax-nmicro+1,nqmax |
|---|
| 1237 | d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime |
|---|
| 1238 | end do |
|---|
| 1239 | #endif |
|---|
| 1240 | end if ! cl_scheme |
|---|
| 1241 | |
|---|
| 1242 | ! update gaseous tracers (chemistry) |
|---|
| 1243 | |
|---|
| 1244 | do iq = 1, nqmax - nmicro |
|---|
| 1245 | tr_seri(:,:,iq) = max(tr_seri(:,:,iq) |
|---|
| 1246 | $ + d_tr_chem(:,:,iq)*zctime,1.e-30) |
|---|
| 1247 | end do |
|---|
| 1248 | |
|---|
| 1249 | ! update condensed tracers (condensation + sedimentation) |
|---|
| 1250 | |
|---|
| 1251 | if (cl_scheme == 1) then |
|---|
| 1252 | tr_seri(:,:,i_h2so4liq) = max(tr_seri(:,:,i_h2so4liq) |
|---|
| 1253 | $ + d_tr_sed(:,:,1)*zctime, 1.e-30) |
|---|
| 1254 | tr_seri(:,:,i_h2oliq) = max(tr_seri(:,:,i_h2oliq) |
|---|
| 1255 | $ + d_tr_sed(:,:,2)*zctime, 1.e-30) |
|---|
| 1256 | else if (cl_scheme == 2) then |
|---|
| 1257 | do iq = nqmax-nmicro+1,nqmax |
|---|
| 1258 | tr_seri(:,:,iq) = max(tr_seri(:,:,iq) |
|---|
| 1259 | $ + d_tr_sed(:,:,iq)*zctime,1.e-30) |
|---|
| 1260 | end do |
|---|
| 1261 | end if ! cl_scheme |
|---|
| 1262 | |
|---|
| 1263 | end if ! ok_sedim |
|---|
| 1264 | end if ! mod(itap,chempas) <------- end of chemistry supercycling |
|---|
| 1265 | |
|---|
| 1266 | !========= |
|---|
| 1267 | ! End Case 3: Full chemistry and/or clouds. |
|---|
| 1268 | !==================================================================== |
|---|
| 1269 | |
|---|
| 1270 | end if ! tr_scheme |
|---|
| 1271 | end if ! iflag_trac |
|---|
| 1272 | |
|---|
| 1273 | c==================================================================== |
|---|
| 1274 | c Appeler la diffusion verticale (programme de couche limite) |
|---|
| 1275 | c==================================================================== |
|---|
| 1276 | |
|---|
| 1277 | c------------------------------- |
|---|
| 1278 | c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force |
|---|
| 1279 | c l'equilibre radiatif du sol |
|---|
| 1280 | if (.not. ok_clmain) then |
|---|
| 1281 | if (debut) then |
|---|
| 1282 | print*,"ATTENTION, CLMAIN SHUNTEE..." |
|---|
| 1283 | endif |
|---|
| 1284 | |
|---|
| 1285 | DO i = 1, klon |
|---|
| 1286 | sens(i) = 0.0e0 ! flux de chaleur sensible au sol |
|---|
| 1287 | fder(i) = 0.0e0 |
|---|
| 1288 | dlw(i) = 0.0e0 |
|---|
| 1289 | ENDDO |
|---|
| 1290 | |
|---|
| 1291 | c Incrementer la temperature du sol |
|---|
| 1292 | c |
|---|
| 1293 | DO i = 1, klon |
|---|
| 1294 | d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200 |
|---|
| 1295 | ftsol(i) = ftsol(i) + d_ts(i) |
|---|
| 1296 | do j=1,nsoilmx |
|---|
| 1297 | ftsoil(i,j)=ftsol(i) |
|---|
| 1298 | enddo |
|---|
| 1299 | ENDDO |
|---|
| 1300 | |
|---|
| 1301 | c------------------------------- |
|---|
| 1302 | else |
|---|
| 1303 | c------------------------------- |
|---|
| 1304 | |
|---|
| 1305 | fder = dlw |
|---|
| 1306 | |
|---|
| 1307 | ! ADAPTATION GCM POUR CP(T) |
|---|
| 1308 | |
|---|
| 1309 | if (physideal) then |
|---|
| 1310 | CALL clmain_ideal(dtime,itap, |
|---|
| 1311 | e t_seri,u_seri,v_seri, |
|---|
| 1312 | e rmu0, |
|---|
| 1313 | e ftsol, |
|---|
| 1314 | $ ftsoil, |
|---|
| 1315 | $ paprs,pplay,ppk,radsol,falbe, |
|---|
| 1316 | e solsw, sollw, sollwdown, fder, |
|---|
| 1317 | e longitude_deg, latitude_deg, dx, dy, |
|---|
| 1318 | e debut, lafin, |
|---|
| 1319 | s d_t_vdf,d_u_vdf,d_v_vdf,d_ts, |
|---|
| 1320 | s fluxt,fluxu,fluxv,cdragh,cdragm, |
|---|
| 1321 | s dsens, |
|---|
| 1322 | s ycoefh,yu1,yv1) |
|---|
| 1323 | else |
|---|
| 1324 | CALL clmain(dtime,itap, |
|---|
| 1325 | e t_seri,u_seri,v_seri, |
|---|
| 1326 | e rmu0, |
|---|
| 1327 | e ftsol, |
|---|
| 1328 | $ ftsoil, |
|---|
| 1329 | $ paprs,pplay,ppk,radsol,falbe, |
|---|
| 1330 | e solsw, sollw, sollwdown, fder, |
|---|
| 1331 | e longitude_deg, latitude_deg, dx, dy, |
|---|
| 1332 | & q2, |
|---|
| 1333 | e debut, lafin, |
|---|
| 1334 | s d_t_vdf,d_u_vdf,d_v_vdf,d_ts, |
|---|
| 1335 | s fluxt,fluxu,fluxv,cdragh,cdragm, |
|---|
| 1336 | s dsens, |
|---|
| 1337 | s ycoefh,yu1,yv1) |
|---|
| 1338 | endif |
|---|
| 1339 | |
|---|
| 1340 | CXXX Incrementation des flux |
|---|
| 1341 | DO i = 1, klon |
|---|
| 1342 | sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol |
|---|
| 1343 | fder(i) = dlw(i) + dsens(i) |
|---|
| 1344 | ENDDO |
|---|
| 1345 | CXXX |
|---|
| 1346 | IF (.not. turb_resolved) then !True only for LES |
|---|
| 1347 | DO k = 1, klev |
|---|
| 1348 | DO i = 1, klon |
|---|
| 1349 | t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k) |
|---|
| 1350 | d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s |
|---|
| 1351 | u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k) |
|---|
| 1352 | d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s |
|---|
| 1353 | v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k) |
|---|
| 1354 | d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s |
|---|
| 1355 | ENDDO |
|---|
| 1356 | ENDDO |
|---|
| 1357 | ENDIF |
|---|
| 1358 | C TRACEURS |
|---|
| 1359 | |
|---|
| 1360 | if (iflag_trac.eq.1) then |
|---|
| 1361 | DO k = 1, klev |
|---|
| 1362 | DO i = 1, klon |
|---|
| 1363 | delp(i,k) = paprs(i,k)-paprs(i,k+1) |
|---|
| 1364 | ENDDO |
|---|
| 1365 | ENDDO |
|---|
| 1366 | |
|---|
| 1367 | if (klon.eq.1) then !For the 1D model |
|---|
| 1368 | !Reading the prescribed profile from deftank |
|---|
| 1369 | open(35,file='kzz_p.txt',form='formatted',status='old') |
|---|
| 1370 | rewind(35) |
|---|
| 1371 | DO k=1,klev |
|---|
| 1372 | read(35,*) kzz_p(k) |
|---|
| 1373 | ENDDO |
|---|
| 1374 | close(35) |
|---|
| 1375 | |
|---|
| 1376 | !Implementing the new profile in m2/s |
|---|
| 1377 | DO k = 1, klev |
|---|
| 1378 | ycoefh(1,k) = kzz_p(k) |
|---|
| 1379 | ENDDO |
|---|
| 1380 | endif |
|---|
| 1381 | |
|---|
| 1382 | DO iq=1, nqmax |
|---|
| 1383 | |
|---|
| 1384 | CALL cltrac(dtime,ycoefh,t_seri, |
|---|
| 1385 | s tr_seri(:,:,iq),source(:,iq), |
|---|
| 1386 | e paprs, pplay,delp, |
|---|
| 1387 | s d_tr_vdf(:,:,iq)) |
|---|
| 1388 | |
|---|
| 1389 | tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq) |
|---|
| 1390 | d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s |
|---|
| 1391 | |
|---|
| 1392 | ENDDO !nqmax |
|---|
| 1393 | |
|---|
| 1394 | endif |
|---|
| 1395 | |
|---|
| 1396 | IF (if_ebil.ge.2) THEN |
|---|
| 1397 | ztit='after clmain' |
|---|
| 1398 | CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime |
|---|
| 1399 | e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay |
|---|
| 1400 | s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) |
|---|
| 1401 | call diagphy(cell_area,ztit,ip_ebil |
|---|
| 1402 | e , zero_v, zero_v, zero_v, zero_v, sens |
|---|
| 1403 | e , zero_v, zero_v, zero_v, ztsol |
|---|
| 1404 | e , d_h_vcol, d_qt, d_ec |
|---|
| 1405 | s , fs_bound, fq_bound ) |
|---|
| 1406 | END IF |
|---|
| 1407 | C |
|---|
| 1408 | c |
|---|
| 1409 | c Incrementer la temperature du sol |
|---|
| 1410 | c |
|---|
| 1411 | DO i = 1, klon |
|---|
| 1412 | ftsol(i) = ftsol(i) + d_ts(i) |
|---|
| 1413 | ENDDO |
|---|
| 1414 | |
|---|
| 1415 | c Calculer la derive du flux infrarouge |
|---|
| 1416 | c |
|---|
| 1417 | DO i = 1, klon |
|---|
| 1418 | dlw(i) = - 4.0*RSIGMA*ftsol(i)**3 |
|---|
| 1419 | ENDDO |
|---|
| 1420 | |
|---|
| 1421 | c------------------------------- |
|---|
| 1422 | endif ! fin du VENUS TEST |
|---|
| 1423 | |
|---|
| 1424 | ! tests: output tendencies |
|---|
| 1425 | ! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev) |
|---|
| 1426 | ! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev) |
|---|
| 1427 | ! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev) |
|---|
| 1428 | ! call writefield_phy('physiq_d_ts',d_ts,1) |
|---|
| 1429 | |
|---|
| 1430 | c=================================================================== |
|---|
| 1431 | c Convection seche |
|---|
| 1432 | c=================================================================== |
|---|
| 1433 | c |
|---|
| 1434 | d_t_ajs(:,:)=0. |
|---|
| 1435 | d_u_ajs(:,:)=0. |
|---|
| 1436 | d_v_ajs(:,:)=0. |
|---|
| 1437 | d_tr_ajs(:,:,:)=0. |
|---|
| 1438 | c |
|---|
| 1439 | IF(prt_level>9)WRITE(lunout,*) |
|---|
| 1440 | . 'AVANT LA CONVECTION SECHE , iflag_ajs=' |
|---|
| 1441 | s ,iflag_ajs |
|---|
| 1442 | |
|---|
| 1443 | if(iflag_ajs.eq.0) then |
|---|
| 1444 | c Rien |
|---|
| 1445 | c ==== |
|---|
| 1446 | IF(prt_level>9)WRITE(lunout,*)'pas de convection' |
|---|
| 1447 | |
|---|
| 1448 | else if(iflag_ajs.eq.1) then |
|---|
| 1449 | |
|---|
| 1450 | c Ajustement sec |
|---|
| 1451 | c ============== |
|---|
| 1452 | IF(prt_level>9)WRITE(lunout,*)'ajsec' |
|---|
| 1453 | |
|---|
| 1454 | ! ADAPTATION GCM POUR CP(T) |
|---|
| 1455 | CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax, |
|---|
| 1456 | . tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs) |
|---|
| 1457 | |
|---|
| 1458 | ! ADAPTATION GCM POUR CP(T) |
|---|
| 1459 | do i=1,klon |
|---|
| 1460 | flux_ajs(i,1) = 0.0 |
|---|
| 1461 | do j=2,klev |
|---|
| 1462 | flux_ajs(i,j) = flux_ajs(i,j-1) |
|---|
| 1463 | . + cpnew(i,j-1)/RG*d_t_ajs(i,j-1)/dtime |
|---|
| 1464 | . *(paprs(i,j-1)-paprs(i,j)) |
|---|
| 1465 | enddo |
|---|
| 1466 | enddo |
|---|
| 1467 | |
|---|
| 1468 | t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:) |
|---|
| 1469 | d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s |
|---|
| 1470 | u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:) |
|---|
| 1471 | d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s |
|---|
| 1472 | v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:) |
|---|
| 1473 | d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s |
|---|
| 1474 | |
|---|
| 1475 | if (iflag_trac.eq.1) then |
|---|
| 1476 | tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:) |
|---|
| 1477 | d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s |
|---|
| 1478 | endif |
|---|
| 1479 | endif |
|---|
| 1480 | |
|---|
| 1481 | ! tests: output tendencies |
|---|
| 1482 | ! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev) |
|---|
| 1483 | ! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev) |
|---|
| 1484 | ! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev) |
|---|
| 1485 | c |
|---|
| 1486 | IF (if_ebil.ge.2) THEN |
|---|
| 1487 | ztit='after dry_adjust' |
|---|
| 1488 | CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime |
|---|
| 1489 | e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay |
|---|
| 1490 | s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) |
|---|
| 1491 | call diagphy(cell_area,ztit,ip_ebil |
|---|
| 1492 | e , zero_v, zero_v, zero_v, zero_v, sens |
|---|
| 1493 | e , zero_v, zero_v, zero_v, ztsol |
|---|
| 1494 | e , d_h_vcol, d_qt, d_ec |
|---|
| 1495 | s , fs_bound, fq_bound ) |
|---|
| 1496 | END IF |
|---|
| 1497 | |
|---|
| 1498 | c==================================================================== |
|---|
| 1499 | c RAYONNEMENT |
|---|
| 1500 | c==================================================================== |
|---|
| 1501 | if (mod(itap,radpas) == 0) then |
|---|
| 1502 | |
|---|
| 1503 | dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s) |
|---|
| 1504 | |
|---|
| 1505 | ! update mmean |
|---|
| 1506 | if (ok_chem) then |
|---|
| 1507 | mmean(:,:) = 0. |
|---|
| 1508 | do iq = 1,nqmax - nmicro |
|---|
| 1509 | mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq) |
|---|
| 1510 | enddo |
|---|
| 1511 | mmean(:,:) = 1./mmean(:,:) |
|---|
| 1512 | rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K |
|---|
| 1513 | endif |
|---|
| 1514 | |
|---|
| 1515 | cc--------------------------------------------- |
|---|
| 1516 | if (callnlte .or. callthermos) then |
|---|
| 1517 | if (ok_chem) then |
|---|
| 1518 | |
|---|
| 1519 | ! nlte : use computed chemical species |
|---|
| 1520 | |
|---|
| 1521 | co2vmr_gcm(:,:) = tr_seri(:,:,i_co2)*mmean(:,:)/m_tr(i_co2) |
|---|
| 1522 | covmr_gcm(:,:) = tr_seri(:,:,i_co)*mmean(:,:)/m_tr(i_co) |
|---|
| 1523 | ovmr_gcm(:,:) = tr_seri(:,:,i_o)*mmean(:,:)/m_tr(i_o) |
|---|
| 1524 | n2vmr_gcm(:,:) = tr_seri(:,:,i_n2)*mmean(:,:)/m_tr(i_n2) |
|---|
| 1525 | |
|---|
| 1526 | else |
|---|
| 1527 | |
|---|
| 1528 | ! nlte : use hedin climatology |
|---|
| 1529 | |
|---|
| 1530 | call compo_hedin83_mod(pplay,rmu0, |
|---|
| 1531 | $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) |
|---|
| 1532 | end if |
|---|
| 1533 | end if |
|---|
| 1534 | |
|---|
| 1535 | c NLTE cooling from CO2 emission |
|---|
| 1536 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1537 | |
|---|
| 1538 | IF(callnlte) THEN |
|---|
| 1539 | if(nltemodel.eq.0.or.nltemodel.eq.1) then |
|---|
| 1540 | c nltecool call not correct... |
|---|
| 1541 | c CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri, |
|---|
| 1542 | c $ tr_seri, d_t_nlte) |
|---|
| 1543 | abort_message='nltemodel=0 or 1 should not be used...' |
|---|
| 1544 | call abort_physic(modname,abort_message,1) |
|---|
| 1545 | else if(nltemodel.eq.2) then |
|---|
| 1546 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1547 | !! HEDIN instead of compo for this calculation |
|---|
| 1548 | ! call compo_hedin83_mod(pplay,rmu0, |
|---|
| 1549 | ! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) |
|---|
| 1550 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1551 | CALL nlte_tcool(klon,klev,pplay*9.869e-6, |
|---|
| 1552 | $ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, |
|---|
| 1553 | $ ovmr_gcm,d_t_nlte,ierr_nlte,varerr ) |
|---|
| 1554 | if(ierr_nlte.gt.0) then |
|---|
| 1555 | write(*,*) |
|---|
| 1556 | $ 'WARNING: nlte_tcool output with error message', |
|---|
| 1557 | $ 'ierr_nlte=',ierr_nlte,'varerr=',varerr |
|---|
| 1558 | write(*,*)'I will continue anyway' |
|---|
| 1559 | endif |
|---|
| 1560 | endif |
|---|
| 1561 | |
|---|
| 1562 | ELSE |
|---|
| 1563 | |
|---|
| 1564 | d_t_nlte(:,:)=0. |
|---|
| 1565 | |
|---|
| 1566 | ENDIF |
|---|
| 1567 | |
|---|
| 1568 | c Find number of layers for LTE radiation calculations |
|---|
| 1569 | |
|---|
| 1570 | IF(callnlte .or. callnirco2) |
|---|
| 1571 | $ CALL nlthermeq(klon, klev, paprs, pplay) |
|---|
| 1572 | |
|---|
| 1573 | cc--------------------------------------------- |
|---|
| 1574 | c LTE radiative transfert / solar / IR matrix |
|---|
| 1575 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1576 | if (physideal) then |
|---|
| 1577 | CALL radlwsw_newtoncool(presnivs,t_seri) |
|---|
| 1578 | else |
|---|
| 1579 | CALL radlwsw |
|---|
| 1580 | e (dist, rmu0, fract, zzlev, |
|---|
| 1581 | e paprs, pplay,ftsol, t_seri) |
|---|
| 1582 | endif |
|---|
| 1583 | |
|---|
| 1584 | c ALBEDO VARIATIONS: test for Yeon Joo Lee |
|---|
| 1585 | c increment to increase it for 20 Vd => +80% |
|---|
| 1586 | c heat(:,:)=heat(:,:)*(1.+0.80*((rjourvrai-356)+gmtime)/20.) |
|---|
| 1587 | c or to decrease it for 20 Vd => 1/1.8 |
|---|
| 1588 | c heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.) |
|---|
| 1589 | |
|---|
| 1590 | c ------------ ALBEDO VARIATIONS: scenarios for VCD |
|---|
| 1591 | c shape of relative variation from Lee et al 2019 (Fig 13b) |
|---|
| 1592 | c between 57 km (4e4 Pa) and 72 km (2.5e3 Pa), peak at 67 km (6e3 Pa) |
|---|
| 1593 | c do j=1,klev |
|---|
| 1594 | c factAlb = 0. |
|---|
| 1595 | c if ((presnivs(j).gt.6.e3).and.(presnivs(j).lt.4.e4)) then |
|---|
| 1596 | c factAlb = (log(presnivs(j))-log(4.e4))/(log(6.e3)-log(4.e4)) |
|---|
| 1597 | c elseif ((presnivs(j).lt.6.e3).and.(presnivs(j).gt.2.5e3)) then |
|---|
| 1598 | c factAlb = (log(presnivs(j))-log(2.5e3))/(log(6.e3)-log(2.5e3)) |
|---|
| 1599 | c endif |
|---|
| 1600 | c Increase by 50% (Minimum albedo) |
|---|
| 1601 | c heat(:,j)=heat(:,j)*(1+factAlb*0.5) |
|---|
| 1602 | c Decrease by 30% (Maximum albedo) |
|---|
| 1603 | c heat(:,j)=heat(:,j)*(1-factAlb*0.3) |
|---|
| 1604 | c enddo |
|---|
| 1605 | c ------------ END ALBEDO VARIATIONS |
|---|
| 1606 | |
|---|
| 1607 | cc--------------------------------------------- |
|---|
| 1608 | c CO2 near infrared absorption |
|---|
| 1609 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1610 | |
|---|
| 1611 | d_t_nirco2(:,:)=0. |
|---|
| 1612 | if (callnirco2) then |
|---|
| 1613 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1614 | !! HEDIN instead of compo for this calculation |
|---|
| 1615 | ! call compo_hedin83_mod(pplay,rmu0, |
|---|
| 1616 | ! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) |
|---|
| 1617 | ! tr_hedin(:,:,:)=tr_seri(:,:,:) |
|---|
| 1618 | ! tr_hedin(:,:,i_co2)=co2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_co2) |
|---|
| 1619 | ! tr_hedin(:,:,i_co) = covmr_gcm(:,:)/mmean(:,:)*m_tr(i_co) |
|---|
| 1620 | ! tr_hedin(:,:,i_o) = ovmr_gcm(:,:)/mmean(:,:)*m_tr(i_o) |
|---|
| 1621 | ! tr_hedin(:,:,i_n2) = n2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_n2) |
|---|
| 1622 | ! call nirco2abs (klon, klev, pplay, dist, nqmax, tr_hedin, |
|---|
| 1623 | ! . rmu0, fract, d_t_nirco2) |
|---|
| 1624 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1625 | call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri, |
|---|
| 1626 | . rmu0, fract, d_t_nirco2) |
|---|
| 1627 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1628 | endif |
|---|
| 1629 | |
|---|
| 1630 | |
|---|
| 1631 | cc--------------------------------------------- |
|---|
| 1632 | c Net atmospheric radiative heating rate (K.s-1) |
|---|
| 1633 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1634 | |
|---|
| 1635 | IF(callnlte.or.callnirco2) THEN |
|---|
| 1636 | CALL blendrad(klon, klev, pplay,heat, |
|---|
| 1637 | & cool, d_t_nirco2,d_t_nlte, dtsw, dtlw) |
|---|
| 1638 | ELSE |
|---|
| 1639 | dtsw(:,:)=heat(:,:) |
|---|
| 1640 | dtlw(:,:)=-1*cool(:,:) |
|---|
| 1641 | ENDIF |
|---|
| 1642 | |
|---|
| 1643 | DO k=1,klev |
|---|
| 1644 | DO i=1,klon |
|---|
| 1645 | d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s |
|---|
| 1646 | ENDDO |
|---|
| 1647 | ENDDO |
|---|
| 1648 | |
|---|
| 1649 | |
|---|
| 1650 | cc--------------------------------------------- |
|---|
| 1651 | c EUV heating rate (K.s-1) |
|---|
| 1652 | c ~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1653 | |
|---|
| 1654 | d_t_euv(:,:)=0. |
|---|
| 1655 | |
|---|
| 1656 | IF (callthermos) THEN |
|---|
| 1657 | |
|---|
| 1658 | call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay, |
|---|
| 1659 | $ rmu0,dtimerad,gmtime, |
|---|
| 1660 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1661 | !! HEDIN instead of compo for this calculation |
|---|
| 1662 | !! cf nlte_tcool and nirco2abs above !! |
|---|
| 1663 | ! $ tr_hedin, d_tr, d_t_euv ) |
|---|
| 1664 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1665 | $ tr_seri, d_tr, d_t_euv ) |
|---|
| 1666 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 1667 | |
|---|
| 1668 | DO k=1,klev |
|---|
| 1669 | DO ig=1,klon |
|---|
| 1670 | d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k) |
|---|
| 1671 | ENDDO |
|---|
| 1672 | ENDDO |
|---|
| 1673 | |
|---|
| 1674 | ENDIF ! callthermos |
|---|
| 1675 | |
|---|
| 1676 | ENDIF ! radpas |
|---|
| 1677 | c==================================================================== |
|---|
| 1678 | c |
|---|
| 1679 | c Ajouter la tendance des rayonnements (tous les pas) |
|---|
| 1680 | c |
|---|
| 1681 | DO k = 1, klev |
|---|
| 1682 | DO i = 1, klon |
|---|
| 1683 | t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime |
|---|
| 1684 | ENDDO |
|---|
| 1685 | ENDDO |
|---|
| 1686 | |
|---|
| 1687 | ! increment physics counter |
|---|
| 1688 | |
|---|
| 1689 | itap = itap + 1 |
|---|
| 1690 | c==================================================================== |
|---|
| 1691 | |
|---|
| 1692 | |
|---|
| 1693 | c============================== |
|---|
| 1694 | ! -- MOLECULAR DIFFUSION --- |
|---|
| 1695 | c============================== |
|---|
| 1696 | |
|---|
| 1697 | d_q_moldif(:,:,:)=0 |
|---|
| 1698 | |
|---|
| 1699 | IF (callthermos .and. ok_chem) THEN |
|---|
| 1700 | |
|---|
| 1701 | call moldiff_red(klon, klev, nqmax, |
|---|
| 1702 | & pplay,paprs,t_seri, tr_seri, pdtphys, |
|---|
| 1703 | & d_t_euv,d_t_conduc,d_q_moldif) |
|---|
| 1704 | |
|---|
| 1705 | |
|---|
| 1706 | ! --- update tendencies tracers --- |
|---|
| 1707 | |
|---|
| 1708 | DO iq = 1, nqmax |
|---|
| 1709 | DO k=1,klev |
|---|
| 1710 | DO ig=1,klon |
|---|
| 1711 | tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+ |
|---|
| 1712 | & d_q_moldif(ig,k,iq)*dtime,1.e-30) |
|---|
| 1713 | ENDDO |
|---|
| 1714 | ENDDO |
|---|
| 1715 | ENDDO |
|---|
| 1716 | |
|---|
| 1717 | ENDIF ! callthermos & ok_chem |
|---|
| 1718 | |
|---|
| 1719 | c==================================================================== |
|---|
| 1720 | |
|---|
| 1721 | c================================== |
|---|
| 1722 | ! -- ION AMBIPOLAR DIFFUSION --- |
|---|
| 1723 | c================================== |
|---|
| 1724 | |
|---|
| 1725 | d_q_iondif(:,:,:)=0 |
|---|
| 1726 | |
|---|
| 1727 | IF (callthermos .and. ok_chem .and. ok_ionchem) THEN |
|---|
| 1728 | IF (ok_iondiff) THEN |
|---|
| 1729 | |
|---|
| 1730 | call iondiff_red(klon,klev,nqmax,pplay,paprs, |
|---|
| 1731 | & t_seri,tr_seri,pphis, |
|---|
| 1732 | & gmtime,latitude_deg,longitude_deg, |
|---|
| 1733 | & pdtphys,d_t_euv,d_t_conduc,d_q_moldif, |
|---|
| 1734 | & d_q_iondif) |
|---|
| 1735 | |
|---|
| 1736 | !write (*,*) 'TITI EST PASSE PAR LA' |
|---|
| 1737 | ! --- update tendencies tracers --- |
|---|
| 1738 | |
|---|
| 1739 | DO iq = 1, nqmax |
|---|
| 1740 | IF (type_tr(iq) .eq. 2) THEN |
|---|
| 1741 | DO k=1,klev |
|---|
| 1742 | DO ig=1,klon |
|---|
| 1743 | tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+ |
|---|
| 1744 | & d_q_iondif(ig,k,iq)*dtime,1.e-30) |
|---|
| 1745 | ENDDO |
|---|
| 1746 | ENDDO |
|---|
| 1747 | ENDIF |
|---|
| 1748 | ENDDO |
|---|
| 1749 | ENDIF ! ok_iondiff |
|---|
| 1750 | ENDIF ! callthermos & ok_chem & ok_ionchem |
|---|
| 1751 | |
|---|
| 1752 | c==================================================================== |
|---|
| 1753 | ! tests: output tendencies |
|---|
| 1754 | ! call writefield_phy('physiq_dtrad',dtrad,klev) |
|---|
| 1755 | |
|---|
| 1756 | IF (if_ebil.ge.2) THEN |
|---|
| 1757 | ztit='after rad' |
|---|
| 1758 | CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime |
|---|
| 1759 | e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay |
|---|
| 1760 | s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) |
|---|
| 1761 | call diagphy(cell_area,ztit,ip_ebil |
|---|
| 1762 | e , topsw, toplw, solsw, sollw, zero_v |
|---|
| 1763 | e , zero_v, zero_v, zero_v, ztsol |
|---|
| 1764 | e , d_h_vcol, d_qt, d_ec |
|---|
| 1765 | s , fs_bound, fq_bound ) |
|---|
| 1766 | END IF |
|---|
| 1767 | c |
|---|
| 1768 | |
|---|
| 1769 | c==================================================================== |
|---|
| 1770 | c Calcul des gravity waves FLOTT |
|---|
| 1771 | c==================================================================== |
|---|
| 1772 | c |
|---|
| 1773 | c if (ok_orodr.or.ok_gw_nonoro) then |
|---|
| 1774 | |
|---|
| 1775 | c CALCUL DE N2 |
|---|
| 1776 | c UTILISE LA RELATION ENTRE N2 ET STABILITE |
|---|
| 1777 | c N2 = RG/T (dT/dz+RG/cp(T)) |
|---|
| 1778 | c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta. |
|---|
| 1779 | |
|---|
| 1780 | do i=1,klon |
|---|
| 1781 | do k=2,klev |
|---|
| 1782 | ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2. |
|---|
| 1783 | enddo |
|---|
| 1784 | enddo |
|---|
| 1785 | do i=1,klon |
|---|
| 1786 | do k=2,klev |
|---|
| 1787 | ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2. |
|---|
| 1788 | zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1) |
|---|
| 1789 | zdzlev(i,k) = (zzlay(i,k)-zzlay(i,k-1)) |
|---|
| 1790 | zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k) |
|---|
| 1791 | . + RG/cpnew(i,k) ) |
|---|
| 1792 | zn2(i,k) = max(zn2(i,k),1.e-12) ! securite |
|---|
| 1793 | enddo |
|---|
| 1794 | zn2(i,1) = 1.e-12 ! securite |
|---|
| 1795 | enddo |
|---|
| 1796 | |
|---|
| 1797 | c endif |
|---|
| 1798 | |
|---|
| 1799 | c ----------------------------ORODRAG |
|---|
| 1800 | IF (ok_orodr) THEN |
|---|
| 1801 | c |
|---|
| 1802 | c selection des points pour lesquels le shema est actif: |
|---|
| 1803 | igwd=0 |
|---|
| 1804 | DO i=1,klon |
|---|
| 1805 | itest(i)=0 |
|---|
| 1806 | c IF ((zstd(i).gt.10.0)) THEN |
|---|
| 1807 | IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN |
|---|
| 1808 | itest(i)=1 |
|---|
| 1809 | igwd=igwd+1 |
|---|
| 1810 | idx(igwd)=i |
|---|
| 1811 | ENDIF |
|---|
| 1812 | ENDDO |
|---|
| 1813 | c igwdim=MAX(1,igwd) |
|---|
| 1814 | c |
|---|
| 1815 | c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ] |
|---|
| 1816 | CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2, |
|---|
| 1817 | e zmea,zstd, zsig, zgam, zthe,zpic,zval, |
|---|
| 1818 | e igwd,idx,itest, |
|---|
| 1819 | e t_seri, u_seri, v_seri, |
|---|
| 1820 | s zulow, zvlow, zustrdr, zvstrdr, |
|---|
| 1821 | s d_t_oro, d_u_oro, d_v_oro, |
|---|
| 1822 | s zublstrdr,zvblstrdr,znlow,zeff,zbl, |
|---|
| 1823 | s ztau,tau0,knu2,kbreak) |
|---|
| 1824 | |
|---|
| 1825 | c print*,"d_u_oro=",d_u_oro(klon/2,:) |
|---|
| 1826 | c ajout des tendances |
|---|
| 1827 | t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:) |
|---|
| 1828 | d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s |
|---|
| 1829 | u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:) |
|---|
| 1830 | d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s |
|---|
| 1831 | v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:) |
|---|
| 1832 | d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s |
|---|
| 1833 | c |
|---|
| 1834 | ELSE |
|---|
| 1835 | d_t_oro = 0. |
|---|
| 1836 | d_u_oro = 0. |
|---|
| 1837 | d_v_oro = 0. |
|---|
| 1838 | zustrdr = 0. |
|---|
| 1839 | zvstrdr = 0. |
|---|
| 1840 | zublstrdr = 0. |
|---|
| 1841 | zvblstrdr = 0. |
|---|
| 1842 | znlow = 0. |
|---|
| 1843 | zeff = 0. |
|---|
| 1844 | zbl = 0 |
|---|
| 1845 | knu2 = 0 |
|---|
| 1846 | kbreak = 0 |
|---|
| 1847 | ztau = 0 |
|---|
| 1848 | tau0 = 0. |
|---|
| 1849 | c |
|---|
| 1850 | ENDIF ! fin de test sur ok_orodr |
|---|
| 1851 | c |
|---|
| 1852 | ! tests: output tendencies |
|---|
| 1853 | ! call writefield_phy('physiq_d_t_oro',d_t_oro,klev) |
|---|
| 1854 | ! call writefield_phy('physiq_d_u_oro',d_u_oro,klev) |
|---|
| 1855 | ! call writefield_phy('physiq_d_v_oro',d_v_oro,klev) |
|---|
| 1856 | |
|---|
| 1857 | c ----------------------------OROLIFT |
|---|
| 1858 | IF (ok_orolf) THEN |
|---|
| 1859 | print*,"ok_orolf NOT IMPLEMENTED !" |
|---|
| 1860 | stop |
|---|
| 1861 | c |
|---|
| 1862 | c selection des points pour lesquels le shema est actif: |
|---|
| 1863 | igwd=0 |
|---|
| 1864 | DO i=1,klon |
|---|
| 1865 | itest(i)=0 |
|---|
| 1866 | IF ((zpic(i)-zmea(i)).GT.100.) THEN |
|---|
| 1867 | itest(i)=1 |
|---|
| 1868 | igwd=igwd+1 |
|---|
| 1869 | idx(igwd)=i |
|---|
| 1870 | ENDIF |
|---|
| 1871 | ENDDO |
|---|
| 1872 | c igwdim=MAX(1,igwd) |
|---|
| 1873 | c |
|---|
| 1874 | c A ADAPTER POUR VENUS!!! |
|---|
| 1875 | c CALL lift_noro(klon,klev,dtime,paprs,pplay, |
|---|
| 1876 | c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, |
|---|
| 1877 | c e igwd,idx,itest, |
|---|
| 1878 | c e t_seri, u_seri, v_seri, |
|---|
| 1879 | c s zulow, zvlow, zustrli, zvstrli, |
|---|
| 1880 | c s d_t_lif, d_u_lif, d_v_lif ) |
|---|
| 1881 | |
|---|
| 1882 | c |
|---|
| 1883 | c ajout des tendances |
|---|
| 1884 | t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:) |
|---|
| 1885 | d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s |
|---|
| 1886 | u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:) |
|---|
| 1887 | d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s |
|---|
| 1888 | v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:) |
|---|
| 1889 | d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s |
|---|
| 1890 | c |
|---|
| 1891 | ELSE |
|---|
| 1892 | d_t_lif = 0. |
|---|
| 1893 | d_u_lif = 0. |
|---|
| 1894 | d_v_lif = 0. |
|---|
| 1895 | zustrli = 0. |
|---|
| 1896 | zvstrli = 0. |
|---|
| 1897 | c |
|---|
| 1898 | ENDIF ! fin de test sur ok_orolf |
|---|
| 1899 | |
|---|
| 1900 | c ---------------------------- NON-ORO GRAVITY WAVES |
|---|
| 1901 | IF(ok_gw_nonoro) then |
|---|
| 1902 | |
|---|
| 1903 | ! Obsolete |
|---|
| 1904 | ! but used for VCD 1.1 |
|---|
| 1905 | ! call flott_gwd_ran(klon,klev,dtime,pplay,zn2, |
|---|
| 1906 | ! e t_seri, u_seri, v_seri, paprs(klon/2+1,:), |
|---|
| 1907 | ! o zustrhi,zvstrhi, |
|---|
| 1908 | ! o d_t_hin, d_u_hin, d_v_hin) |
|---|
| 1909 | |
|---|
| 1910 | ! New routine based on Generic |
|---|
| 1911 | ! used after VCD 1.1, for VCD 2.0 |
|---|
| 1912 | call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs, |
|---|
| 1913 | e t_seri, u_seri, v_seri, |
|---|
| 1914 | o zustrhi,zvstrhi, |
|---|
| 1915 | o d_t_hin, d_u_hin, d_v_hin) |
|---|
| 1916 | |
|---|
| 1917 | c ajout des tendances |
|---|
| 1918 | |
|---|
| 1919 | t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:) |
|---|
| 1920 | d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s |
|---|
| 1921 | u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:) |
|---|
| 1922 | d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s |
|---|
| 1923 | v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:) |
|---|
| 1924 | d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s |
|---|
| 1925 | |
|---|
| 1926 | ELSE |
|---|
| 1927 | d_t_hin = 0. |
|---|
| 1928 | d_u_hin = 0. |
|---|
| 1929 | d_v_hin = 0. |
|---|
| 1930 | zustrhi = 0. |
|---|
| 1931 | zvstrhi = 0. |
|---|
| 1932 | |
|---|
| 1933 | ENDIF ! fin de test sur ok_gw_nonoro |
|---|
| 1934 | |
|---|
| 1935 | ! tests: output tendencies |
|---|
| 1936 | ! call writefield_phy('physiq_d_t_hin',d_t_hin,klev) |
|---|
| 1937 | ! call writefield_phy('physiq_d_u_hin',d_u_hin,klev) |
|---|
| 1938 | ! call writefield_phy('physiq_d_v_hin',d_v_hin,klev) |
|---|
| 1939 | |
|---|
| 1940 | c==================================================================== |
|---|
| 1941 | c Transport de ballons |
|---|
| 1942 | c==================================================================== |
|---|
| 1943 | if (ballons.eq.1) then |
|---|
| 1944 | CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY, |
|---|
| 1945 | & latitude_deg,longitude_deg, |
|---|
| 1946 | c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM) |
|---|
| 1947 | C t,pplay,u,v,zphi) ! alt above planet average radius |
|---|
| 1948 | endif !ballons |
|---|
| 1949 | |
|---|
| 1950 | c==================================================================== |
|---|
| 1951 | c Bilan de mmt angulaire |
|---|
| 1952 | c==================================================================== |
|---|
| 1953 | if (bilansmc.eq.1) then |
|---|
| 1954 | CMODDEB FLOTT |
|---|
| 1955 | C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE) |
|---|
| 1956 | C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE |
|---|
| 1957 | |
|---|
| 1958 | DO i = 1, klon |
|---|
| 1959 | zustrph(i)=0. |
|---|
| 1960 | zvstrph(i)=0. |
|---|
| 1961 | zustrcl(i)=0. |
|---|
| 1962 | zvstrcl(i)=0. |
|---|
| 1963 | ENDDO |
|---|
| 1964 | DO k = 1, klev |
|---|
| 1965 | DO i = 1, klon |
|---|
| 1966 | zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* |
|---|
| 1967 | c (paprs(i,k)-paprs(i,k+1))/rg |
|---|
| 1968 | zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* |
|---|
| 1969 | c (paprs(i,k)-paprs(i,k+1))/rg |
|---|
| 1970 | zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)* |
|---|
| 1971 | c (paprs(i,k)-paprs(i,k+1))/rg |
|---|
| 1972 | zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)* |
|---|
| 1973 | c (paprs(i,k)-paprs(i,k+1))/rg |
|---|
| 1974 | ENDDO |
|---|
| 1975 | ENDDO |
|---|
| 1976 | |
|---|
| 1977 | CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY, |
|---|
| 1978 | C ra,rg,romega, |
|---|
| 1979 | C latitude_deg,longitude_deg,pphis, |
|---|
| 1980 | C zustrdr,zustrli,zustrcl, |
|---|
| 1981 | C zvstrdr,zvstrli,zvstrcl, |
|---|
| 1982 | C paprs,u,v) |
|---|
| 1983 | |
|---|
| 1984 | CCMODFIN FLOTT |
|---|
| 1985 | endif !bilansmc |
|---|
| 1986 | |
|---|
| 1987 | c==================================================================== |
|---|
| 1988 | c==================================================================== |
|---|
| 1989 | c Calculer le transport de l'eau et de l'energie (diagnostique) |
|---|
| 1990 | c |
|---|
| 1991 | c A REVOIR POUR VENUS... |
|---|
| 1992 | c |
|---|
| 1993 | c CALL transp (paprs,ftsol, |
|---|
| 1994 | c e t_seri, q_seri, u_seri, v_seri, zphi, |
|---|
| 1995 | c s ve, vq, ue, uq) |
|---|
| 1996 | c |
|---|
| 1997 | c==================================================================== |
|---|
| 1998 | c+jld ec_conser |
|---|
| 1999 | DO k = 1, klev |
|---|
| 2000 | DO i = 1, klon |
|---|
| 2001 | d_t_ec(i,k)=0.5/cpnew(i,k) |
|---|
| 2002 | $ *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2) |
|---|
| 2003 | t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k) |
|---|
| 2004 | d_t_ec(i,k) = d_t_ec(i,k)/dtime |
|---|
| 2005 | END DO |
|---|
| 2006 | END DO |
|---|
| 2007 | do i=1,klon |
|---|
| 2008 | flux_ec(i,1) = 0.0 |
|---|
| 2009 | do j=2,klev |
|---|
| 2010 | flux_ec(i,j) = flux_ec(i,j-1) |
|---|
| 2011 | . +cpnew(i,j-1)/RG*d_t_ec(i,j-1)*(paprs(i,j-1)-paprs(i,j)) |
|---|
| 2012 | enddo |
|---|
| 2013 | enddo |
|---|
| 2014 | |
|---|
| 2015 | c-jld ec_conser |
|---|
| 2016 | c==================================================================== |
|---|
| 2017 | IF (if_ebil.ge.1) THEN |
|---|
| 2018 | ztit='after physic' |
|---|
| 2019 | CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime |
|---|
| 2020 | e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay |
|---|
| 2021 | s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) |
|---|
| 2022 | C Comme les tendances de la physique sont ajoute dans la dynamique, |
|---|
| 2023 | C on devrait avoir que la variation d'entalpie par la dynamique |
|---|
| 2024 | C est egale a la variation de la physique au pas de temps precedent. |
|---|
| 2025 | C Donc la somme de ces 2 variations devrait etre nulle. |
|---|
| 2026 | call diagphy(cell_area,ztit,ip_ebil |
|---|
| 2027 | e , topsw, toplw, solsw, sollw, sens |
|---|
| 2028 | e , zero_v, zero_v, zero_v, ztsol |
|---|
| 2029 | e , d_h_vcol, d_qt, d_ec |
|---|
| 2030 | s , fs_bound, fq_bound ) |
|---|
| 2031 | C |
|---|
| 2032 | d_h_vcol_phy=d_h_vcol |
|---|
| 2033 | C |
|---|
| 2034 | END IF |
|---|
| 2035 | C |
|---|
| 2036 | c======================================================================= |
|---|
| 2037 | c SORTIES |
|---|
| 2038 | c======================================================================= |
|---|
| 2039 | |
|---|
| 2040 | c Convertir les incrementations en tendances |
|---|
| 2041 | c |
|---|
| 2042 | DO k = 1, klev |
|---|
| 2043 | DO i = 1, klon |
|---|
| 2044 | d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime |
|---|
| 2045 | d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime |
|---|
| 2046 | d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime |
|---|
| 2047 | ENDDO |
|---|
| 2048 | ENDDO |
|---|
| 2049 | c |
|---|
| 2050 | DO iq = 1, nqmax |
|---|
| 2051 | DO k = 1, klev |
|---|
| 2052 | DO i = 1, klon |
|---|
| 2053 | d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime |
|---|
| 2054 | ENDDO |
|---|
| 2055 | ENDDO |
|---|
| 2056 | ENDDO |
|---|
| 2057 | |
|---|
| 2058 | c mise à jour rho,mmean pour sorties |
|---|
| 2059 | if(callthermos) then |
|---|
| 2060 | call concentrations2(pplay,t_seri,tr_seri, nqmax) |
|---|
| 2061 | endif |
|---|
| 2062 | |
|---|
| 2063 | c calcul vitesse verticale en m/s |
|---|
| 2064 | DO k = 1, klev |
|---|
| 2065 | DO i = 1, klon |
|---|
| 2066 | vertwind(i,k) = -omega(i,k)/(rho(i,k)*RG) |
|---|
| 2067 | END DO |
|---|
| 2068 | END DO |
|---|
| 2069 | |
|---|
| 2070 | c------------------------ |
|---|
| 2071 | c Calcul moment cinetique |
|---|
| 2072 | c------------------------ |
|---|
| 2073 | c TEST VENUS... |
|---|
| 2074 | c mangtot = 0.0 |
|---|
| 2075 | c DO k = 1, klev |
|---|
| 2076 | c DO i = 1, klon |
|---|
| 2077 | c mang(i,k) = RA*cos(latitude(i)) |
|---|
| 2078 | c . *(u_seri(i,k)+RA*cos(latitude(i))*ROMEGA) |
|---|
| 2079 | c . *cell_area(i)*(paprs(i,k)-paprs(i,k+1))/RG |
|---|
| 2080 | c mangtot=mangtot+mang(i,k) |
|---|
| 2081 | c ENDDO |
|---|
| 2082 | c ENDDO |
|---|
| 2083 | c print*,"Moment cinetique total = ",mangtot |
|---|
| 2084 | c |
|---|
| 2085 | c------------------------ |
|---|
| 2086 | c |
|---|
| 2087 | c Sauvegarder les valeurs de t et u a la fin de la physique: |
|---|
| 2088 | c |
|---|
| 2089 | DO k = 1, klev |
|---|
| 2090 | DO i = 1, klon |
|---|
| 2091 | u_ancien(i,k) = u_seri(i,k) |
|---|
| 2092 | t_ancien(i,k) = t_seri(i,k) |
|---|
| 2093 | ENDDO |
|---|
| 2094 | ENDDO |
|---|
| 2095 | c |
|---|
| 2096 | c============================================================= |
|---|
| 2097 | c Ecriture des sorties |
|---|
| 2098 | c============================================================= |
|---|
| 2099 | #ifndef MESOSCALE |
|---|
| 2100 | #ifdef CPP_IOIPSL |
|---|
| 2101 | |
|---|
| 2102 | #ifdef histhf |
|---|
| 2103 | #include "write_histhf.h" |
|---|
| 2104 | #endif |
|---|
| 2105 | |
|---|
| 2106 | #ifdef histday |
|---|
| 2107 | #include "write_histday.h" |
|---|
| 2108 | #endif |
|---|
| 2109 | |
|---|
| 2110 | #ifdef histmth |
|---|
| 2111 | #include "write_histmth.h" |
|---|
| 2112 | #endif |
|---|
| 2113 | |
|---|
| 2114 | #ifdef histins |
|---|
| 2115 | #include "write_histins.h" |
|---|
| 2116 | #endif |
|---|
| 2117 | |
|---|
| 2118 | #endif |
|---|
| 2119 | |
|---|
| 2120 | ! XIOS outputs |
|---|
| 2121 | ! This can be done ANYWHERE in the physics routines ! |
|---|
| 2122 | |
|---|
| 2123 | #ifdef CPP_XIOS |
|---|
| 2124 | ! Send fields to XIOS: (NB these fields must also be defined as |
|---|
| 2125 | ! <field id="..." /> in context_lmdz_physics.xml to be correctly used) |
|---|
| 2126 | |
|---|
| 2127 | ! 2D fields |
|---|
| 2128 | |
|---|
| 2129 | CALL send_xios_field("phis",pphis) |
|---|
| 2130 | cell_area_out(:)=cell_area(:) |
|---|
| 2131 | if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon |
|---|
| 2132 | if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon |
|---|
| 2133 | CALL send_xios_field("aire",cell_area_out) |
|---|
| 2134 | CALL send_xios_field("tsol",ftsol) |
|---|
| 2135 | CALL send_xios_field("psol",paprs(:,1)) |
|---|
| 2136 | CALL send_xios_field("cdragh",cdragh) |
|---|
| 2137 | CALL send_xios_field("cdragm",cdragm) |
|---|
| 2138 | |
|---|
| 2139 | CALL send_xios_field("tops",topsw) |
|---|
| 2140 | CALL send_xios_field("topl",toplw) |
|---|
| 2141 | CALL send_xios_field("sols",solsw) |
|---|
| 2142 | CALL send_xios_field("soll",sollw) |
|---|
| 2143 | |
|---|
| 2144 | ! 3D fields |
|---|
| 2145 | |
|---|
| 2146 | CALL send_xios_field("temp",t_seri) |
|---|
| 2147 | CALL send_xios_field("pres",pplay) |
|---|
| 2148 | CALL send_xios_field("geop",zphi) |
|---|
| 2149 | CALL send_xios_field("vitu",u_seri) |
|---|
| 2150 | c VENUS: regardee a l envers!!!!!!!!!!!!!!! |
|---|
| 2151 | CALL send_xios_field("vitv",-1.*v_seri) |
|---|
| 2152 | CALL send_xios_field("vitw",omega) |
|---|
| 2153 | CALL send_xios_field("vitwz",vertwind) |
|---|
| 2154 | CALL send_xios_field("Kz",ycoefh) |
|---|
| 2155 | CALL send_xios_field("mmean",mmean) |
|---|
| 2156 | CALL send_xios_field("Cp",cpnew) |
|---|
| 2157 | CALL send_xios_field("rho",rho) |
|---|
| 2158 | CALL send_xios_field("BV2",zn2) |
|---|
| 2159 | |
|---|
| 2160 | CALL send_xios_field("dudyn",d_u_dyn) |
|---|
| 2161 | CALL send_xios_field("duvdf",d_u_vdf) |
|---|
| 2162 | c VENUS: regardee a l envers!!!!!!!!!!!!!!! |
|---|
| 2163 | CALL send_xios_field("dvvdf",-1.*d_v_vdf) |
|---|
| 2164 | CALL send_xios_field("duajs",d_u_ajs) |
|---|
| 2165 | CALL send_xios_field("dugwo",d_u_oro) |
|---|
| 2166 | CALL send_xios_field("dugwno",d_u_hin) |
|---|
| 2167 | CALL send_xios_field("dvgwno",-1.*d_v_hin) |
|---|
| 2168 | CALL send_xios_field("dumolvis",d_u_molvis) |
|---|
| 2169 | c VENUS: regardee a l envers!!!!!!!!!!!!!!! |
|---|
| 2170 | CALL send_xios_field("dvmolvis",-1.*d_v_molvis) |
|---|
| 2171 | CALL send_xios_field("dtdyn",d_t_dyn) |
|---|
| 2172 | CALL send_xios_field("dtphy",d_t) |
|---|
| 2173 | CALL send_xios_field("dtvdf",d_t_vdf) |
|---|
| 2174 | CALL send_xios_field("dtajs",d_t_ajs) |
|---|
| 2175 | CALL send_xios_field("dtswr",dtsw) |
|---|
| 2176 | CALL send_xios_field("dtswrNLTE",d_t_nirco2) |
|---|
| 2177 | CALL send_xios_field("dtswrLTE",heat) |
|---|
| 2178 | CALL send_xios_field("dtlwr",dtlw) |
|---|
| 2179 | CALL send_xios_field("dtlwrNLTE",d_t_nlte) |
|---|
| 2180 | CALL send_xios_field("dtlwrLTE",-1.*cool) |
|---|
| 2181 | CALL send_xios_field("dteuv",d_t_euv) |
|---|
| 2182 | CALL send_xios_field("dtcond",d_t_conduc) |
|---|
| 2183 | CALL send_xios_field("dtec",d_t_ec) |
|---|
| 2184 | |
|---|
| 2185 | CALL send_xios_field("SWnet",swnet(:,1:klev)) |
|---|
| 2186 | CALL send_xios_field("LWnet",lwnet(:,1:klev)) |
|---|
| 2187 | CALL send_xios_field("fluxvdf",fluxt) |
|---|
| 2188 | CALL send_xios_field("fluxdyn",flux_dyn) |
|---|
| 2189 | CALL send_xios_field("fluxajs",flux_ajs) |
|---|
| 2190 | CALL send_xios_field("fluxec",flux_ec) |
|---|
| 2191 | |
|---|
| 2192 | ! when using tracers |
|---|
| 2193 | |
|---|
| 2194 | if (iflag_trac == 1) then |
|---|
| 2195 | |
|---|
| 2196 | if (ok_aoa) then |
|---|
| 2197 | call send_xios_field("age",age) |
|---|
| 2198 | call send_xios_field("aoa",tr_seri(:,:,i_aoa)) |
|---|
| 2199 | endif |
|---|
| 2200 | |
|---|
| 2201 | ! production and destruction rate, cm-3.s-1 |
|---|
| 2202 | ! Beware of the context*.xml file !! |
|---|
| 2203 | if ((tr_scheme == 3) .and. (ok_chem)) THEN |
|---|
| 2204 | do iq = 1,nqmax - nmicro |
|---|
| 2205 | if ((iq.eq.i_o).or.(iq.eq.i_co)) THEN |
|---|
| 2206 | call send_xios_field("prod_"//tname(iq), |
|---|
| 2207 | $ prod_tr(:,:,iq)) |
|---|
| 2208 | call send_xios_field("loss_"//tname(iq), |
|---|
| 2209 | $ loss_tr(:,:,iq)) |
|---|
| 2210 | end if |
|---|
| 2211 | !if (iq.eq.i_o) then |
|---|
| 2212 | ! call send_xios_field('prod_o', prod_tr(:,:,iq)) |
|---|
| 2213 | ! call send_xios_field('loss_o', loss_tr(:,:,iq)) |
|---|
| 2214 | !end if |
|---|
| 2215 | !if (iq.eq.i_co) then |
|---|
| 2216 | ! call send_xios_field('prod_co', prod_tr(:,:,iq)) |
|---|
| 2217 | ! call send_xios_field('loss_co', loss_tr(:,:,iq)) |
|---|
| 2218 | !end if |
|---|
| 2219 | end do |
|---|
| 2220 | end if |
|---|
| 2221 | |
|---|
| 2222 | ! tracers in gas phase, volume mixing ratio |
|---|
| 2223 | |
|---|
| 2224 | do iq = 1,nqmax - nmicro |
|---|
| 2225 | call send_xios_field(tname(iq), |
|---|
| 2226 | $ tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq)) |
|---|
| 2227 | end do |
|---|
| 2228 | |
|---|
| 2229 | ! tracers in gas phase, column densities |
|---|
| 2230 | |
|---|
| 2231 | do iq = 1,nqmax - nmicro |
|---|
| 2232 | col_dens_tr(:,iq)=0. |
|---|
| 2233 | if (type_tr(iq).eq.1) THEN |
|---|
| 2234 | do k = 1, klev |
|---|
| 2235 | col_dens_tr(:,iq) = col_dens_tr(:,iq) + |
|---|
| 2236 | $ tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG |
|---|
| 2237 | end do |
|---|
| 2238 | call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq)) |
|---|
| 2239 | end if |
|---|
| 2240 | end do |
|---|
| 2241 | |
|---|
| 2242 | ! tracers in liquid phase, volume mixing ratio |
|---|
| 2243 | |
|---|
| 2244 | if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN |
|---|
| 2245 | call send_xios_field(tname(i_h2oliq), |
|---|
| 2246 | $ tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq)) |
|---|
| 2247 | call send_xios_field(tname(i_h2so4liq), |
|---|
| 2248 | $ tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq)) |
|---|
| 2249 | if (ok_sedim) then |
|---|
| 2250 | call send_xios_field("Fsedim",fsedim(:,1:klev)) |
|---|
| 2251 | end if |
|---|
| 2252 | end if |
|---|
| 2253 | |
|---|
| 2254 | ! aeronomical emissions |
|---|
| 2255 | |
|---|
| 2256 | call send_xios_field("no_emis",no_emission) |
|---|
| 2257 | call send_xios_field("o2_emis",o2_emission) |
|---|
| 2258 | |
|---|
| 2259 | ! chemical iterations |
|---|
| 2260 | |
|---|
| 2261 | if (tr_scheme.eq.3) call send_xios_field("iter",real(iter)) |
|---|
| 2262 | |
|---|
| 2263 | end if |
|---|
| 2264 | |
|---|
| 2265 | IF (callthermos .and. ok_chem) THEN |
|---|
| 2266 | CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2)) |
|---|
| 2267 | CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o)) |
|---|
| 2268 | CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2)) |
|---|
| 2269 | ENDIF |
|---|
| 2270 | |
|---|
| 2271 | !! DEBUG |
|---|
| 2272 | ! if (is_master) print*,"itauphy=",itap |
|---|
| 2273 | ! if (itap.eq.10) lafin=.true. |
|---|
| 2274 | |
|---|
| 2275 | if (lafin.and.is_omp_master) then |
|---|
| 2276 | write(*,*) "physiq: call xios_context_finalize" |
|---|
| 2277 | call xios_context_finalize |
|---|
| 2278 | endif |
|---|
| 2279 | |
|---|
| 2280 | #endif |
|---|
| 2281 | #else |
|---|
| 2282 | ! Outputs MESOSCALE |
|---|
| 2283 | CALL allocate_comm_wrf(klon,klev) |
|---|
| 2284 | comm_HR_SW(1:klon,1:klev) = dtsw(1:klon,1:klev) |
|---|
| 2285 | comm_HR_LW(1:klon,1:klev) = dtlw(1:klon,1:klev) |
|---|
| 2286 | comm_DT_RAD(1:klon,1:klev) = d_t_rad(1:klon,1:klev) |
|---|
| 2287 | IF (turb_resolved) THEN |
|---|
| 2288 | open(17,file='hrdyn.txt',form='formatted',status='old') |
|---|
| 2289 | rewind(17) |
|---|
| 2290 | DO k=1,klev |
|---|
| 2291 | read(17,*) dt_dyn(k) |
|---|
| 2292 | ENDDO |
|---|
| 2293 | close(17) |
|---|
| 2294 | |
|---|
| 2295 | do i=1,klon |
|---|
| 2296 | d_t(i,:)=d_t(i,:)+dt_dyn(:) |
|---|
| 2297 | comm_HR_DYN(i,:) = dt_dyn(:) |
|---|
| 2298 | enddo |
|---|
| 2299 | ELSE |
|---|
| 2300 | comm_HR_DYN(1:klon,1:klev) = d_t_dyn(1:klon,1:klev) |
|---|
| 2301 | comm_DT_VDF(1:klon,1:klev) = d_t_vdf(1:klon,1:klev) |
|---|
| 2302 | comm_DT_AJS(1:klon,1:klev) = d_t_ajs(1:klon,1:klev) |
|---|
| 2303 | ENDIF |
|---|
| 2304 | comm_DT(1:klon,1:klev)=d_t(1:klon,1:klev) |
|---|
| 2305 | #endif |
|---|
| 2306 | |
|---|
| 2307 | |
|---|
| 2308 | c==================================================================== |
|---|
| 2309 | c Si c'est la fin, il faut conserver l'etat de redemarrage |
|---|
| 2310 | c==================================================================== |
|---|
| 2311 | c |
|---|
| 2312 | IF (lafin) THEN |
|---|
| 2313 | itau_phy = itau_phy + itap |
|---|
| 2314 | CALL phyredem ("restartphy.nc") |
|---|
| 2315 | |
|---|
| 2316 | c--------------FLOTT |
|---|
| 2317 | CMODEB LOTT |
|---|
| 2318 | C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES |
|---|
| 2319 | C DU BILAN DE MOMENT ANGULAIRE. |
|---|
| 2320 | if (bilansmc.eq.1) then |
|---|
| 2321 | write(*,*)'Fermeture de aaam_bud.out (FL Vous parle)' |
|---|
| 2322 | close(27) |
|---|
| 2323 | close(28) |
|---|
| 2324 | endif !bilansmc |
|---|
| 2325 | CMODFIN |
|---|
| 2326 | c------------- |
|---|
| 2327 | c--------------SLEBONNOIS |
|---|
| 2328 | C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES |
|---|
| 2329 | C DES BALLONS |
|---|
| 2330 | if (ballons.eq.1) then |
|---|
| 2331 | write(*,*)'Fermeture des ballons*.out' |
|---|
| 2332 | close(30) |
|---|
| 2333 | close(31) |
|---|
| 2334 | close(32) |
|---|
| 2335 | close(33) |
|---|
| 2336 | close(34) |
|---|
| 2337 | endif !ballons |
|---|
| 2338 | c------------- |
|---|
| 2339 | ENDIF |
|---|
| 2340 | |
|---|
| 2341 | END SUBROUTINE physiq |
|---|
| 2342 | |
|---|
| 2343 | END MODULE physiq_mod |
|---|
| 2344 | |
|---|
