- Timestamp:
- Jul 23, 2024, 5:57:06 PM (3 months ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_1dutils.f90
r5103 r5104 1 ! $Id$ 2 3 INCLUDE "conf_gcm.f90" 4 5 SUBROUTINE conf_unicol 6 7 use IOIPSL 8 USE print_control_mod, ONLY: lunout 9 IMPLICIT NONE 10 !----------------------------------------------------------------------- 11 ! Auteurs : A. Lahellec . 12 13 ! Declarations : 14 ! -------------- 15 16 include "compar1d.h" 17 include "flux_arp.h" 18 include "tsoilnudge.h" 19 include "fcg_gcssold.h" 20 #include "fcg_racmo.h" 21 include "fcg_racmo.h" 22 23 24 ! local: 25 ! ------ 26 27 ! CHARACTER ch1*72,ch2*72,ch3*72,ch4*12 28 29 ! ------------------------------------------------------------------- 30 31 ! ......... Initilisation parametres du lmdz1D .......... 32 33 !--------------------------------------------------------------------- 34 ! initialisations: 35 ! ---------------- 36 37 !Config Key = lunout 38 !Config Desc = unite de fichier pour les impressions 39 !Config Def = 6 40 !Config Help = unite de fichier pour les impressions 41 !Config (defaut sortie standard = 6) 42 lunout = 6 43 ! CALL getin('lunout', lunout) 44 IF (lunout /= 5 .and. lunout /= 6) THEN 45 OPEN(lunout, FILE = 'lmdz.out') 46 ENDIF 47 48 !Config Key = prt_level 49 !Config Desc = niveau d'impressions de debogage 50 !Config Def = 0 51 !Config Help = Niveau d'impression pour le debogage 52 !Config (0 = minimum d'impression) 53 ! prt_level = 0 54 ! CALL getin('prt_level',prt_level) 55 56 !----------------------------------------------------------------------- 57 ! Parametres de controle du run: 58 !----------------------------------------------------------------------- 59 60 !Config Key = restart 61 !Config Desc = on repart des startphy et start1dyn 62 !Config Def = false 63 !Config Help = les fichiers restart doivent etre renomme en start 64 restart = .FALSE. 65 CALL getin('restart', restart) 66 67 !Config Key = forcing_type 68 !Config Desc = defines the way the SCM is forced: 69 !Config Def = 0 70 !!Config Help = 0 ==> forcing_les = .TRUE. 71 ! initial profiles from file prof.inp.001 72 ! no forcing by LS convergence ; 73 ! surface temperature imposed ; 74 ! radiative cooling may be imposed (iflag_radia=0 in physiq.def) 75 ! = 1 ==> forcing_radconv = .TRUE. 76 ! idem forcing_type = 0, but the imposed radiative cooling 77 ! is set to 0 (hence, if iflag_radia=0 in physiq.def, 78 ! then there is no radiative cooling at all) 79 ! = 2 ==> forcing_toga = .TRUE. 80 ! initial profiles from TOGA-COARE IFA files 81 ! LS convergence and SST imposed from TOGA-COARE IFA files 82 ! = 3 ==> forcing_GCM2SCM = .TRUE. 83 ! initial profiles from the GCM output 84 ! LS convergence imposed from the GCM output 85 ! = 4 ==> forcing_twpi = .TRUE. 86 ! initial profiles from TWPICE nc files 87 ! LS convergence and SST imposed from TWPICE nc files 88 ! = 5 ==> forcing_rico = .TRUE. 89 ! initial profiles from RICO idealized 90 ! LS convergence imposed from RICO (cst) 91 ! = 6 ==> forcing_amma = .TRUE. 92 ! = 10 ==> forcing_case = .TRUE. 93 ! initial profiles from case.nc file 94 ! = 40 ==> forcing_GCSSold = .TRUE. 95 ! initial profile from GCSS file 96 ! LS convergence imposed from GCSS file 97 ! = 50 ==> forcing_fire = .TRUE. 98 ! = 59 ==> forcing_sandu = .TRUE. 99 ! initial profiles from sanduref file: see prof.inp.001 100 ! SST varying with time and divergence constante: see ifa_sanduref.txt file 101 ! Radiation has to be computed interactively 102 ! = 60 ==> forcing_astex = .TRUE. 103 ! initial profiles from file: see prof.inp.001 104 ! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file 105 ! Radiation has to be computed interactively 106 ! = 61 ==> forcing_armcu = .TRUE. 107 ! initial profiles from file: see prof.inp.001 108 ! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt 109 ! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt 110 ! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 111 ! Radiation to be switched off 112 ! > 100 ==> forcing_case = .TRUE. or forcing_case2 = .TRUE. 113 ! initial profiles from case.nc file 114 115 forcing_type = 0 116 CALL getin('forcing_type', forcing_type) 117 imp_fcg_gcssold = .FALSE. 118 ts_fcg_gcssold = .FALSE. 119 Tp_fcg_gcssold = .FALSE. 120 Tp_ini_gcssold = .FALSE. 121 xTurb_fcg_gcssold = .FALSE. 122 IF (forcing_type ==40) THEN 123 CALL getin('imp_fcg', imp_fcg_gcssold) 124 CALL getin('ts_fcg', ts_fcg_gcssold) 125 CALL getin('tp_fcg', Tp_fcg_gcssold) 126 CALL getin('tp_ini', Tp_ini_gcssold) 127 CALL getin('turb_fcg', xTurb_fcg_gcssold) 128 ENDIF 129 130 !Parametres de forcage 131 !Config Key = tend_t 132 !Config Desc = forcage ou non par advection de T 133 !Config Def = false 134 !Config Help = forcage ou non par advection de T 135 tend_t = 0 136 CALL getin('tend_t', tend_t) 137 138 !Config Key = tend_q 139 !Config Desc = forcage ou non par advection de q 140 !Config Def = false 141 !Config Help = forcage ou non par advection de q 142 tend_q = 0 143 CALL getin('tend_q', tend_q) 144 145 !Config Key = tend_u 146 !Config Desc = forcage ou non par advection de u 147 !Config Def = false 148 !Config Help = forcage ou non par advection de u 149 tend_u = 0 150 CALL getin('tend_u', tend_u) 151 152 !Config Key = tend_v 153 !Config Desc = forcage ou non par advection de v 154 !Config Def = false 155 !Config Help = forcage ou non par advection de v 156 tend_v = 0 157 CALL getin('tend_v', tend_v) 158 159 !Config Key = tend_w 160 !Config Desc = forcage ou non par vitesse verticale 161 !Config Def = false 162 !Config Help = forcage ou non par vitesse verticale 163 tend_w = 0 164 CALL getin('tend_w', tend_w) 165 166 !Config Key = tend_rayo 167 !Config Desc = forcage ou non par dtrad 168 !Config Def = false 169 !Config Help = forcage ou non par dtrad 170 tend_rayo = 0 171 CALL getin('tend_rayo', tend_rayo) 172 173 174 !Config Key = nudge_t 175 !Config Desc = constante de nudging de T 176 !Config Def = false 177 !Config Help = constante de nudging de T 178 nudge_t = 0. 179 CALL getin('nudge_t', nudge_t) 180 181 !Config Key = nudge_q 182 !Config Desc = constante de nudging de q 183 !Config Def = false 184 !Config Help = constante de nudging de q 185 nudge_q = 0. 186 CALL getin('nudge_q', nudge_q) 187 188 !Config Key = nudge_u 189 !Config Desc = constante de nudging de u 190 !Config Def = false 191 !Config Help = constante de nudging de u 192 nudge_u = 0. 193 CALL getin('nudge_u', nudge_u) 194 195 !Config Key = nudge_v 196 !Config Desc = constante de nudging de v 197 !Config Def = false 198 !Config Help = constante de nudging de v 199 nudge_v = 0. 200 CALL getin('nudge_v', nudge_v) 201 202 !Config Key = nudge_w 203 !Config Desc = constante de nudging de w 204 !Config Def = false 205 !Config Help = constante de nudging de w 206 nudge_w = 0. 207 CALL getin('nudge_w', nudge_w) 208 209 210 !Config Key = iflag_nudge 211 !Config Desc = atmospheric nudging ttype (decimal code) 212 !Config Def = 0 213 !Config Help = 0 ==> no nudging 214 ! If digit number n of iflag_nudge is set, then nudging of type n is on 215 ! If digit number n of iflag_nudge is not set, then nudging of type n is off 216 ! (digits are numbered from the right) 217 iflag_nudge = 0 218 CALL getin('iflag_nudge', iflag_nudge) 219 220 !Config Key = ok_flux_surf 221 !Config Desc = forcage ou non par les flux de surface 222 !Config Def = false 223 !Config Help = forcage ou non par les flux de surface 224 ok_flux_surf = .FALSE. 225 CALL getin('ok_flux_surf', ok_flux_surf) 226 227 !Config Key = ok_forc_tsurf 228 !Config Desc = forcage ou non par la Ts 229 !Config Def = false 230 !Config Help = forcage ou non par la Ts 231 ok_forc_tsurf = .FALSE. 232 CALL getin('ok_forc_tsurf', ok_forc_tsurf) 233 234 !Config Key = ok_prescr_ust 235 !Config Desc = ustar impose ou non 236 !Config Def = false 237 !Config Help = ustar impose ou non 238 ok_prescr_ust = .FALSE. 239 CALL getin('ok_prescr_ust', ok_prescr_ust) 240 241 242 !Config Key = ok_prescr_beta 243 !Config Desc = betaevap impose ou non 244 !Config Def = false 245 !Config Help = betaevap impose ou non 246 ok_prescr_beta = .FALSE. 247 CALL getin('ok_prescr_beta', ok_prescr_beta) 248 249 !Config Key = ok_old_disvert 250 !Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) 251 !Config Def = false 252 !Config Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) 253 ok_old_disvert = .FALSE. 254 CALL getin('ok_old_disvert', ok_old_disvert) 255 256 !Config Key = time_ini 257 !Config Desc = meaningless in this case 258 !Config Def = 0. 259 !Config Help = 260 time_ini = 0. 261 CALL getin('time_ini', time_ini) 262 263 !Config Key = rlat et rlon 264 !Config Desc = latitude et longitude 265 !Config Def = 0.0 0.0 266 !Config Help = fixe la position de la colonne 267 xlat = 0. 268 xlon = 0. 269 CALL getin('rlat', xlat) 270 CALL getin('rlon', xlon) 271 272 !Config Key = airephy 273 !Config Desc = Grid cell area 274 !Config Def = 1.e11 275 !Config Help = 276 airefi = 1.e11 277 CALL getin('airephy', airefi) 278 279 !Config Key = nat_surf 280 !Config Desc = surface type 281 !Config Def = 0 (ocean) 282 !Config Help = 0=ocean,1=land,2=glacier,3=banquise 283 nat_surf = 0. 284 CALL getin('nat_surf', nat_surf) 285 286 !Config Key = tsurf 287 !Config Desc = surface temperature 288 !Config Def = 290. 289 !Config Help = surface temperature 290 tsurf = 290. 291 CALL getin('tsurf', tsurf) 292 293 !Config Key = psurf 294 !Config Desc = surface pressure 295 !Config Def = 102400. 296 !Config Help = 297 psurf = 102400. 298 CALL getin('psurf', psurf) 299 300 !Config Key = zsurf 301 !Config Desc = surface altitude 302 !Config Def = 0. 303 !Config Help = 304 zsurf = 0. 305 CALL getin('zsurf', zsurf) 306 ! EV pour accord avec format standard 307 CALL getin('zorog', zsurf) 308 309 310 !Config Key = rugos 311 !Config Desc = coefficient de frottement 312 !Config Def = 0.0001 313 !Config Help = calcul du Cdrag 314 rugos = 0.0001 315 CALL getin('rugos', rugos) 316 ! FH/2020/04/08/confinement: Pour le nouveau format standard, la rugosite s'appelle z0 317 CALL getin('z0', rugos) 318 319 !Config Key = rugosh 320 !Config Desc = coefficient de frottement 321 !Config Def = rugos 322 !Config Help = calcul du Cdrag 323 rugosh = rugos 324 CALL getin('rugosh', rugosh) 325 326 327 328 !Config Key = snowmass 329 !Config Desc = mass de neige de la surface en kg/m2 330 !Config Def = 0.0000 331 !Config Help = snowmass 332 snowmass = 0.0000 333 CALL getin('snowmass', snowmass) 334 335 !Config Key = wtsurf et wqsurf 336 !Config Desc = ??? 337 !Config Def = 0.0 0.0 338 !Config Help = 339 wtsurf = 0.0 340 wqsurf = 0.0 341 CALL getin('wtsurf', wtsurf) 342 CALL getin('wqsurf', wqsurf) 343 344 !Config Key = albedo 345 !Config Desc = albedo 346 !Config Def = 0.09 347 !Config Help = 348 albedo = 0.09 349 CALL getin('albedo', albedo) 350 351 !Config Key = agesno 352 !Config Desc = age de la neige 353 !Config Def = 30.0 354 !Config Help = 355 xagesno = 30.0 356 CALL getin('agesno', xagesno) 357 358 !Config Key = restart_runoff 359 !Config Desc = age de la neige 360 !Config Def = 30.0 361 !Config Help = 362 restart_runoff = 0.0 363 CALL getin('restart_runoff', restart_runoff) 364 365 !Config Key = qsolinp 366 !Config Desc = initial bucket water content (kg/m2) when land (5std) 367 !Config Def = 30.0 368 !Config Help = 369 qsolinp = 1. 370 CALL getin('qsolinp', qsolinp) 371 372 373 374 !Config Key = betaevap 375 !Config Desc = beta for actual evaporation when prescribed 376 !Config Def = 1.0 377 !Config Help = 378 betaevap = 1. 379 CALL getin('betaevap', betaevap) 380 381 !Config Key = zpicinp 382 !Config Desc = denivellation orographie 383 !Config Def = 0. 384 !Config Help = input brise 385 zpicinp = 0. 386 CALL getin('zpicinp', zpicinp) 387 !Config key = nudge_tsoil 388 !Config Desc = activation of soil temperature nudging 389 !Config Def = .FALSE. 390 !Config Help = ... 391 392 nudge_tsoil = .FALSE. 393 CALL getin('nudge_tsoil', nudge_tsoil) 394 395 !Config key = isoil_nudge 396 !Config Desc = level number where soil temperature is nudged 397 !Config Def = 3 398 !Config Help = ... 399 400 isoil_nudge = 3 401 CALL getin('isoil_nudge', isoil_nudge) 402 403 !Config key = Tsoil_nudge 404 !Config Desc = target temperature for tsoil(isoil_nudge) 405 !Config Def = 300. 406 !Config Help = ... 407 408 Tsoil_nudge = 300. 409 CALL getin('Tsoil_nudge', Tsoil_nudge) 410 411 !Config key = tau_soil_nudge 412 !Config Desc = nudging relaxation time for tsoil 413 !Config Def = 3600. 414 !Config Help = ... 415 416 tau_soil_nudge = 3600. 417 CALL getin('tau_soil_nudge', tau_soil_nudge) 418 419 !---------------------------------------------------------- 420 ! Param??tres de for??age pour les forcages communs: 421 ! Pour les forcages communs: ces entiers valent 0 ou 1 422 ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale 423 ! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale 424 ! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv) 425 ! forcages en omega, w, vent geostrophique ou ustar 426 ! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging 427 !---------------------------------------------------------- 428 429 !Config Key = tadv 430 !Config Desc = forcage ou non par advection totale de T 431 !Config Def = false 432 !Config Help = forcage ou non par advection totale de T 433 tadv = 0 434 CALL getin('tadv', tadv) 435 436 !Config Key = tadvv 437 !Config Desc = forcage ou non par advection verticale de T 438 !Config Def = false 439 !Config Help = forcage ou non par advection verticale de T 440 tadvv = 0 441 CALL getin('tadvv', tadvv) 442 443 !Config Key = tadvh 444 !Config Desc = forcage ou non par advection horizontale de T 445 !Config Def = false 446 !Config Help = forcage ou non par advection horizontale de T 447 tadvh = 0 448 CALL getin('tadvh', tadvh) 449 450 !Config Key = thadv 451 !Config Desc = forcage ou non par advection totale de Theta 452 !Config Def = false 453 !Config Help = forcage ou non par advection totale de Theta 454 thadv = 0 455 CALL getin('thadv', thadv) 456 457 !Config Key = thadvv 458 !Config Desc = forcage ou non par advection verticale de Theta 459 !Config Def = false 460 !Config Help = forcage ou non par advection verticale de Theta 461 thadvv = 0 462 CALL getin('thadvv', thadvv) 463 464 !Config Key = thadvh 465 !Config Desc = forcage ou non par advection horizontale de Theta 466 !Config Def = false 467 !Config Help = forcage ou non par advection horizontale de Theta 468 thadvh = 0 469 CALL getin('thadvh', thadvh) 470 471 !Config Key = qadv 472 !Config Desc = forcage ou non par advection totale de Q 473 !Config Def = false 474 !Config Help = forcage ou non par advection totale de Q 475 qadv = 0 476 CALL getin('qadv', qadv) 477 478 !Config Key = qadvv 479 !Config Desc = forcage ou non par advection verticale de Q 480 !Config Def = false 481 !Config Help = forcage ou non par advection verticale de Q 482 qadvv = 0 483 CALL getin('qadvv', qadvv) 484 485 !Config Key = qadvh 486 !Config Desc = forcage ou non par advection horizontale de Q 487 !Config Def = false 488 !Config Help = forcage ou non par advection horizontale de Q 489 qadvh = 0 490 CALL getin('qadvh', qadvh) 491 492 !Config Key = trad 493 !Config Desc = forcage ou non par tendance radiative 494 !Config Def = false 495 !Config Help = forcage ou non par tendance radiative 496 trad = 0 497 CALL getin('trad', trad) 498 499 !Config Key = forc_omega 500 !Config Desc = forcage ou non par omega 501 !Config Def = false 502 !Config Help = forcage ou non par omega 503 forc_omega = 0 504 CALL getin('forc_omega', forc_omega) 505 506 !Config Key = forc_u 507 !Config Desc = forcage ou non par u 508 !Config Def = false 509 !Config Help = forcage ou non par u 510 forc_u = 0 511 CALL getin('forc_u', forc_u) 512 513 !Config Key = forc_v 514 !Config Desc = forcage ou non par v 515 !Config Def = false 516 !Config Help = forcage ou non par v 517 forc_v = 0 518 CALL getin('forc_v', forc_v) 519 !Config Key = forc_w 520 !Config Desc = forcage ou non par w 521 !Config Def = false 522 !Config Help = forcage ou non par w 523 forc_w = 0 524 CALL getin('forc_w', forc_w) 525 526 !Config Key = forc_geo 527 !Config Desc = forcage ou non par geo 528 !Config Def = false 529 !Config Help = forcage ou non par geo 530 forc_geo = 0 531 CALL getin('forc_geo', forc_geo) 532 533 ! Meme chose que ok_precr_ust 534 !Config Key = forc_ustar 535 !Config Desc = forcage ou non par ustar 536 !Config Def = false 537 !Config Help = forcage ou non par ustar 538 forc_ustar = 0 539 CALL getin('forc_ustar', forc_ustar) 540 IF (forc_ustar == 1) ok_prescr_ust = .TRUE. 541 542 543 !Config Key = nudging_u 544 !Config Desc = forcage ou non par nudging sur u 545 !Config Def = false 546 !Config Help = forcage ou non par nudging sur u 547 nudging_u = 0 548 CALL getin('nudging_u', nudging_u) 549 550 !Config Key = nudging_v 551 !Config Desc = forcage ou non par nudging sur v 552 !Config Def = false 553 !Config Help = forcage ou non par nudging sur v 554 nudging_v = 0 555 CALL getin('nudging_v', nudging_v) 556 557 !Config Key = nudging_w 558 !Config Desc = forcage ou non par nudging sur w 559 !Config Def = false 560 !Config Help = forcage ou non par nudging sur w 561 nudging_w = 0 562 CALL getin('nudging_w', nudging_w) 563 564 ! RELIQUE ANCIENS FORMAT. ECRASE PAR LE SUIVANT 565 !Config Key = nudging_q 566 !Config Desc = forcage ou non par nudging sur q 567 !Config Def = false 568 !Config Help = forcage ou non par nudging sur q 569 nudging_qv = 0 570 CALL getin('nudging_q', nudging_qv) 571 CALL getin('nudging_qv', nudging_qv) 572 573 p_nudging_u = 11000. 574 p_nudging_v = 11000. 575 p_nudging_t = 11000. 576 p_nudging_qv = 11000. 577 CALL getin('p_nudging_u', p_nudging_u) 578 CALL getin('p_nudging_v', p_nudging_v) 579 CALL getin('p_nudging_t', p_nudging_t) 580 CALL getin('p_nudging_qv', p_nudging_qv) 581 582 !Config Key = nudging_t 583 !Config Desc = forcage ou non par nudging sur t 584 !Config Def = false 585 !Config Help = forcage ou non par nudging sur t 586 nudging_t = 0 587 CALL getin('nudging_t', nudging_t) 588 589 write(lunout, *)' +++++++++++++++++++++++++++++++++++++++' 590 write(lunout, *)' Configuration des parametres du gcm1D: ' 591 write(lunout, *)' +++++++++++++++++++++++++++++++++++++++' 592 write(lunout, *)' restart = ', restart 593 write(lunout, *)' forcing_type = ', forcing_type 594 write(lunout, *)' time_ini = ', time_ini 595 write(lunout, *)' rlat = ', xlat 596 write(lunout, *)' rlon = ', xlon 597 write(lunout, *)' airephy = ', airefi 598 write(lunout, *)' nat_surf = ', nat_surf 599 write(lunout, *)' tsurf = ', tsurf 600 write(lunout, *)' psurf = ', psurf 601 write(lunout, *)' zsurf = ', zsurf 602 write(lunout, *)' rugos = ', rugos 603 write(lunout, *)' snowmass=', snowmass 604 write(lunout, *)' wtsurf = ', wtsurf 605 write(lunout, *)' wqsurf = ', wqsurf 606 write(lunout, *)' albedo = ', albedo 607 write(lunout, *)' xagesno = ', xagesno 608 write(lunout, *)' restart_runoff = ', restart_runoff 609 write(lunout, *)' qsolinp = ', qsolinp 610 write(lunout, *)' zpicinp = ', zpicinp 611 write(lunout, *)' nudge_tsoil = ', nudge_tsoil 612 write(lunout, *)' isoil_nudge = ', isoil_nudge 613 write(lunout, *)' Tsoil_nudge = ', Tsoil_nudge 614 write(lunout, *)' tau_soil_nudge = ', tau_soil_nudge 615 write(lunout, *)' tadv = ', tadv 616 write(lunout, *)' tadvv = ', tadvv 617 write(lunout, *)' tadvh = ', tadvh 618 write(lunout, *)' thadv = ', thadv 619 write(lunout, *)' thadvv = ', thadvv 620 write(lunout, *)' thadvh = ', thadvh 621 write(lunout, *)' qadv = ', qadv 622 write(lunout, *)' qadvv = ', qadvv 623 write(lunout, *)' qadvh = ', qadvh 624 write(lunout, *)' trad = ', trad 625 write(lunout, *)' forc_omega = ', forc_omega 626 write(lunout, *)' forc_w = ', forc_w 627 write(lunout, *)' forc_geo = ', forc_geo 628 write(lunout, *)' forc_ustar = ', forc_ustar 629 write(lunout, *)' nudging_u = ', nudging_u 630 write(lunout, *)' nudging_v = ', nudging_v 631 write(lunout, *)' nudging_t = ', nudging_t 632 write(lunout, *)' nudging_qv = ', nudging_qv 633 IF (forcing_type ==40) THEN 634 write(lunout, *) '--- Forcing type GCSS Old --- with:' 635 write(lunout, *)'imp_fcg', imp_fcg_gcssold 636 write(lunout, *)'ts_fcg', ts_fcg_gcssold 637 write(lunout, *)'tp_fcg', Tp_fcg_gcssold 638 write(lunout, *)'tp_ini', Tp_ini_gcssold 639 write(lunout, *)'xturb_fcg', xTurb_fcg_gcssold 640 ENDIF 641 642 write(lunout, *)' +++++++++++++++++++++++++++++++++++++++' 643 write(lunout, *) 644 645 RETURN 646 END 647 648 ! $Id: dyn1deta0.F 1279 2010/07/30 A Lahellec$ 649 650 651 SUBROUTINE dyn1deta0(fichnom, plev, play, phi, phis, presnivs, & 652 & ucov, vcov, temp, q, omega2) 653 USE dimphy 654 USE mod_grid_phy_lmdz 655 USE mod_phys_lmdz_para 656 USE iophy 657 USE phys_state_var_mod 658 USE iostart 659 USE write_field_phy 660 USE infotrac 661 use control_mod 662 USE comconst_mod, ONLY: im, jm, lllm 663 USE logic_mod, ONLY: fxyhypb, ysinus 664 USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn 665 USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr 666 667 IMPLICIT NONE 668 !======================================================= 669 ! Ecriture du fichier de redemarrage sous format NetCDF 670 !======================================================= 671 ! Declarations: 672 ! ------------- 673 include "dimensions.h" 674 !!#include "control.h" 675 676 ! Arguments: 677 ! ---------- 678 CHARACTER*(*) fichnom 679 !Al1 plev tronque pour .nc mais plev(klev+1):=0 680 real :: plev(klon, klev + 1), play (klon, klev), phi(klon, klev) 681 real :: presnivs(klon, klev) 682 real :: ucov(klon, klev), vcov(klon, klev), temp(klon, klev) 683 real :: q(klon, klev, nqtot), omega2(klon, klev) 684 ! real :: ug(klev),vg(klev),fcoriolis 685 real :: phis(klon) 686 687 ! Variables locales pour NetCDF: 688 ! ------------------------------ 689 INTEGER iq 690 INTEGER length 691 PARAMETER (length = 100) 692 REAL tab_cntrl(length) ! tableau des parametres du run 693 character*4 nmq(nqtot) 694 character*12 modname 695 character*80 abort_message 696 LOGICAL found 697 698 modname = 'dyn1deta0 : ' 699 !! nmq(1)="vap" 700 !! nmq(2)="cond" 701 !! do iq=3,nqtot 702 !! write(nmq(iq),'("tra",i1)') iq-2 703 !! enddo 704 DO iq = 1, nqtot 705 nmq(iq) = trim(tracers(iq)%name) 706 ENDDO 707 PRINT*, 'in dyn1deta0 ', fichnom, klon, klev, nqtot 708 CALL open_startphy(fichnom) 709 PRINT*, 'after open startphy ', fichnom, nmq 710 711 ! Lecture des parametres de controle: 712 713 CALL get_var("controle", tab_cntrl) 714 715 im = tab_cntrl(1) 716 jm = tab_cntrl(2) 717 lllm = tab_cntrl(3) 718 day_ref = tab_cntrl(4) 719 annee_ref = tab_cntrl(5) 720 ! rad = tab_cntrl(6) 721 ! omeg = tab_cntrl(7) 722 ! g = tab_cntrl(8) 723 ! cpp = tab_cntrl(9) 724 ! kappa = tab_cntrl(10) 725 ! daysec = tab_cntrl(11) 726 ! dtvr = tab_cntrl(12) 727 ! etot0 = tab_cntrl(13) 728 ! ptot0 = tab_cntrl(14) 729 ! ztot0 = tab_cntrl(15) 730 ! stot0 = tab_cntrl(16) 731 ! ang0 = tab_cntrl(17) 732 ! pa = tab_cntrl(18) 733 ! preff = tab_cntrl(19) 734 735 ! clon = tab_cntrl(20) 736 ! clat = tab_cntrl(21) 737 ! grossismx = tab_cntrl(22) 738 ! grossismy = tab_cntrl(23) 739 740 IF (tab_cntrl(24)==1.) THEN 741 fxyhypb = .TRUE. 742 ! dzoomx = tab_cntrl(25) 743 ! dzoomy = tab_cntrl(26) 744 ! taux = tab_cntrl(28) 745 ! tauy = tab_cntrl(29) 746 ELSE 747 fxyhypb = .FALSE. 748 ysinus = .FALSE. 749 IF(tab_cntrl(27)==1.) ysinus = .TRUE. 750 ENDIF 751 752 day_ini = tab_cntrl(30) 753 itau_dyn = tab_cntrl(31) 754 ! ................................................................. 755 756 757 ! PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa 758 !Al1 759 Print*, 'day_ref,annee_ref,day_ini,itau_dyn', & 760 & day_ref, annee_ref, day_ini, itau_dyn 761 762 ! Lecture des champs 763 764 CALL get_field("play", play, found) 765 IF (.NOT. found) PRINT*, modname // 'Le champ <Play> est absent' 766 CALL get_field("phi", phi, found) 767 IF (.NOT. found) PRINT*, modname // 'Le champ <Phi> est absent' 768 CALL get_field("phis", phis, found) 769 IF (.NOT. found) PRINT*, modname // 'Le champ <Phis> est absent' 770 CALL get_field("presnivs", presnivs, found) 771 IF (.NOT. found) PRINT*, modname // 'Le champ <Presnivs> est absent' 772 CALL get_field("ucov", ucov, found) 773 IF (.NOT. found) PRINT*, modname // 'Le champ <ucov> est absent' 774 CALL get_field("vcov", vcov, found) 775 IF (.NOT. found) PRINT*, modname // 'Le champ <vcov> est absent' 776 CALL get_field("temp", temp, found) 777 IF (.NOT. found) PRINT*, modname // 'Le champ <temp> est absent' 778 CALL get_field("omega2", omega2, found) 779 IF (.NOT. found) PRINT*, modname // 'Le champ <omega2> est absent' 780 plev(1, klev + 1) = 0. 781 CALL get_field("plev", plev(:, 1:klev), found) 782 IF (.NOT. found) PRINT*, modname // 'Le champ <Plev> est absent' 783 784 Do iq = 1, nqtot 785 CALL get_field("q" // nmq(iq), q(:, :, iq), found) 786 IF (.NOT.found)PRINT*, modname // 'Le champ <q' // nmq // '> est absent' 787 EndDo 788 789 CALL close_startphy 790 PRINT*, ' close startphy', fichnom, play(1, 1), play(1, klev), temp(1, klev) 791 792 RETURN 793 END 794 795 ! $Id: dyn1dredem.F 1279 2010/07/29 A Lahellec$ 796 797 798 SUBROUTINE dyn1dredem(fichnom, plev, play, phi, phis, presnivs, & 799 & ucov, vcov, temp, q, omega2) 800 USE dimphy 801 USE mod_grid_phy_lmdz 802 USE mod_phys_lmdz_para 803 USE phys_state_var_mod 804 USE iostart 805 USE infotrac 806 use control_mod 807 USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad 808 USE logic_mod, ONLY: fxyhypb, ysinus 809 USE temps_mod, ONLY: annee_ref, day_end, day_ref, itau_dyn, itaufin 810 USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr 811 812 IMPLICIT NONE 813 !======================================================= 814 ! Ecriture du fichier de redemarrage sous format NetCDF 815 !======================================================= 816 ! Declarations: 817 ! ------------- 818 include "dimensions.h" 819 !!#include "control.h" 820 821 ! Arguments: 822 ! ---------- 823 CHARACTER*(*) fichnom 824 !Al1 plev tronque pour .nc mais plev(klev+1):=0 825 real :: plev(klon, klev), play (klon, klev), phi(klon, klev) 826 real :: presnivs(klon, klev) 827 real :: ucov(klon, klev), vcov(klon, klev), temp(klon, klev) 828 real :: q(klon, klev, nqtot) 829 real :: omega2(klon, klev), rho(klon, klev + 1) 830 ! real :: ug(klev),vg(klev),fcoriolis 831 real :: phis(klon) 832 833 ! Variables locales pour NetCDF: 834 ! ------------------------------ 835 INTEGER nid 836 INTEGER ierr 837 INTEGER iq, l 838 INTEGER length 839 PARAMETER (length = 100) 840 REAL tab_cntrl(length) ! tableau des parametres du run 841 character*4 nmq(nqtot) 842 character*20 modname 843 character*80 abort_message 844 845 INTEGER pass 846 847 CALL open_restartphy(fichnom) 848 PRINT*, 'redm1 ', fichnom, klon, klev, nqtot 849 !! nmq(1)="vap" 850 !! nmq(2)="cond" 851 !! nmq(3)="tra1" 852 !! nmq(4)="tra2" 853 DO iq = 1, nqtot 854 nmq(iq) = trim(tracers(iq)%name) 855 ENDDO 856 857 ! modname = 'dyn1dredem' 858 ! ierr = nf90_open(fichnom, nf90_write, nid) 859 ! IF (ierr .NE. nf90_noerr) THEN 860 ! abort_message="Pb. d ouverture "//fichnom 861 ! CALL abort_gcm('Modele 1D',abort_message,1) 862 ! ENDIF 863 864 DO l = 1, length 865 tab_cntrl(l) = 0. 866 ENDDO 867 tab_cntrl(1) = FLOAT(iim) 868 tab_cntrl(2) = FLOAT(jjm) 869 tab_cntrl(3) = FLOAT(llm) 870 tab_cntrl(4) = FLOAT(day_ref) 871 tab_cntrl(5) = FLOAT(annee_ref) 872 tab_cntrl(6) = rad 873 tab_cntrl(7) = omeg 874 tab_cntrl(8) = g 875 tab_cntrl(9) = cpp 876 tab_cntrl(10) = kappa 877 tab_cntrl(11) = daysec 878 tab_cntrl(12) = dtvr 879 ! tab_cntrl(13) = etot0 880 ! tab_cntrl(14) = ptot0 881 ! tab_cntrl(15) = ztot0 882 ! tab_cntrl(16) = stot0 883 ! tab_cntrl(17) = ang0 884 ! tab_cntrl(18) = pa 885 ! tab_cntrl(19) = preff 886 887 ! ..... parametres pour le zoom ...... 888 889 ! tab_cntrl(20) = clon 890 ! tab_cntrl(21) = clat 891 ! tab_cntrl(22) = grossismx 892 ! tab_cntrl(23) = grossismy 893 894 IF (fxyhypb) THEN 895 tab_cntrl(24) = 1. 896 ! tab_cntrl(25) = dzoomx 897 ! tab_cntrl(26) = dzoomy 898 tab_cntrl(27) = 0. 899 ! tab_cntrl(28) = taux 900 ! tab_cntrl(29) = tauy 901 ELSE 902 tab_cntrl(24) = 0. 903 ! tab_cntrl(25) = dzoomx 904 ! tab_cntrl(26) = dzoomy 905 tab_cntrl(27) = 0. 906 tab_cntrl(28) = 0. 907 tab_cntrl(29) = 0. 908 IF(ysinus) tab_cntrl(27) = 1. 909 ENDIF 910 !Al1 iday_end -> day_end 911 tab_cntrl(30) = FLOAT(day_end) 912 tab_cntrl(31) = FLOAT(itau_dyn + itaufin) 913 914 DO pass = 1, 2 915 CALL put_var(pass, "controle", "Param. de controle Dyn1D", tab_cntrl) 916 917 ! Ecriture/extension de la coordonnee temps 918 919 920 ! Ecriture des champs 921 922 CALL put_field(pass, "plev", "p interfaces sauf la nulle", plev) 923 CALL put_field(pass, "play", "", play) 924 CALL put_field(pass, "phi", "geopotentielle", phi) 925 CALL put_field(pass, "phis", "geopotentiell de surface", phis) 926 CALL put_field(pass, "presnivs", "", presnivs) 927 CALL put_field(pass, "ucov", "", ucov) 928 CALL put_field(pass, "vcov", "", vcov) 929 CALL put_field(pass, "temp", "", temp) 930 CALL put_field(pass, "omega2", "", omega2) 1 MODULE lmdz_1dutils 2 IMPLICIT NONE; PRIVATE 3 PUBLIC fq_sat, conf_unicol, dyn1deta0, dyn1dredem, gr_fi_dyn, abort_gcm, gr_dyn_fi, & 4 disvert0, advect_vert, advect_va, lstendh, nudge_rht_init, nudge_uv_init, & 5 nudge_rht, nudge_uv, interp2_case_vertical 6 CONTAINS 7 REAL FUNCTION fq_sat(kelvin, millibar) 8 IMPLICIT none 9 !====================================================================== 10 ! Autheur(s): Z.X. Li (LMD/CNRS) 11 ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.) 12 !====================================================================== 13 ! Arguments: 14 ! kelvin---input-R: temperature en Kelvin 15 ! millibar--input-R: pression en mb 16 17 ! fq_sat----output-R: vapeur d'eau saturante en kg/kg 18 !====================================================================== 19 20 REAL, INTENT(IN) :: kelvin, millibar 21 22 REAL r2es 23 PARAMETER (r2es = 611.14 * 18.0153 / 28.9644) 24 REAL r3les, r3ies, r3es 25 PARAMETER (R3LES = 17.269) 26 PARAMETER (R3IES = 21.875) 27 28 REAL r4les, r4ies, r4es 29 PARAMETER (R4LES = 35.86) 30 PARAMETER (R4IES = 7.66) 31 32 REAL rtt 33 PARAMETER (rtt = 273.16) 34 35 REAL retv 36 PARAMETER (retv = 28.9644 / 18.0153 - 1.0) 37 38 REAL zqsat 39 REAL temp, pres 40 ! ------------------------------------------------------------------ 41 42 temp = kelvin 43 pres = millibar * 100.0 44 ! write(*,*)'kelvin,millibar=',kelvin,millibar 45 ! write(*,*)'temp,pres=',temp,pres 46 47 IF (temp <= rtt) THEN 48 r3es = r3ies 49 r4es = r4ies 50 ELSE 51 r3es = r3les 52 r4es = r4les 53 ENDIF 54 55 zqsat = r2es / pres * EXP (r3es * (temp - rtt) / (temp - r4es)) 56 zqsat = MIN(0.5, ZQSAT) 57 zqsat = zqsat / (1. - retv * zqsat) 58 59 fq_sat = zqsat 60 END FUNCTION fq_sat 61 62 SUBROUTINE conf_unicol 63 64 use IOIPSL 65 USE print_control_mod, ONLY: lunout 66 !----------------------------------------------------------------------- 67 ! Auteurs : A. Lahellec . 68 69 ! Declarations : 70 ! -------------- 71 72 include "compar1d.h" 73 include "flux_arp.h" 74 include "tsoilnudge.h" 75 include "fcg_gcssold.h" 76 include "fcg_racmo.h" 77 78 79 ! local: 80 ! ------ 81 82 ! CHARACTER ch1*72,ch2*72,ch3*72,ch4*12 83 84 ! ------------------------------------------------------------------- 85 86 ! ......... Initilisation parametres du lmdz1D .......... 87 88 !--------------------------------------------------------------------- 89 ! initialisations: 90 ! ---------------- 91 92 !Config Key = lunout 93 !Config Desc = unite de fichier pour les impressions 94 !Config Def = 6 95 !Config Help = unite de fichier pour les impressions 96 !Config (defaut sortie standard = 6) 97 lunout = 6 98 ! CALL getin('lunout', lunout) 99 IF (lunout /= 5 .and. lunout /= 6) THEN 100 OPEN(lunout, FILE = 'lmdz.out') 101 ENDIF 102 103 !Config Key = prt_level 104 !Config Desc = niveau d'impressions de debogage 105 !Config Def = 0 106 !Config Help = Niveau d'impression pour le debogage 107 !Config (0 = minimum d'impression) 108 ! prt_level = 0 109 ! CALL getin('prt_level',prt_level) 110 111 !----------------------------------------------------------------------- 112 ! Parametres de controle du run: 113 !----------------------------------------------------------------------- 114 115 !Config Key = restart 116 !Config Desc = on repart des startphy et start1dyn 117 !Config Def = false 118 !Config Help = les fichiers restart doivent etre renomme en start 119 restart = .FALSE. 120 CALL getin('restart', restart) 121 122 !Config Key = forcing_type 123 !Config Desc = defines the way the SCM is forced: 124 !Config Def = 0 125 !!Config Help = 0 ==> forcing_les = .TRUE. 126 ! initial profiles from file prof.inp.001 127 ! no forcing by LS convergence ; 128 ! surface temperature imposed ; 129 ! radiative cooling may be imposed (iflag_radia=0 in physiq.def) 130 ! = 1 ==> forcing_radconv = .TRUE. 131 ! idem forcing_type = 0, but the imposed radiative cooling 132 ! is set to 0 (hence, if iflag_radia=0 in physiq.def, 133 ! then there is no radiative cooling at all) 134 ! = 2 ==> forcing_toga = .TRUE. 135 ! initial profiles from TOGA-COARE IFA files 136 ! LS convergence and SST imposed from TOGA-COARE IFA files 137 ! = 3 ==> forcing_GCM2SCM = .TRUE. 138 ! initial profiles from the GCM output 139 ! LS convergence imposed from the GCM output 140 ! = 4 ==> forcing_twpi = .TRUE. 141 ! initial profiles from TWPICE nc files 142 ! LS convergence and SST imposed from TWPICE nc files 143 ! = 5 ==> forcing_rico = .TRUE. 144 ! initial profiles from RICO idealized 145 ! LS convergence imposed from RICO (cst) 146 ! = 6 ==> forcing_amma = .TRUE. 147 ! = 10 ==> forcing_case = .TRUE. 148 ! initial profiles from case.nc file 149 ! = 40 ==> forcing_GCSSold = .TRUE. 150 ! initial profile from GCSS file 151 ! LS convergence imposed from GCSS file 152 ! = 50 ==> forcing_fire = .TRUE. 153 ! = 59 ==> forcing_sandu = .TRUE. 154 ! initial profiles from sanduref file: see prof.inp.001 155 ! SST varying with time and divergence constante: see ifa_sanduref.txt file 156 ! Radiation has to be computed interactively 157 ! = 60 ==> forcing_astex = .TRUE. 158 ! initial profiles from file: see prof.inp.001 159 ! SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file 160 ! Radiation has to be computed interactively 161 ! = 61 ==> forcing_armcu = .TRUE. 162 ! initial profiles from file: see prof.inp.001 163 ! sensible and latent heat flux imposed: see ifa_arm_cu_1.txt 164 ! large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt 165 ! use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 166 ! Radiation to be switched off 167 ! > 100 ==> forcing_case = .TRUE. or forcing_case2 = .TRUE. 168 ! initial profiles from case.nc file 169 170 forcing_type = 0 171 CALL getin('forcing_type', forcing_type) 172 imp_fcg_gcssold = .FALSE. 173 ts_fcg_gcssold = .FALSE. 174 Tp_fcg_gcssold = .FALSE. 175 Tp_ini_gcssold = .FALSE. 176 xTurb_fcg_gcssold = .FALSE. 177 IF (forcing_type ==40) THEN 178 CALL getin('imp_fcg', imp_fcg_gcssold) 179 CALL getin('ts_fcg', ts_fcg_gcssold) 180 CALL getin('tp_fcg', Tp_fcg_gcssold) 181 CALL getin('tp_ini', Tp_ini_gcssold) 182 CALL getin('turb_fcg', xTurb_fcg_gcssold) 183 ENDIF 184 185 !Parametres de forcage 186 !Config Key = tend_t 187 !Config Desc = forcage ou non par advection de T 188 !Config Def = false 189 !Config Help = forcage ou non par advection de T 190 tend_t = 0 191 CALL getin('tend_t', tend_t) 192 193 !Config Key = tend_q 194 !Config Desc = forcage ou non par advection de q 195 !Config Def = false 196 !Config Help = forcage ou non par advection de q 197 tend_q = 0 198 CALL getin('tend_q', tend_q) 199 200 !Config Key = tend_u 201 !Config Desc = forcage ou non par advection de u 202 !Config Def = false 203 !Config Help = forcage ou non par advection de u 204 tend_u = 0 205 CALL getin('tend_u', tend_u) 206 207 !Config Key = tend_v 208 !Config Desc = forcage ou non par advection de v 209 !Config Def = false 210 !Config Help = forcage ou non par advection de v 211 tend_v = 0 212 CALL getin('tend_v', tend_v) 213 214 !Config Key = tend_w 215 !Config Desc = forcage ou non par vitesse verticale 216 !Config Def = false 217 !Config Help = forcage ou non par vitesse verticale 218 tend_w = 0 219 CALL getin('tend_w', tend_w) 220 221 !Config Key = tend_rayo 222 !Config Desc = forcage ou non par dtrad 223 !Config Def = false 224 !Config Help = forcage ou non par dtrad 225 tend_rayo = 0 226 CALL getin('tend_rayo', tend_rayo) 227 228 229 !Config Key = nudge_t 230 !Config Desc = constante de nudging de T 231 !Config Def = false 232 !Config Help = constante de nudging de T 233 nudge_t = 0. 234 CALL getin('nudge_t', nudge_t) 235 236 !Config Key = nudge_q 237 !Config Desc = constante de nudging de q 238 !Config Def = false 239 !Config Help = constante de nudging de q 240 nudge_q = 0. 241 CALL getin('nudge_q', nudge_q) 242 243 !Config Key = nudge_u 244 !Config Desc = constante de nudging de u 245 !Config Def = false 246 !Config Help = constante de nudging de u 247 nudge_u = 0. 248 CALL getin('nudge_u', nudge_u) 249 250 !Config Key = nudge_v 251 !Config Desc = constante de nudging de v 252 !Config Def = false 253 !Config Help = constante de nudging de v 254 nudge_v = 0. 255 CALL getin('nudge_v', nudge_v) 256 257 !Config Key = nudge_w 258 !Config Desc = constante de nudging de w 259 !Config Def = false 260 !Config Help = constante de nudging de w 261 nudge_w = 0. 262 CALL getin('nudge_w', nudge_w) 263 264 265 !Config Key = iflag_nudge 266 !Config Desc = atmospheric nudging ttype (decimal code) 267 !Config Def = 0 268 !Config Help = 0 ==> no nudging 269 ! If digit number n of iflag_nudge is set, then nudging of type n is on 270 ! If digit number n of iflag_nudge is not set, then nudging of type n is off 271 ! (digits are numbered from the right) 272 iflag_nudge = 0 273 CALL getin('iflag_nudge', iflag_nudge) 274 275 !Config Key = ok_flux_surf 276 !Config Desc = forcage ou non par les flux de surface 277 !Config Def = false 278 !Config Help = forcage ou non par les flux de surface 279 ok_flux_surf = .FALSE. 280 CALL getin('ok_flux_surf', ok_flux_surf) 281 282 !Config Key = ok_forc_tsurf 283 !Config Desc = forcage ou non par la Ts 284 !Config Def = false 285 !Config Help = forcage ou non par la Ts 286 ok_forc_tsurf = .FALSE. 287 CALL getin('ok_forc_tsurf', ok_forc_tsurf) 288 289 !Config Key = ok_prescr_ust 290 !Config Desc = ustar impose ou non 291 !Config Def = false 292 !Config Help = ustar impose ou non 293 ok_prescr_ust = .FALSE. 294 CALL getin('ok_prescr_ust', ok_prescr_ust) 295 296 297 !Config Key = ok_prescr_beta 298 !Config Desc = betaevap impose ou non 299 !Config Def = false 300 !Config Help = betaevap impose ou non 301 ok_prescr_beta = .FALSE. 302 CALL getin('ok_prescr_beta', ok_prescr_beta) 303 304 !Config Key = ok_old_disvert 305 !Config Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) 306 !Config Def = false 307 !Config Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h) 308 ok_old_disvert = .FALSE. 309 CALL getin('ok_old_disvert', ok_old_disvert) 310 311 !Config Key = time_ini 312 !Config Desc = meaningless in this case 313 !Config Def = 0. 314 !Config Help = 315 time_ini = 0. 316 CALL getin('time_ini', time_ini) 317 318 !Config Key = rlat et rlon 319 !Config Desc = latitude et longitude 320 !Config Def = 0.0 0.0 321 !Config Help = fixe la position de la colonne 322 xlat = 0. 323 xlon = 0. 324 CALL getin('rlat', xlat) 325 CALL getin('rlon', xlon) 326 327 !Config Key = airephy 328 !Config Desc = Grid cell area 329 !Config Def = 1.e11 330 !Config Help = 331 airefi = 1.e11 332 CALL getin('airephy', airefi) 333 334 !Config Key = nat_surf 335 !Config Desc = surface type 336 !Config Def = 0 (ocean) 337 !Config Help = 0=ocean,1=land,2=glacier,3=banquise 338 nat_surf = 0. 339 CALL getin('nat_surf', nat_surf) 340 341 !Config Key = tsurf 342 !Config Desc = surface temperature 343 !Config Def = 290. 344 !Config Help = surface temperature 345 tsurf = 290. 346 CALL getin('tsurf', tsurf) 347 348 !Config Key = psurf 349 !Config Desc = surface pressure 350 !Config Def = 102400. 351 !Config Help = 352 psurf = 102400. 353 CALL getin('psurf', psurf) 354 355 !Config Key = zsurf 356 !Config Desc = surface altitude 357 !Config Def = 0. 358 !Config Help = 359 zsurf = 0. 360 CALL getin('zsurf', zsurf) 361 ! EV pour accord avec format standard 362 CALL getin('zorog', zsurf) 363 364 365 !Config Key = rugos 366 !Config Desc = coefficient de frottement 367 !Config Def = 0.0001 368 !Config Help = calcul du Cdrag 369 rugos = 0.0001 370 CALL getin('rugos', rugos) 371 ! FH/2020/04/08/confinement: Pour le nouveau format standard, la rugosite s'appelle z0 372 CALL getin('z0', rugos) 373 374 !Config Key = rugosh 375 !Config Desc = coefficient de frottement 376 !Config Def = rugos 377 !Config Help = calcul du Cdrag 378 rugosh = rugos 379 CALL getin('rugosh', rugosh) 380 381 382 383 !Config Key = snowmass 384 !Config Desc = mass de neige de la surface en kg/m2 385 !Config Def = 0.0000 386 !Config Help = snowmass 387 snowmass = 0.0000 388 CALL getin('snowmass', snowmass) 389 390 !Config Key = wtsurf et wqsurf 391 !Config Desc = ??? 392 !Config Def = 0.0 0.0 393 !Config Help = 394 wtsurf = 0.0 395 wqsurf = 0.0 396 CALL getin('wtsurf', wtsurf) 397 CALL getin('wqsurf', wqsurf) 398 399 !Config Key = albedo 400 !Config Desc = albedo 401 !Config Def = 0.09 402 !Config Help = 403 albedo = 0.09 404 CALL getin('albedo', albedo) 405 406 !Config Key = agesno 407 !Config Desc = age de la neige 408 !Config Def = 30.0 409 !Config Help = 410 xagesno = 30.0 411 CALL getin('agesno', xagesno) 412 413 !Config Key = restart_runoff 414 !Config Desc = age de la neige 415 !Config Def = 30.0 416 !Config Help = 417 restart_runoff = 0.0 418 CALL getin('restart_runoff', restart_runoff) 419 420 !Config Key = qsolinp 421 !Config Desc = initial bucket water content (kg/m2) when land (5std) 422 !Config Def = 30.0 423 !Config Help = 424 qsolinp = 1. 425 CALL getin('qsolinp', qsolinp) 426 427 428 429 !Config Key = betaevap 430 !Config Desc = beta for actual evaporation when prescribed 431 !Config Def = 1.0 432 !Config Help = 433 betaevap = 1. 434 CALL getin('betaevap', betaevap) 435 436 !Config Key = zpicinp 437 !Config Desc = denivellation orographie 438 !Config Def = 0. 439 !Config Help = input brise 440 zpicinp = 0. 441 CALL getin('zpicinp', zpicinp) 442 !Config key = nudge_tsoil 443 !Config Desc = activation of soil temperature nudging 444 !Config Def = .FALSE. 445 !Config Help = ... 446 447 nudge_tsoil = .FALSE. 448 CALL getin('nudge_tsoil', nudge_tsoil) 449 450 !Config key = isoil_nudge 451 !Config Desc = level number where soil temperature is nudged 452 !Config Def = 3 453 !Config Help = ... 454 455 isoil_nudge = 3 456 CALL getin('isoil_nudge', isoil_nudge) 457 458 !Config key = Tsoil_nudge 459 !Config Desc = target temperature for tsoil(isoil_nudge) 460 !Config Def = 300. 461 !Config Help = ... 462 463 Tsoil_nudge = 300. 464 CALL getin('Tsoil_nudge', Tsoil_nudge) 465 466 !Config key = tau_soil_nudge 467 !Config Desc = nudging relaxation time for tsoil 468 !Config Def = 3600. 469 !Config Help = ... 470 471 tau_soil_nudge = 3600. 472 CALL getin('tau_soil_nudge', tau_soil_nudge) 473 474 !---------------------------------------------------------- 475 ! Param??tres de for??age pour les forcages communs: 476 ! Pour les forcages communs: ces entiers valent 0 ou 1 477 ! tadv= advection tempe, tadvv= adv tempe verticale, tadvh= adv tempe horizontale 478 ! qadv= advection q, qadvv= adv q verticale, qadvh= adv q horizontale 479 ! trad= 0 (rayonnement actif) ou 1 (prescrit par tend_rad) ou adv (prescir et contenu dans les tadv) 480 ! forcages en omega, w, vent geostrophique ou ustar 481 ! Parametres de nudging en u,v,t,q valent 0 ou 1 ou le temps de nudging 482 !---------------------------------------------------------- 483 484 !Config Key = tadv 485 !Config Desc = forcage ou non par advection totale de T 486 !Config Def = false 487 !Config Help = forcage ou non par advection totale de T 488 tadv = 0 489 CALL getin('tadv', tadv) 490 491 !Config Key = tadvv 492 !Config Desc = forcage ou non par advection verticale de T 493 !Config Def = false 494 !Config Help = forcage ou non par advection verticale de T 495 tadvv = 0 496 CALL getin('tadvv', tadvv) 497 498 !Config Key = tadvh 499 !Config Desc = forcage ou non par advection horizontale de T 500 !Config Def = false 501 !Config Help = forcage ou non par advection horizontale de T 502 tadvh = 0 503 CALL getin('tadvh', tadvh) 504 505 !Config Key = thadv 506 !Config Desc = forcage ou non par advection totale de Theta 507 !Config Def = false 508 !Config Help = forcage ou non par advection totale de Theta 509 thadv = 0 510 CALL getin('thadv', thadv) 511 512 !Config Key = thadvv 513 !Config Desc = forcage ou non par advection verticale de Theta 514 !Config Def = false 515 !Config Help = forcage ou non par advection verticale de Theta 516 thadvv = 0 517 CALL getin('thadvv', thadvv) 518 519 !Config Key = thadvh 520 !Config Desc = forcage ou non par advection horizontale de Theta 521 !Config Def = false 522 !Config Help = forcage ou non par advection horizontale de Theta 523 thadvh = 0 524 CALL getin('thadvh', thadvh) 525 526 !Config Key = qadv 527 !Config Desc = forcage ou non par advection totale de Q 528 !Config Def = false 529 !Config Help = forcage ou non par advection totale de Q 530 qadv = 0 531 CALL getin('qadv', qadv) 532 533 !Config Key = qadvv 534 !Config Desc = forcage ou non par advection verticale de Q 535 !Config Def = false 536 !Config Help = forcage ou non par advection verticale de Q 537 qadvv = 0 538 CALL getin('qadvv', qadvv) 539 540 !Config Key = qadvh 541 !Config Desc = forcage ou non par advection horizontale de Q 542 !Config Def = false 543 !Config Help = forcage ou non par advection horizontale de Q 544 qadvh = 0 545 CALL getin('qadvh', qadvh) 546 547 !Config Key = trad 548 !Config Desc = forcage ou non par tendance radiative 549 !Config Def = false 550 !Config Help = forcage ou non par tendance radiative 551 trad = 0 552 CALL getin('trad', trad) 553 554 !Config Key = forc_omega 555 !Config Desc = forcage ou non par omega 556 !Config Def = false 557 !Config Help = forcage ou non par omega 558 forc_omega = 0 559 CALL getin('forc_omega', forc_omega) 560 561 !Config Key = forc_u 562 !Config Desc = forcage ou non par u 563 !Config Def = false 564 !Config Help = forcage ou non par u 565 forc_u = 0 566 CALL getin('forc_u', forc_u) 567 568 !Config Key = forc_v 569 !Config Desc = forcage ou non par v 570 !Config Def = false 571 !Config Help = forcage ou non par v 572 forc_v = 0 573 CALL getin('forc_v', forc_v) 574 !Config Key = forc_w 575 !Config Desc = forcage ou non par w 576 !Config Def = false 577 !Config Help = forcage ou non par w 578 forc_w = 0 579 CALL getin('forc_w', forc_w) 580 581 !Config Key = forc_geo 582 !Config Desc = forcage ou non par geo 583 !Config Def = false 584 !Config Help = forcage ou non par geo 585 forc_geo = 0 586 CALL getin('forc_geo', forc_geo) 587 588 ! Meme chose que ok_precr_ust 589 !Config Key = forc_ustar 590 !Config Desc = forcage ou non par ustar 591 !Config Def = false 592 !Config Help = forcage ou non par ustar 593 forc_ustar = 0 594 CALL getin('forc_ustar', forc_ustar) 595 IF (forc_ustar == 1) ok_prescr_ust = .TRUE. 596 597 598 !Config Key = nudging_u 599 !Config Desc = forcage ou non par nudging sur u 600 !Config Def = false 601 !Config Help = forcage ou non par nudging sur u 602 nudging_u = 0 603 CALL getin('nudging_u', nudging_u) 604 605 !Config Key = nudging_v 606 !Config Desc = forcage ou non par nudging sur v 607 !Config Def = false 608 !Config Help = forcage ou non par nudging sur v 609 nudging_v = 0 610 CALL getin('nudging_v', nudging_v) 611 612 !Config Key = nudging_w 613 !Config Desc = forcage ou non par nudging sur w 614 !Config Def = false 615 !Config Help = forcage ou non par nudging sur w 616 nudging_w = 0 617 CALL getin('nudging_w', nudging_w) 618 619 ! RELIQUE ANCIENS FORMAT. ECRASE PAR LE SUIVANT 620 !Config Key = nudging_q 621 !Config Desc = forcage ou non par nudging sur q 622 !Config Def = false 623 !Config Help = forcage ou non par nudging sur q 624 nudging_qv = 0 625 CALL getin('nudging_q', nudging_qv) 626 CALL getin('nudging_qv', nudging_qv) 627 628 p_nudging_u = 11000. 629 p_nudging_v = 11000. 630 p_nudging_t = 11000. 631 p_nudging_qv = 11000. 632 CALL getin('p_nudging_u', p_nudging_u) 633 CALL getin('p_nudging_v', p_nudging_v) 634 CALL getin('p_nudging_t', p_nudging_t) 635 CALL getin('p_nudging_qv', p_nudging_qv) 636 637 !Config Key = nudging_t 638 !Config Desc = forcage ou non par nudging sur t 639 !Config Def = false 640 !Config Help = forcage ou non par nudging sur t 641 nudging_t = 0 642 CALL getin('nudging_t', nudging_t) 643 644 write(lunout, *)' +++++++++++++++++++++++++++++++++++++++' 645 write(lunout, *)' Configuration des parametres du gcm1D: ' 646 write(lunout, *)' +++++++++++++++++++++++++++++++++++++++' 647 write(lunout, *)' restart = ', restart 648 write(lunout, *)' forcing_type = ', forcing_type 649 write(lunout, *)' time_ini = ', time_ini 650 write(lunout, *)' rlat = ', xlat 651 write(lunout, *)' rlon = ', xlon 652 write(lunout, *)' airephy = ', airefi 653 write(lunout, *)' nat_surf = ', nat_surf 654 write(lunout, *)' tsurf = ', tsurf 655 write(lunout, *)' psurf = ', psurf 656 write(lunout, *)' zsurf = ', zsurf 657 write(lunout, *)' rugos = ', rugos 658 write(lunout, *)' snowmass=', snowmass 659 write(lunout, *)' wtsurf = ', wtsurf 660 write(lunout, *)' wqsurf = ', wqsurf 661 write(lunout, *)' albedo = ', albedo 662 write(lunout, *)' xagesno = ', xagesno 663 write(lunout, *)' restart_runoff = ', restart_runoff 664 write(lunout, *)' qsolinp = ', qsolinp 665 write(lunout, *)' zpicinp = ', zpicinp 666 write(lunout, *)' nudge_tsoil = ', nudge_tsoil 667 write(lunout, *)' isoil_nudge = ', isoil_nudge 668 write(lunout, *)' Tsoil_nudge = ', Tsoil_nudge 669 write(lunout, *)' tau_soil_nudge = ', tau_soil_nudge 670 write(lunout, *)' tadv = ', tadv 671 write(lunout, *)' tadvv = ', tadvv 672 write(lunout, *)' tadvh = ', tadvh 673 write(lunout, *)' thadv = ', thadv 674 write(lunout, *)' thadvv = ', thadvv 675 write(lunout, *)' thadvh = ', thadvh 676 write(lunout, *)' qadv = ', qadv 677 write(lunout, *)' qadvv = ', qadvv 678 write(lunout, *)' qadvh = ', qadvh 679 write(lunout, *)' trad = ', trad 680 write(lunout, *)' forc_omega = ', forc_omega 681 write(lunout, *)' forc_w = ', forc_w 682 write(lunout, *)' forc_geo = ', forc_geo 683 write(lunout, *)' forc_ustar = ', forc_ustar 684 write(lunout, *)' nudging_u = ', nudging_u 685 write(lunout, *)' nudging_v = ', nudging_v 686 write(lunout, *)' nudging_t = ', nudging_t 687 write(lunout, *)' nudging_qv = ', nudging_qv 688 IF (forcing_type ==40) THEN 689 write(lunout, *) '--- Forcing type GCSS Old --- with:' 690 write(lunout, *)'imp_fcg', imp_fcg_gcssold 691 write(lunout, *)'ts_fcg', ts_fcg_gcssold 692 write(lunout, *)'tp_fcg', Tp_fcg_gcssold 693 write(lunout, *)'tp_ini', Tp_ini_gcssold 694 write(lunout, *)'xturb_fcg', xTurb_fcg_gcssold 695 ENDIF 696 697 write(lunout, *)' +++++++++++++++++++++++++++++++++++++++' 698 write(lunout, *) 699 700 END SUBROUTINE conf_unicol 701 702 703 SUBROUTINE dyn1deta0(fichnom, plev, play, phi, phis, presnivs, & 704 & ucov, vcov, temp, q, omega2) 705 USE dimphy 706 USE mod_grid_phy_lmdz 707 USE mod_phys_lmdz_para 708 USE iophy 709 USE phys_state_var_mod 710 USE iostart 711 USE write_field_phy 712 USE infotrac 713 use control_mod 714 USE comconst_mod, ONLY: im, jm, lllm 715 USE logic_mod, ONLY: fxyhypb, ysinus 716 USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn 717 USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr 718 719 IMPLICIT NONE 720 !======================================================= 721 ! Ecriture du fichier de redemarrage sous format NetCDF 722 !======================================================= 723 ! Declarations: 724 ! ------------- 725 include "dimensions.h" 726 727 ! Arguments: 728 ! ---------- 729 CHARACTER*(*) fichnom 730 !Al1 plev tronque pour .nc mais plev(klev+1):=0 731 real :: plev(klon, klev + 1), play (klon, klev), phi(klon, klev) 732 real :: presnivs(klon, klev) 733 real :: ucov(klon, klev), vcov(klon, klev), temp(klon, klev) 734 real :: q(klon, klev, nqtot), omega2(klon, klev) 735 ! real :: ug(klev),vg(klev),fcoriolis 736 real :: phis(klon) 737 738 ! Variables locales pour NetCDF: 739 ! ------------------------------ 740 INTEGER iq 741 INTEGER length 742 PARAMETER (length = 100) 743 REAL tab_cntrl(length) ! tableau des parametres du run 744 character*4 nmq(nqtot) 745 character*12 modname 746 character*80 abort_message 747 LOGICAL found 748 749 modname = 'dyn1deta0 : ' 750 !! nmq(1)="vap" 751 !! nmq(2)="cond" 752 !! do iq=3,nqtot 753 !! write(nmq(iq),'("tra",i1)') iq-2 754 !! enddo 755 DO iq = 1, nqtot 756 nmq(iq) = trim(tracers(iq)%name) 757 ENDDO 758 PRINT*, 'in dyn1deta0 ', fichnom, klon, klev, nqtot 759 CALL open_startphy(fichnom) 760 PRINT*, 'after open startphy ', fichnom, nmq 761 762 ! Lecture des parametres de controle: 763 CALL get_var("controle", tab_cntrl) 764 765 im = tab_cntrl(1) 766 jm = tab_cntrl(2) 767 lllm = tab_cntrl(3) 768 day_ref = tab_cntrl(4) 769 annee_ref = tab_cntrl(5) 770 ! rad = tab_cntrl(6) 771 ! omeg = tab_cntrl(7) 772 ! g = tab_cntrl(8) 773 ! cpp = tab_cntrl(9) 774 ! kappa = tab_cntrl(10) 775 ! daysec = tab_cntrl(11) 776 ! dtvr = tab_cntrl(12) 777 ! etot0 = tab_cntrl(13) 778 ! ptot0 = tab_cntrl(14) 779 ! ztot0 = tab_cntrl(15) 780 ! stot0 = tab_cntrl(16) 781 ! ang0 = tab_cntrl(17) 782 ! pa = tab_cntrl(18) 783 ! preff = tab_cntrl(19) 784 785 ! clon = tab_cntrl(20) 786 ! clat = tab_cntrl(21) 787 ! grossismx = tab_cntrl(22) 788 ! grossismy = tab_cntrl(23) 789 790 IF (tab_cntrl(24)==1.) THEN 791 fxyhypb = .TRUE. 792 ! dzoomx = tab_cntrl(25) 793 ! dzoomy = tab_cntrl(26) 794 ! taux = tab_cntrl(28) 795 ! tauy = tab_cntrl(29) 796 ELSE 797 fxyhypb = .FALSE. 798 ysinus = .FALSE. 799 IF(tab_cntrl(27)==1.) ysinus = .TRUE. 800 ENDIF 801 802 day_ini = tab_cntrl(30) 803 itau_dyn = tab_cntrl(31) 804 805 Print*, 'day_ref,annee_ref,day_ini,itau_dyn', day_ref, annee_ref, day_ini, itau_dyn 806 807 ! Lecture des champs 808 CALL get_field("play", play, found) 809 IF (.NOT. found) PRINT*, modname // 'Le champ <Play> est absent' 810 CALL get_field("phi", phi, found) 811 IF (.NOT. found) PRINT*, modname // 'Le champ <Phi> est absent' 812 CALL get_field("phis", phis, found) 813 IF (.NOT. found) PRINT*, modname // 'Le champ <Phis> est absent' 814 CALL get_field("presnivs", presnivs, found) 815 IF (.NOT. found) PRINT*, modname // 'Le champ <Presnivs> est absent' 816 CALL get_field("ucov", ucov, found) 817 IF (.NOT. found) PRINT*, modname // 'Le champ <ucov> est absent' 818 CALL get_field("vcov", vcov, found) 819 IF (.NOT. found) PRINT*, modname // 'Le champ <vcov> est absent' 820 CALL get_field("temp", temp, found) 821 IF (.NOT. found) PRINT*, modname // 'Le champ <temp> est absent' 822 CALL get_field("omega2", omega2, found) 823 IF (.NOT. found) PRINT*, modname // 'Le champ <omega2> est absent' 824 plev(1, klev + 1) = 0. 825 CALL get_field("plev", plev(:, 1:klev), found) 826 IF (.NOT. found) PRINT*, modname // 'Le champ <Plev> est absent' 931 827 932 828 Do iq = 1, nqtot 933 CALL put_field(pass, "q" // nmq(iq), "eau vap ou condens et traceurs", &934 & q(:, :, iq))829 CALL get_field("q" // nmq(iq), q(:, :, iq), found) 830 IF (.NOT.found)PRINT*, modname // 'Le champ <q' // nmq // '> est absent' 935 831 EndDo 936 IF (pass==1) CALL enddef_restartphy 937 IF (pass==2) CALL close_restartphy 938 939 ENDDO 940 941 RETURN 942 END 943 SUBROUTINE gr_fi_dyn(nfield, ngrid, im, jm, pfi, pdyn) 944 IMPLICIT NONE 945 !======================================================================= 946 ! passage d'un champ de la grille scalaire a la grille physique 947 !======================================================================= 948 949 !----------------------------------------------------------------------- 950 ! declarations: 951 ! ------------- 952 953 INTEGER im, jm, ngrid, nfield 954 REAL pdyn(im, jm, nfield) 955 REAL pfi(ngrid, nfield) 956 957 INTEGER i, j, ifield, ig 958 959 !----------------------------------------------------------------------- 960 ! calcul: 961 ! ------- 962 963 DO ifield = 1, nfield 832 833 CALL close_startphy 834 PRINT*, ' close startphy', fichnom, play(1, 1), play(1, klev), temp(1, klev) 835 END SUBROUTINE dyn1deta0 836 837 838 SUBROUTINE dyn1dredem(fichnom, plev, play, phi, phis, presnivs, & 839 & ucov, vcov, temp, q, omega2) 840 USE dimphy 841 USE mod_grid_phy_lmdz 842 USE mod_phys_lmdz_para 843 USE phys_state_var_mod 844 USE iostart 845 USE infotrac 846 use control_mod 847 USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad 848 USE logic_mod, ONLY: fxyhypb, ysinus 849 USE temps_mod, ONLY: annee_ref, day_end, day_ref, itau_dyn, itaufin 850 USE netcdf, ONLY: nf90_open, nf90_write, nf90_noerr 851 852 IMPLICIT NONE 853 !======================================================= 854 ! Ecriture du fichier de redemarrage sous format NetCDF 855 !======================================================= 856 ! Declarations: 857 ! ------------- 858 include "dimensions.h" 859 860 ! Arguments: 861 ! ---------- 862 CHARACTER*(*) fichnom 863 !Al1 plev tronque pour .nc mais plev(klev+1):=0 864 real :: plev(klon, klev), play (klon, klev), phi(klon, klev) 865 real :: presnivs(klon, klev) 866 real :: ucov(klon, klev), vcov(klon, klev), temp(klon, klev) 867 real :: q(klon, klev, nqtot) 868 real :: omega2(klon, klev), rho(klon, klev + 1) 869 ! real :: ug(klev),vg(klev),fcoriolis 870 real :: phis(klon) 871 872 ! Variables locales pour NetCDF: 873 ! ------------------------------ 874 INTEGER nid 875 INTEGER ierr 876 INTEGER iq, l 877 INTEGER length 878 PARAMETER (length = 100) 879 REAL tab_cntrl(length) ! tableau des parametres du run 880 character*4 nmq(nqtot) 881 character*20 modname 882 character*80 abort_message 883 884 INTEGER pass 885 886 CALL open_restartphy(fichnom) 887 PRINT*, 'redm1 ', fichnom, klon, klev, nqtot 888 !! nmq(1)="vap" 889 !! nmq(2)="cond" 890 !! nmq(3)="tra1" 891 !! nmq(4)="tra2" 892 DO iq = 1, nqtot 893 nmq(iq) = trim(tracers(iq)%name) 894 ENDDO 895 896 ! modname = 'dyn1dredem' 897 ! ierr = nf90_open(fichnom, nf90_write, nid) 898 ! IF (ierr .NE. nf90_noerr) THEN 899 ! abort_message="Pb. d ouverture "//fichnom 900 ! CALL abort_gcm('Modele 1D',abort_message,1) 901 ! ENDIF 902 903 DO l = 1, length 904 tab_cntrl(l) = 0. 905 ENDDO 906 tab_cntrl(1) = FLOAT(iim) 907 tab_cntrl(2) = FLOAT(jjm) 908 tab_cntrl(3) = FLOAT(llm) 909 tab_cntrl(4) = FLOAT(day_ref) 910 tab_cntrl(5) = FLOAT(annee_ref) 911 tab_cntrl(6) = rad 912 tab_cntrl(7) = omeg 913 tab_cntrl(8) = g 914 tab_cntrl(9) = cpp 915 tab_cntrl(10) = kappa 916 tab_cntrl(11) = daysec 917 tab_cntrl(12) = dtvr 918 ! tab_cntrl(13) = etot0 919 ! tab_cntrl(14) = ptot0 920 ! tab_cntrl(15) = ztot0 921 ! tab_cntrl(16) = stot0 922 ! tab_cntrl(17) = ang0 923 ! tab_cntrl(18) = pa 924 ! tab_cntrl(19) = preff 925 926 ! ..... parametres pour le zoom ...... 927 928 ! tab_cntrl(20) = clon 929 ! tab_cntrl(21) = clat 930 ! tab_cntrl(22) = grossismx 931 ! tab_cntrl(23) = grossismy 932 933 IF (fxyhypb) THEN 934 tab_cntrl(24) = 1. 935 ! tab_cntrl(25) = dzoomx 936 ! tab_cntrl(26) = dzoomy 937 tab_cntrl(27) = 0. 938 ! tab_cntrl(28) = taux 939 ! tab_cntrl(29) = tauy 940 ELSE 941 tab_cntrl(24) = 0. 942 ! tab_cntrl(25) = dzoomx 943 ! tab_cntrl(26) = dzoomy 944 tab_cntrl(27) = 0. 945 tab_cntrl(28) = 0. 946 tab_cntrl(29) = 0. 947 IF(ysinus) tab_cntrl(27) = 1. 948 ENDIF 949 !Al1 iday_end -> day_end 950 tab_cntrl(30) = FLOAT(day_end) 951 tab_cntrl(31) = FLOAT(itau_dyn + itaufin) 952 953 DO pass = 1, 2 954 CALL put_var(pass, "controle", "Param. de controle Dyn1D", tab_cntrl) 955 956 ! Ecriture/extension de la coordonnee temps 957 958 959 ! Ecriture des champs 960 961 CALL put_field(pass, "plev", "p interfaces sauf la nulle", plev) 962 CALL put_field(pass, "play", "", play) 963 CALL put_field(pass, "phi", "geopotentielle", phi) 964 CALL put_field(pass, "phis", "geopotentiell de surface", phis) 965 CALL put_field(pass, "presnivs", "", presnivs) 966 CALL put_field(pass, "ucov", "", ucov) 967 CALL put_field(pass, "vcov", "", vcov) 968 CALL put_field(pass, "temp", "", temp) 969 CALL put_field(pass, "omega2", "", omega2) 970 971 Do iq = 1, nqtot 972 CALL put_field(pass, "q" // nmq(iq), "eau vap ou condens et traceurs", & 973 & q(:, :, iq)) 974 EndDo 975 IF (pass==1) CALL enddef_restartphy 976 IF (pass==2) CALL close_restartphy 977 978 ENDDO 979 980 RETURN 981 END SUBROUTINE dyn1dredem 982 983 984 SUBROUTINE gr_fi_dyn(nfield, ngrid, im, jm, pfi, pdyn) 985 IMPLICIT NONE 986 !======================================================================= 987 ! passage d'un champ de la grille scalaire a la grille physique 988 !======================================================================= 989 990 !----------------------------------------------------------------------- 991 ! declarations: 992 ! ------------- 993 994 INTEGER im, jm, ngrid, nfield 995 REAL pdyn(im, jm, nfield) 996 REAL pfi(ngrid, nfield) 997 998 INTEGER i, j, ifield, ig 999 1000 !----------------------------------------------------------------------- 1001 ! calcul: 1002 ! ------- 1003 1004 DO ifield = 1, nfield 1005 ! traitement des poles 1006 DO i = 1, im 1007 pdyn(i, 1, ifield) = pfi(1, ifield) 1008 pdyn(i, jm, ifield) = pfi(ngrid, ifield) 1009 ENDDO 1010 1011 ! traitement des point normaux 1012 DO j = 2, jm - 1 1013 ig = 2 + (j - 2) * (im - 1) 1014 CALL SCOPY(im - 1, pfi(ig, ifield), 1, pdyn(1, j, ifield), 1) 1015 pdyn(im, j, ifield) = pdyn(1, j, ifield) 1016 ENDDO 1017 ENDDO 1018 1019 RETURN 1020 END SUBROUTINE gr_fi_dyn 1021 1022 1023 SUBROUTINE abort_gcm(modname, message, ierr) 1024 USE IOIPSL 1025 1026 ! Stops the simulation cleanly, closing files and printing various 1027 ! comments 1028 1029 ! Input: modname = name of calling program 1030 ! message = stuff to print 1031 ! ierr = severity of situation ( = 0 normal ) 1032 1033 character(len = *) modname 1034 integer ierr 1035 character(len = *) message 1036 1037 write(*, *) 'in abort_gcm' 1038 CALL histclo 1039 ! CALL histclo(2) 1040 ! CALL histclo(3) 1041 ! CALL histclo(4) 1042 ! CALL histclo(5) 1043 write(*, *) 'out of histclo' 1044 write(*, *) 'Stopping in ', modname 1045 write(*, *) 'Reason = ', message 1046 CALL getin_dump 1047 1048 if (ierr == 0) then 1049 write(*, *) 'Everything is cool' 1050 else 1051 write(*, *) 'Houston, we have a problem ', ierr 1052 endif 1053 STOP 1054 END SUBROUTINE abort_gcm 1055 1056 1057 SUBROUTINE gr_dyn_fi(nfield, im, jm, ngrid, pdyn, pfi) 1058 IMPLICIT NONE 1059 !======================================================================= 1060 ! passage d'un champ de la grille scalaire a la grille physique 1061 !======================================================================= 1062 1063 !----------------------------------------------------------------------- 1064 ! declarations: 1065 ! ------------- 1066 1067 INTEGER im, jm, ngrid, nfield 1068 REAL pdyn(im, jm, nfield) 1069 REAL pfi(ngrid, nfield) 1070 1071 INTEGER j, ifield, ig 1072 1073 !----------------------------------------------------------------------- 1074 ! calcul: 1075 ! ------- 1076 1077 IF(ngrid/=2 + (jm - 2) * (im - 1).AND.ngrid/=1) & 1078 & STOP 'probleme de dim' 964 1079 ! traitement des poles 965 DO i = 1, im 966 pdyn(i, 1, ifield) = pfi(1, ifield) 967 pdyn(i, jm, ifield) = pfi(ngrid, ifield) 1080 CALL SCOPY(nfield, pdyn, im * jm, pfi, ngrid) 1081 CALL SCOPY(nfield, pdyn(1, jm, 1), im * jm, pfi(ngrid, 1), ngrid) 1082 1083 ! traitement des point normaux 1084 DO ifield = 1, nfield 1085 DO j = 2, jm - 1 1086 ig = 2 + (j - 2) * (im - 1) 1087 CALL SCOPY(im - 1, pdyn(1, j, ifield), 1, pfi(ig, ifield), 1) 1088 ENDDO 968 1089 ENDDO 969 970 ! traitement des point normaux 971 DO j = 2, jm - 1 972 ig = 2 + (j - 2) * (im - 1) 973 CALL SCOPY(im - 1, pfi(ig, ifield), 1, pdyn(1, j, ifield), 1) 974 pdyn(im, j, ifield) = pdyn(1, j, ifield) 1090 END SUBROUTINE gr_dyn_fi 1091 1092 1093 SUBROUTINE disvert0(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) 1094 1095 ! Ancienne version disvert dont on a modifie nom pour utiliser 1096 ! le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes) 1097 ! (MPL 18092012) 1098 1099 ! Auteur : P. Le Van . 1100 1101 IMPLICIT NONE 1102 1103 include "dimensions.h" 1104 include "paramet.h" 1105 1106 !======================================================================= 1107 1108 1109 ! s = sigma ** kappa : coordonnee verticale 1110 ! dsig(l) : epaisseur de la couche l ds la coord. s 1111 ! sig(l) : sigma a l'interface des couches l et l-1 1112 ! ds(l) : distance entre les couches l et l-1 en coord.s 1113 1114 !======================================================================= 1115 1116 REAL pa, preff 1117 REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) 1118 REAL presnivs(llm) 1119 1120 ! declarations: 1121 ! ------------- 1122 1123 REAL sig(llm + 1), dsig(llm) 1124 1125 INTEGER l 1126 REAL snorm 1127 REAL alpha, beta, gama, delta, deltaz, h 1128 INTEGER np, ierr 1129 REAL pi, x 1130 1131 !----------------------------------------------------------------------- 1132 1133 pi = 2. * ASIN(1.) 1134 1135 OPEN(99, file = 'sigma.def', status = 'old', form = 'formatted', & 1136 & iostat = ierr) 1137 1138 !----------------------------------------------------------------------- 1139 ! cas 1 on lit les options dans sigma.def: 1140 ! ---------------------------------------- 1141 1142 IF (ierr==0) THEN 1143 1144 PRINT*, 'WARNING!!! on lit les options dans sigma.def' 1145 READ(99, *) deltaz 1146 READ(99, *) h 1147 READ(99, *) beta 1148 READ(99, *) gama 1149 READ(99, *) delta 1150 READ(99, *) np 1151 CLOSE(99) 1152 alpha = deltaz / (llm * h) 1153 1154 DO l = 1, llm 1155 dsig(l) = (alpha + (1. - alpha) * exp(-beta * (llm - l))) * & 1156 & ((tanh(gama * l) / tanh(gama * llm))**np + & 1157 & (1. - l / FLOAT(llm)) * delta) 1158 END DO 1159 1160 sig(1) = 1. 1161 DO l = 1, llm - 1 1162 sig(l + 1) = sig(l) * (1. - dsig(l)) / (1. + dsig(l)) 1163 END DO 1164 sig(llm + 1) = 0. 1165 1166 DO l = 1, llm 1167 dsig(l) = sig(l) - sig(l + 1) 1168 END DO 1169 1170 ELSE 1171 !----------------------------------------------------------------------- 1172 ! cas 2 ancienne discretisation (LMD5...): 1173 ! ---------------------------------------- 1174 1175 PRINT*, 'WARNING!!! Ancienne discretisation verticale' 1176 1177 h = 7. 1178 snorm = 0. 1179 DO l = 1, llm 1180 x = 2. * asin(1.) * (FLOAT(l) - 0.5) / float(llm + 1) 1181 dsig(l) = 1.0 + 7.0 * SIN(x)**2 1182 snorm = snorm + dsig(l) 1183 ENDDO 1184 snorm = 1. / snorm 1185 DO l = 1, llm 1186 dsig(l) = dsig(l) * snorm 1187 ENDDO 1188 sig(llm + 1) = 0. 1189 DO l = llm, 1, -1 1190 sig(l) = sig(l + 1) + dsig(l) 1191 ENDDO 1192 1193 ENDIF 1194 1195 DO l = 1, llm 1196 nivsigs(l) = FLOAT(l) 975 1197 ENDDO 976 ENDDO 977 978 RETURN 979 END 980 981 982 SUBROUTINE abort_gcm(modname, message, ierr) 983 984 USE IOIPSL 985 986 ! Stops the simulation cleanly, closing files and printing various 987 ! comments 988 989 ! Input: modname = name of calling program 990 ! message = stuff to print 991 ! ierr = severity of situation ( = 0 normal ) 992 993 character(len = *) modname 994 integer ierr 995 character(len = *) message 996 997 write(*, *) 'in abort_gcm' 998 CALL histclo 999 ! CALL histclo(2) 1000 ! CALL histclo(3) 1001 ! CALL histclo(4) 1002 ! CALL histclo(5) 1003 write(*, *) 'out of histclo' 1004 write(*, *) 'Stopping in ', modname 1005 write(*, *) 'Reason = ', message 1006 CALL getin_dump 1007 1008 if (ierr == 0) then 1009 write(*, *) 'Everything is cool' 1010 else 1011 write(*, *) 'Houston, we have a problem ', ierr 1012 endif 1013 STOP 1014 END 1015 REAL FUNCTION fq_sat(kelvin, millibar) 1016 1017 IMPLICIT none 1018 !====================================================================== 1019 ! Autheur(s): Z.X. Li (LMD/CNRS) 1020 ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.) 1021 !====================================================================== 1022 ! Arguments: 1023 ! kelvin---input-R: temperature en Kelvin 1024 ! millibar--input-R: pression en mb 1025 1026 ! fq_sat----output-R: vapeur d'eau saturante en kg/kg 1027 !====================================================================== 1028 1029 REAL kelvin, millibar 1030 1031 REAL r2es 1032 PARAMETER (r2es = 611.14 * 18.0153 / 28.9644) 1033 1034 REAL r3les, r3ies, r3es 1035 PARAMETER (R3LES = 17.269) 1036 PARAMETER (R3IES = 21.875) 1037 1038 REAL r4les, r4ies, r4es 1039 PARAMETER (R4LES = 35.86) 1040 PARAMETER (R4IES = 7.66) 1041 1042 REAL rtt 1043 PARAMETER (rtt = 273.16) 1044 1045 REAL retv 1046 PARAMETER (retv = 28.9644 / 18.0153 - 1.0) 1047 1048 REAL zqsat 1049 REAL temp, pres 1050 ! ------------------------------------------------------------------ 1051 1052 temp = kelvin 1053 pres = millibar * 100.0 1054 ! write(*,*)'kelvin,millibar=',kelvin,millibar 1055 ! write(*,*)'temp,pres=',temp,pres 1056 1057 IF (temp <= rtt) THEN 1058 r3es = r3ies 1059 r4es = r4ies 1060 ELSE 1061 r3es = r3les 1062 r4es = r4les 1063 ENDIF 1064 1065 zqsat = r2es / pres * EXP (r3es * (temp - rtt) / (temp - r4es)) 1066 zqsat = MIN(0.5, ZQSAT) 1067 zqsat = zqsat / (1. - retv * zqsat) 1068 1069 fq_sat = zqsat 1070 1071 RETURN 1072 END 1073 1074 SUBROUTINE gr_dyn_fi(nfield, im, jm, ngrid, pdyn, pfi) 1075 IMPLICIT NONE 1076 !======================================================================= 1077 ! passage d'un champ de la grille scalaire a la grille physique 1078 !======================================================================= 1079 1080 !----------------------------------------------------------------------- 1081 ! declarations: 1082 ! ------------- 1083 1084 INTEGER im, jm, ngrid, nfield 1085 REAL pdyn(im, jm, nfield) 1086 REAL pfi(ngrid, nfield) 1087 1088 INTEGER j, ifield, ig 1089 1090 !----------------------------------------------------------------------- 1091 ! calcul: 1092 ! ------- 1093 1094 IF(ngrid/=2 + (jm - 2) * (im - 1).AND.ngrid/=1) & 1095 & STOP 'probleme de dim' 1096 ! traitement des poles 1097 CALL SCOPY(nfield, pdyn, im * jm, pfi, ngrid) 1098 CALL SCOPY(nfield, pdyn(1, jm, 1), im * jm, pfi(ngrid, 1), ngrid) 1099 1100 ! traitement des point normaux 1101 DO ifield = 1, nfield 1102 DO j = 2, jm - 1 1103 ig = 2 + (j - 2) * (im - 1) 1104 CALL SCOPY(im - 1, pdyn(1, j, ifield), 1, pfi(ig, ifield), 1) 1198 1199 DO l = 1, llmp1 1200 nivsig(l) = FLOAT(l) 1105 1201 ENDDO 1106 ENDDO 1107 1108 RETURN 1109 END 1110 1111 SUBROUTINE disvert0(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) 1112 1113 ! Ancienne version disvert dont on a modifie nom pour utiliser 1114 ! le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes) 1115 ! (MPL 18092012) 1116 1117 ! Auteur : P. Le Van . 1118 1119 IMPLICIT NONE 1120 1121 include "dimensions.h" 1122 include "paramet.h" 1123 1124 !======================================================================= 1125 1126 1127 ! s = sigma ** kappa : coordonnee verticale 1128 ! dsig(l) : epaisseur de la couche l ds la coord. s 1129 ! sig(l) : sigma a l'interface des couches l et l-1 1130 ! ds(l) : distance entre les couches l et l-1 en coord.s 1131 1132 !======================================================================= 1133 1134 REAL pa, preff 1135 REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) 1136 REAL presnivs(llm) 1137 1138 ! declarations: 1139 ! ------------- 1140 1141 REAL sig(llm + 1), dsig(llm) 1142 1143 INTEGER l 1144 REAL snorm 1145 REAL alpha, beta, gama, delta, deltaz, h 1146 INTEGER np, ierr 1147 REAL pi, x 1148 1149 !----------------------------------------------------------------------- 1150 1151 pi = 2. * ASIN(1.) 1152 1153 OPEN(99, file = 'sigma.def', status = 'old', form = 'formatted', & 1154 & iostat = ierr) 1155 1156 !----------------------------------------------------------------------- 1157 ! cas 1 on lit les options dans sigma.def: 1158 ! ---------------------------------------- 1159 1160 IF (ierr==0) THEN 1161 1162 PRINT*, 'WARNING!!! on lit les options dans sigma.def' 1163 READ(99, *) deltaz 1164 READ(99, *) h 1165 READ(99, *) beta 1166 READ(99, *) gama 1167 READ(99, *) delta 1168 READ(99, *) np 1169 CLOSE(99) 1170 alpha = deltaz / (llm * h) 1202 1203 ! .... Calculs de ap(l) et de bp(l) .... 1204 ! ......................................... 1205 1206 ! ..... pa et preff sont lus sur les fichiers start par lectba ..... 1207 1208 bp(llmp1) = 0. 1171 1209 1172 1210 DO l = 1, llm 1173 dsig(l) = (alpha + (1. - alpha) * exp(-beta * (llm - l))) * & 1174 & ((tanh(gama * l) / tanh(gama * llm))**np + & 1175 & (1. - l / FLOAT(llm)) * delta) 1176 END DO 1177 1178 sig(1) = 1. 1179 DO l = 1, llm - 1 1180 sig(l + 1) = sig(l) * (1. - dsig(l)) / (1. + dsig(l)) 1181 END DO 1182 sig(llm + 1) = 0. 1211 !c 1212 !cc ap(l) = 0. 1213 !cc bp(l) = sig(l) 1214 1215 bp(l) = EXP(1. - 1. / (sig(l) * sig(l))) 1216 ap(l) = pa * (sig(l) - bp(l)) 1217 1218 ENDDO 1219 ap(llmp1) = pa * (sig(llmp1) - bp(llmp1)) 1220 1221 PRINT *, ' BP ' 1222 PRINT *, bp 1223 PRINT *, ' AP ' 1224 PRINT *, ap 1183 1225 1184 1226 DO l = 1, llm 1185 dsig(l) = sig(l) - sig(l + 1) 1186 END DO 1187 1188 ELSE 1189 !----------------------------------------------------------------------- 1190 ! cas 2 ancienne discretisation (LMD5...): 1191 ! ---------------------------------------- 1192 1193 PRINT*, 'WARNING!!! Ancienne discretisation verticale' 1194 1195 h = 7. 1196 snorm = 0. 1197 DO l = 1, llm 1198 x = 2. * asin(1.) * (FLOAT(l) - 0.5) / float(llm + 1) 1199 dsig(l) = 1.0 + 7.0 * SIN(x)**2 1200 snorm = snorm + dsig(l) 1227 dpres(l) = bp(l) - bp(l + 1) 1228 presnivs(l) = 0.5 * (ap(l) + bp(l) * preff + ap(l + 1) + bp(l + 1) * preff) 1201 1229 ENDDO 1202 snorm = 1. / snorm 1203 DO l = 1, llm 1204 dsig(l) = dsig(l) * snorm 1205 ENDDO 1206 sig(llm + 1) = 0. 1207 DO l = llm, 1, -1 1208 sig(l) = sig(l + 1) + dsig(l) 1209 ENDDO 1210 1211 ENDIF 1212 1213 DO l = 1, llm 1214 nivsigs(l) = FLOAT(l) 1215 ENDDO 1216 1217 DO l = 1, llmp1 1218 nivsig(l) = FLOAT(l) 1219 ENDDO 1220 1221 ! .... Calculs de ap(l) et de bp(l) .... 1222 ! ......................................... 1223 1224 ! ..... pa et preff sont lus sur les fichiers start par lectba ..... 1225 1226 bp(llmp1) = 0. 1227 1228 DO l = 1, llm 1229 !c 1230 !cc ap(l) = 0. 1231 !cc bp(l) = sig(l) 1232 1233 bp(l) = EXP(1. - 1. / (sig(l) * sig(l))) 1234 ap(l) = pa * (sig(l) - bp(l)) 1235 1236 ENDDO 1237 ap(llmp1) = pa * (sig(llmp1) - bp(llmp1)) 1238 1239 PRINT *, ' BP ' 1240 PRINT *, bp 1241 PRINT *, ' AP ' 1242 PRINT *, ap 1243 1244 DO l = 1, llm 1245 dpres(l) = bp(l) - bp(l + 1) 1246 presnivs(l) = 0.5 * (ap(l) + bp(l) * preff + ap(l + 1) + bp(l + 1) * preff) 1247 ENDDO 1248 1249 PRINT *, ' PRESNIVS ' 1250 PRINT *, presnivs 1251 1252 RETURN 1253 END 1254 1255 !!====================================================================== 1256 ! SUBROUTINE read_tsurf1d(knon,sst_out) 1257 1258 !! This subroutine specifies the surface temperature to be used in 1D simulations 1259 1260 ! USE dimphy, ONLY: klon 1261 1262 ! INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid 1263 ! REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model 1264 1265 ! INTEGER :: i 1266 !! COMMON defined in lmdz1d.F: 1267 ! real ts_cur 1268 ! common /sst_forcing/ts_cur 1269 1270 ! DO i = 1, knon 1271 ! sst_out(i) = ts_cur 1272 ! ENDDO 1273 1274 ! END SUBROUTINE read_tsurf1d 1275 1276 !=============================================================== 1277 subroutine advect_vert(llm, w, dt, q, plev) 1278 !=============================================================== 1279 ! Schema amont pour l'advection verticale en 1D 1280 ! w est la vitesse verticale dp/dt en Pa/s 1281 ! Traitement en volumes finis 1282 ! d / dt ( zm q ) = delta_z ( omega q ) 1283 ! d / dt ( zm ) = delta_z ( omega ) 1284 ! avec zm = delta_z ( p ) 1285 ! si * designe la valeur au pas de temps t+dt 1286 ! zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l) 1287 ! zm*(l) -zm(l) = w(l+1) - w(l) 1288 ! avec w=omega * dt 1289 !--------------------------------------------------------------- 1290 implicit none 1291 ! arguments 1292 integer llm 1293 real w(llm + 1), q(llm), plev(llm + 1), dt 1294 1295 ! local 1296 integer l 1297 real zwq(llm + 1), zm(llm + 1), zw(llm + 1) 1298 real qold 1299 1300 !--------------------------------------------------------------- 1301 1302 do l = 1, llm 1303 zw(l) = dt * w(l) 1304 zm(l) = plev(l) - plev(l + 1) 1305 zwq(l) = q(l) * zw(l) 1306 enddo 1307 zwq(llm + 1) = 0. 1308 zw(llm + 1) = 0. 1309 1310 do l = 1, llm 1311 qold = q(l) 1312 q(l) = (q(l) * zm(l) + zwq(l + 1) - zwq(l)) / (zm(l) + zw(l + 1) - zw(l)) 1313 PRINT*, 'ADV Q ', zm(l), zw(l), zwq(l), qold, q(l) 1314 enddo 1315 1316 return 1317 end 1318 1319 !=============================================================== 1320 1321 1322 SUBROUTINE advect_va(llm, omega, d_t_va, d_q_va, d_u_va, d_v_va, & 1323 & q, temp, u, v, play) 1324 !itlmd 1325 !---------------------------------------------------------------------- 1326 ! Calcul de l'advection verticale (ascendance et subsidence) de 1327 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur 1328 ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 1329 ! sans WTG rajouter une advection horizontale 1330 !---------------------------------------------------------------------- 1331 implicit none 1332 include "YOMCST.h" 1333 ! argument 1334 integer llm 1335 real omega(llm + 1), d_t_va(llm), d_q_va(llm, 3) 1336 real d_u_va(llm), d_v_va(llm) 1337 real q(llm, 3), temp(llm) 1338 real u(llm), v(llm) 1339 real play(llm) 1340 ! interne 1341 integer l 1342 real alpha, omgdown, omgup 1343 1344 do l = 1, llm 1345 if(l==1) then 1346 !si omgup pour la couche 1, alors tendance nulle 1347 omgdown = max(omega(2), 0.0) 1348 alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1))) 1349 d_t_va(l) = alpha * (omgdown) - omgdown * (temp(l) - temp(l + 1)) & 1350 & / (play(l) - play(l + 1)) 1351 1352 d_q_va(l, :) = -omgdown * (q(l, :) - q(l + 1, :)) / (play(l) - play(l + 1)) 1353 1354 d_u_va(l) = -omgdown * (u(l) - u(l + 1)) / (play(l) - play(l + 1)) 1355 d_v_va(l) = -omgdown * (v(l) - v(l + 1)) / (play(l) - play(l + 1)) 1356 1357 elseif(l==llm) then 1358 omgup = min(omega(l), 0.0) 1359 alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1))) 1360 d_t_va(l) = alpha * (omgup) - & 1361 1362 !bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l)) 1363 & omgup * (temp(l - 1) - temp(l)) / (play(l - 1) - play(l)) 1364 d_q_va(l, :) = -omgup * (q(l - 1, :) - q(l, :)) / (play(l - 1) - play(l)) 1365 d_u_va(l) = -omgup * (u(l - 1) - u(l)) / (play(l - 1) - play(l)) 1366 d_v_va(l) = -omgup * (v(l - 1) - v(l)) / (play(l - 1) - play(l)) 1367 1368 else 1369 omgup = min(omega(l), 0.0) 1370 omgdown = max(omega(l + 1), 0.0) 1371 alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1))) 1372 d_t_va(l) = alpha * (omgup + omgdown) - omgdown * (temp(l) - temp(l + 1)) & 1373 & / (play(l) - play(l + 1)) - & 1374 !bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l)) 1375 & omgup * (temp(l - 1) - temp(l)) / (play(l - 1) - play(l)) 1376 ! PRINT*, ' ??? ' 1377 1378 d_q_va(l, :) = -omgdown * (q(l, :) - q(l + 1, :)) & 1379 & / (play(l) - play(l + 1)) - & 1380 & omgup * (q(l - 1, :) - q(l, :)) / (play(l - 1) - play(l)) 1381 d_u_va(l) = -omgdown * (u(l) - u(l + 1)) & 1382 & / (play(l) - play(l + 1)) - & 1383 & omgup * (u(l - 1) - u(l)) / (play(l - 1) - play(l)) 1384 d_v_va(l) = -omgdown * (v(l) - v(l + 1)) & 1385 & / (play(l) - play(l + 1)) - & 1386 & omgup * (v(l - 1) - v(l)) / (play(l - 1) - play(l)) 1387 1388 endif 1389 1390 enddo 1391 !fin itlmd 1392 return 1393 end 1394 ! SUBROUTINE lstendH(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va, 1395 SUBROUTINE lstendH(llm, nqtot, omega, d_t_va, d_q_va, & 1396 & q, temp, u, v, play) 1397 !itlmd 1398 !---------------------------------------------------------------------- 1399 ! Calcul de l'advection verticale (ascendance et subsidence) de 1400 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur 1401 ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 1402 ! sans WTG rajouter une advection horizontale 1403 !---------------------------------------------------------------------- 1404 implicit none 1405 include "YOMCST.h" 1406 ! argument 1407 integer llm, nqtot 1408 real omega(llm + 1), d_t_va(llm), d_q_va(llm, nqtot) 1409 ! real d_u_va(llm), d_v_va(llm) 1410 real q(llm, nqtot), temp(llm) 1411 real u(llm), v(llm) 1412 real play(llm) 1413 real cor(llm) 1414 ! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm) 1415 real dph(llm), dqdp(llm), dtdp(llm) 1416 ! interne 1417 integer k 1418 real omdn, omup 1419 1420 ! dudp=0. 1421 ! dvdp=0. 1422 dqdp = 0. 1423 dtdp = 0. 1424 ! d_u_va=0. 1425 ! d_v_va=0. 1426 1427 cor(:) = rkappa * temp * (1. + q(:, 1) * rv / rd) / (play * (1. + q(:, 1))) 1428 1429 do k = 2, llm - 1 1430 1431 dph (k - 1) = (play(k) - play(k - 1)) 1432 ! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1) 1433 ! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1) 1434 dqdp (k - 1) = (q (k, 1) - q (k - 1, 1)) / dph(k - 1) 1435 dtdp (k - 1) = (temp(k) - temp(k - 1)) / dph(k - 1) 1436 1437 enddo 1438 1439 ! dudp ( llm ) = dudp ( llm-1 ) 1440 ! dvdp ( llm ) = dvdp ( llm-1 ) 1441 dqdp (llm) = dqdp (llm - 1) 1442 dtdp (llm) = dtdp (llm - 1) 1443 1444 do k = 2, llm - 1 1445 omdn = max(0.0, omega(k + 1)) 1446 omup = min(0.0, omega(k)) 1447 1448 ! d_u_va(k) = -omdn*dudp(k)-omup*dudp(k-1) 1449 ! d_v_va(k) = -omdn*dvdp(k)-omup*dvdp(k-1) 1450 d_q_va(k, 1) = -omdn * dqdp(k) - omup * dqdp(k - 1) 1451 d_t_va(k) = -omdn * dtdp(k) - omup * dtdp(k - 1) + (omup + omdn) * cor(k) 1452 enddo 1453 1454 omdn = max(0.0, omega(2)) 1455 omup = min(0.0, omega(llm)) 1456 ! d_u_va( 1 ) = -omdn*dudp( 1 ) 1457 ! d_u_va(llm) = -omup*dudp(llm) 1458 ! d_v_va( 1 ) = -omdn*dvdp( 1 ) 1459 ! d_v_va(llm) = -omup*dvdp(llm) 1460 d_q_va(1, 1) = -omdn * dqdp(1) 1461 d_q_va(llm, 1) = -omup * dqdp(llm) 1462 d_t_va(1) = -omdn * dtdp(1) + omdn * cor(1) 1463 d_t_va(llm) = -omup * dtdp(llm)!+omup*cor(llm) 1464 1465 ! if(abs(rlat(1))>10.) then 1466 ! Calculate the tendency due agestrophic motions 1467 ! du_age = fcoriolis*(v-vg) 1468 ! dv_age = fcoriolis*(ug-u) 1469 ! endif 1470 1471 ! CALL writefield_phy('d_t_va',d_t_va,llm) 1472 1473 return 1474 end 1475 1476 !====================================================================== 1477 1478 ! Subroutines for nudging 1479 1480 Subroutine Nudge_RHT_init (paprs, pplay, t, q, t_targ, rh_targ) 1481 ! ======================================================== 1482 USE dimphy 1483 1484 implicit none 1485 1486 ! ======================================================== 1487 REAL paprs(klon, klevp1) 1488 REAL pplay(klon, klev) 1489 1490 ! Variables d'etat 1491 REAL t(klon, klev) 1492 REAL q(klon, klev) 1493 1494 ! Profiles cible 1495 REAL t_targ(klon, klev) 1496 REAL rh_targ(klon, klev) 1497 1498 INTEGER k, i 1499 REAL zx_qs 1500 1501 ! Declaration des constantes et des fonctions thermodynamiques 1502 1503 include "YOMCST.h" 1504 include "YOETHF.h" 1505 1506 ! ---------------------------------------- 1507 ! Statement functions 1508 include "FCTTRE.h" 1509 ! ---------------------------------------- 1510 1511 DO k = 1, klev 1512 DO i = 1, klon 1513 t_targ(i, k) = t(i, k) 1514 IF (t(i, k)<RTT) THEN 1515 zx_qs = qsats(t(i, k)) / (pplay(i, k)) 1516 ELSE 1517 zx_qs = qsatl(t(i, k)) / (pplay(i, k)) 1518 ENDIF 1519 rh_targ(i, k) = q(i, k) / zx_qs 1520 ENDDO 1521 ENDDO 1522 print *, 't_targ', t_targ 1523 print *, 'rh_targ', rh_targ 1524 1525 RETURN 1526 END 1527 1528 Subroutine Nudge_UV_init (paprs, pplay, u, v, u_targ, v_targ) 1529 ! ======================================================== 1530 USE dimphy 1531 1532 implicit none 1533 1534 ! ======================================================== 1535 REAL paprs(klon, klevp1) 1536 REAL pplay(klon, klev) 1537 1538 ! Variables d'etat 1539 REAL u(klon, klev) 1540 REAL v(klon, klev) 1541 1542 ! Profiles cible 1543 REAL u_targ(klon, klev) 1544 REAL v_targ(klon, klev) 1545 1546 INTEGER k, i 1547 1548 DO k = 1, klev 1549 DO i = 1, klon 1550 u_targ(i, k) = u(i, k) 1551 v_targ(i, k) = v(i, k) 1552 ENDDO 1553 ENDDO 1554 print *, 'u_targ', u_targ 1555 print *, 'v_targ', v_targ 1556 1557 RETURN 1558 END 1559 1560 Subroutine Nudge_RHT (dtime, paprs, pplay, t_targ, rh_targ, t, q, & 1561 & d_t, d_q) 1562 ! ======================================================== 1563 USE dimphy 1564 1565 implicit none 1566 1567 ! ======================================================== 1568 REAL dtime 1569 REAL paprs(klon, klevp1) 1570 REAL pplay(klon, klev) 1571 1572 ! Variables d'etat 1573 REAL t(klon, klev) 1574 REAL q(klon, klev) 1575 1576 ! Tendances 1577 REAL d_t(klon, klev) 1578 REAL d_q(klon, klev) 1579 1580 ! Profiles cible 1581 REAL t_targ(klon, klev) 1582 REAL rh_targ(klon, klev) 1583 1584 ! Temps de relaxation 1585 REAL tau 1586 !c DATA tau /3600./ 1587 !! DATA tau /5400./ 1588 DATA tau /1800./ 1589 1590 INTEGER k, i 1591 REAL zx_qs, rh, tnew, d_rh, rhnew 1592 1593 ! Declaration des constantes et des fonctions thermodynamiques 1594 1595 include "YOMCST.h" 1596 include "YOETHF.h" 1597 1598 ! ---------------------------------------- 1599 ! Statement functions 1600 include "FCTTRE.h" 1601 ! ---------------------------------------- 1602 1603 print *, 'dtime, tau ', dtime, tau 1604 print *, 't_targ', t_targ 1605 print *, 'rh_targ', rh_targ 1606 print *, 'temp ', t 1607 print *, 'hum ', q 1608 1609 DO k = 1, klev 1610 DO i = 1, klon 1611 IF (paprs(i, 1) - pplay(i, k) > 10000.) THEN 1230 1231 PRINT *, ' PRESNIVS ' 1232 PRINT *, presnivs 1233 END SUBROUTINE disvert0 1234 1235 subroutine advect_vert(llm, w, dt, q, plev) 1236 !=============================================================== 1237 ! Schema amont pour l'advection verticale en 1D 1238 ! w est la vitesse verticale dp/dt en Pa/s 1239 ! Traitement en volumes finis 1240 ! d / dt ( zm q ) = delta_z ( omega q ) 1241 ! d / dt ( zm ) = delta_z ( omega ) 1242 ! avec zm = delta_z ( p ) 1243 ! si * designe la valeur au pas de temps t+dt 1244 ! zm*(l) q*(l) - zm(l) q(l) = w(l+1) q(l+1) - w(l) q(l) 1245 ! zm*(l) -zm(l) = w(l+1) - w(l) 1246 ! avec w=omega * dt 1247 !--------------------------------------------------------------- 1248 implicit none 1249 ! arguments 1250 integer llm 1251 real w(llm + 1), q(llm), plev(llm + 1), dt 1252 1253 ! local 1254 integer l 1255 real zwq(llm + 1), zm(llm + 1), zw(llm + 1) 1256 real qold 1257 1258 !--------------------------------------------------------------- 1259 1260 do l = 1, llm 1261 zw(l) = dt * w(l) 1262 zm(l) = plev(l) - plev(l + 1) 1263 zwq(l) = q(l) * zw(l) 1264 enddo 1265 zwq(llm + 1) = 0. 1266 zw(llm + 1) = 0. 1267 1268 do l = 1, llm 1269 qold = q(l) 1270 q(l) = (q(l) * zm(l) + zwq(l + 1) - zwq(l)) / (zm(l) + zw(l + 1) - zw(l)) 1271 PRINT*, 'ADV Q ', zm(l), zw(l), zwq(l), qold, q(l) 1272 enddo 1273 end SUBROUTINE advect_vert 1274 1275 SUBROUTINE advect_va(llm, omega, d_t_va, d_q_va, d_u_va, d_v_va, & 1276 & q, temp, u, v, play) 1277 !itlmd 1278 !---------------------------------------------------------------------- 1279 ! Calcul de l'advection verticale (ascendance et subsidence) de 1280 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur 1281 ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 1282 ! sans WTG rajouter une advection horizontale 1283 !---------------------------------------------------------------------- 1284 implicit none 1285 include "YOMCST.h" 1286 ! argument 1287 integer llm 1288 real omega(llm + 1), d_t_va(llm), d_q_va(llm, 3) 1289 real d_u_va(llm), d_v_va(llm) 1290 real q(llm, 3), temp(llm) 1291 real u(llm), v(llm) 1292 real play(llm) 1293 ! interne 1294 integer l 1295 real alpha, omgdown, omgup 1296 1297 do l = 1, llm 1298 if(l==1) then 1299 !si omgup pour la couche 1, alors tendance nulle 1300 omgdown = max(omega(2), 0.0) 1301 alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1))) 1302 d_t_va(l) = alpha * (omgdown) - omgdown * (temp(l) - temp(l + 1)) & 1303 & / (play(l) - play(l + 1)) 1304 1305 d_q_va(l, :) = -omgdown * (q(l, :) - q(l + 1, :)) / (play(l) - play(l + 1)) 1306 1307 d_u_va(l) = -omgdown * (u(l) - u(l + 1)) / (play(l) - play(l + 1)) 1308 d_v_va(l) = -omgdown * (v(l) - v(l + 1)) / (play(l) - play(l + 1)) 1309 1310 elseif(l==llm) then 1311 omgup = min(omega(l), 0.0) 1312 alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1))) 1313 d_t_va(l) = alpha * (omgup) - & 1314 1315 !bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l)) 1316 & omgup * (temp(l - 1) - temp(l)) / (play(l - 1) - play(l)) 1317 d_q_va(l, :) = -omgup * (q(l - 1, :) - q(l, :)) / (play(l - 1) - play(l)) 1318 d_u_va(l) = -omgup * (u(l - 1) - u(l)) / (play(l - 1) - play(l)) 1319 d_v_va(l) = -omgup * (v(l - 1) - v(l)) / (play(l - 1) - play(l)) 1320 1321 else 1322 omgup = min(omega(l), 0.0) 1323 omgdown = max(omega(l + 1), 0.0) 1324 alpha = rkappa * temp(l) * (1. + q(l, 1) * rv / rd) / (play(l) * (1. + q(l, 1))) 1325 d_t_va(l) = alpha * (omgup + omgdown) - omgdown * (temp(l) - temp(l + 1)) & 1326 & / (play(l) - play(l + 1)) - & 1327 !bug? & omgup*(temp(l-1)-temp(l))/(play(l-1)-plev(l)) 1328 & omgup * (temp(l - 1) - temp(l)) / (play(l - 1) - play(l)) 1329 ! PRINT*, ' ??? ' 1330 1331 d_q_va(l, :) = -omgdown * (q(l, :) - q(l + 1, :)) & 1332 & / (play(l) - play(l + 1)) - & 1333 & omgup * (q(l - 1, :) - q(l, :)) / (play(l - 1) - play(l)) 1334 d_u_va(l) = -omgdown * (u(l) - u(l + 1)) & 1335 & / (play(l) - play(l + 1)) - & 1336 & omgup * (u(l - 1) - u(l)) / (play(l - 1) - play(l)) 1337 d_v_va(l) = -omgdown * (v(l) - v(l + 1)) & 1338 & / (play(l) - play(l + 1)) - & 1339 & omgup * (v(l - 1) - v(l)) / (play(l - 1) - play(l)) 1340 1341 endif 1342 1343 enddo 1344 !fin itlmd 1345 end SUBROUTINE advect_va 1346 1347 1348 SUBROUTINE lstendH(llm, nqtot, omega, d_t_va, d_q_va, q, temp, u, v, play) 1349 !itlmd 1350 !---------------------------------------------------------------------- 1351 ! Calcul de l'advection verticale (ascendance et subsidence) de 1352 ! temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur 1353 ! a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 1354 ! sans WTG rajouter une advection horizontale 1355 !---------------------------------------------------------------------- 1356 implicit none 1357 include "YOMCST.h" 1358 ! argument 1359 integer llm, nqtot 1360 real omega(llm + 1), d_t_va(llm), d_q_va(llm, nqtot) 1361 ! real d_u_va(llm), d_v_va(llm) 1362 real q(llm, nqtot), temp(llm) 1363 real u(llm), v(llm) 1364 real play(llm) 1365 real cor(llm) 1366 ! real dph(llm),dudp(llm),dvdp(llm),dqdp(llm),dtdp(llm) 1367 real dph(llm), dqdp(llm), dtdp(llm) 1368 ! interne 1369 integer k 1370 real omdn, omup 1371 1372 ! dudp=0. 1373 ! dvdp=0. 1374 dqdp = 0. 1375 dtdp = 0. 1376 ! d_u_va=0. 1377 ! d_v_va=0. 1378 1379 cor(:) = rkappa * temp * (1. + q(:, 1) * rv / rd) / (play * (1. + q(:, 1))) 1380 1381 do k = 2, llm - 1 1382 1383 dph (k - 1) = (play(k) - play(k - 1)) 1384 ! dudp (k-1) = (u (k )- u (k-1 ))/dph(k-1) 1385 ! dvdp (k-1) = (v (k )- v (k-1 ))/dph(k-1) 1386 dqdp (k - 1) = (q (k, 1) - q (k - 1, 1)) / dph(k - 1) 1387 dtdp (k - 1) = (temp(k) - temp(k - 1)) / dph(k - 1) 1388 1389 enddo 1390 1391 ! dudp ( llm ) = dudp ( llm-1 ) 1392 ! dvdp ( llm ) = dvdp ( llm-1 ) 1393 dqdp (llm) = dqdp (llm - 1) 1394 dtdp (llm) = dtdp (llm - 1) 1395 1396 do k = 2, llm - 1 1397 omdn = max(0.0, omega(k + 1)) 1398 omup = min(0.0, omega(k)) 1399 1400 ! d_u_va(k) = -omdn*dudp(k)-omup*dudp(k-1) 1401 ! d_v_va(k) = -omdn*dvdp(k)-omup*dvdp(k-1) 1402 d_q_va(k, 1) = -omdn * dqdp(k) - omup * dqdp(k - 1) 1403 d_t_va(k) = -omdn * dtdp(k) - omup * dtdp(k - 1) + (omup + omdn) * cor(k) 1404 enddo 1405 1406 omdn = max(0.0, omega(2)) 1407 omup = min(0.0, omega(llm)) 1408 ! d_u_va( 1 ) = -omdn*dudp( 1 ) 1409 ! d_u_va(llm) = -omup*dudp(llm) 1410 ! d_v_va( 1 ) = -omdn*dvdp( 1 ) 1411 ! d_v_va(llm) = -omup*dvdp(llm) 1412 d_q_va(1, 1) = -omdn * dqdp(1) 1413 d_q_va(llm, 1) = -omup * dqdp(llm) 1414 d_t_va(1) = -omdn * dtdp(1) + omdn * cor(1) 1415 d_t_va(llm) = -omup * dtdp(llm)!+omup*cor(llm) 1416 1417 ! if(abs(rlat(1))>10.) then 1418 ! Calculate the tendency due agestrophic motions 1419 ! du_age = fcoriolis*(v-vg) 1420 ! dv_age = fcoriolis*(ug-u) 1421 ! endif 1422 1423 ! CALL writefield_phy('d_t_va',d_t_va,llm) 1424 end SUBROUTINE lstendH 1425 1426 1427 Subroutine Nudge_RHT_init (paprs, pplay, t, q, t_targ, rh_targ) 1428 ! ======================================================== 1429 USE dimphy 1430 1431 implicit none 1432 1433 ! ======================================================== 1434 REAL paprs(klon, klevp1) 1435 REAL pplay(klon, klev) 1436 1437 ! Variables d'etat 1438 REAL t(klon, klev) 1439 REAL q(klon, klev) 1440 1441 ! Profiles cible 1442 REAL t_targ(klon, klev) 1443 REAL rh_targ(klon, klev) 1444 1445 INTEGER k, i 1446 REAL zx_qs 1447 1448 ! Declaration des constantes et des fonctions thermodynamiques 1449 1450 include "YOMCST.h" 1451 include "YOETHF.h" 1452 1453 ! ---------------------------------------- 1454 ! Statement functions 1455 include "FCTTRE.h" 1456 ! ---------------------------------------- 1457 1458 DO k = 1, klev 1459 DO i = 1, klon 1460 t_targ(i, k) = t(i, k) 1612 1461 IF (t(i, k)<RTT) THEN 1613 1462 zx_qs = qsats(t(i, k)) / (pplay(i, k)) … … 1615 1464 zx_qs = qsatl(t(i, k)) / (pplay(i, k)) 1616 1465 ENDIF 1617 rh = q(i, k) / zx_qs 1618 1619 d_t(i, k) = d_t(i, k) + 1. / tau * (t_targ(i, k) - t(i, k)) 1620 d_rh = 1. / tau * (rh_targ(i, k) - rh) 1621 1622 tnew = t(i, k) + d_t(i, k) * dtime 1623 !jyg< 1624 ! Formule pour q : 1625 ! d_q = (1/tau) [rh_targ*qsat(T_new) - q] 1626 1627 ! Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new) 1628 ! qui n'etait pas correcte. 1629 1630 IF (tnew<RTT) THEN 1631 zx_qs = qsats(tnew) / (pplay(i, k)) 1632 ELSE 1633 zx_qs = qsatl(tnew) / (pplay(i, k)) 1466 rh_targ(i, k) = q(i, k) / zx_qs 1467 ENDDO 1468 ENDDO 1469 print *, 't_targ', t_targ 1470 print *, 'rh_targ', rh_targ 1471 1472 RETURN 1473 END SUBROUTINE nudge_rht_init 1474 1475 Subroutine Nudge_UV_init (paprs, pplay, u, v, u_targ, v_targ) 1476 ! ======================================================== 1477 USE dimphy 1478 1479 implicit none 1480 1481 ! ======================================================== 1482 REAL paprs(klon, klevp1) 1483 REAL pplay(klon, klev) 1484 1485 ! Variables d'etat 1486 REAL u(klon, klev) 1487 REAL v(klon, klev) 1488 1489 ! Profiles cible 1490 REAL u_targ(klon, klev) 1491 REAL v_targ(klon, klev) 1492 1493 INTEGER k, i 1494 1495 DO k = 1, klev 1496 DO i = 1, klon 1497 u_targ(i, k) = u(i, k) 1498 v_targ(i, k) = v(i, k) 1499 ENDDO 1500 ENDDO 1501 print *, 'u_targ', u_targ 1502 print *, 'v_targ', v_targ 1503 1504 RETURN 1505 END 1506 1507 Subroutine Nudge_RHT (dtime, paprs, pplay, t_targ, rh_targ, t, q, & 1508 & d_t, d_q) 1509 ! ======================================================== 1510 USE dimphy 1511 1512 implicit none 1513 1514 ! ======================================================== 1515 REAL dtime 1516 REAL paprs(klon, klevp1) 1517 REAL pplay(klon, klev) 1518 1519 ! Variables d'etat 1520 REAL t(klon, klev) 1521 REAL q(klon, klev) 1522 1523 ! Tendances 1524 REAL d_t(klon, klev) 1525 REAL d_q(klon, klev) 1526 1527 ! Profiles cible 1528 REAL t_targ(klon, klev) 1529 REAL rh_targ(klon, klev) 1530 1531 ! Temps de relaxation 1532 REAL tau 1533 !c DATA tau /3600./ 1534 !! DATA tau /5400./ 1535 DATA tau /1800./ 1536 1537 INTEGER k, i 1538 REAL zx_qs, rh, tnew, d_rh, rhnew 1539 1540 ! Declaration des constantes et des fonctions thermodynamiques 1541 1542 include "YOMCST.h" 1543 include "YOETHF.h" 1544 1545 ! ---------------------------------------- 1546 ! Statement functions 1547 include "FCTTRE.h" 1548 ! ---------------------------------------- 1549 1550 print *, 'dtime, tau ', dtime, tau 1551 print *, 't_targ', t_targ 1552 print *, 'rh_targ', rh_targ 1553 print *, 'temp ', t 1554 print *, 'hum ', q 1555 1556 DO k = 1, klev 1557 DO i = 1, klon 1558 IF (paprs(i, 1) - pplay(i, k) > 10000.) THEN 1559 IF (t(i, k)<RTT) THEN 1560 zx_qs = qsats(t(i, k)) / (pplay(i, k)) 1561 ELSE 1562 zx_qs = qsatl(t(i, k)) / (pplay(i, k)) 1563 ENDIF 1564 rh = q(i, k) / zx_qs 1565 1566 d_t(i, k) = d_t(i, k) + 1. / tau * (t_targ(i, k) - t(i, k)) 1567 d_rh = 1. / tau * (rh_targ(i, k) - rh) 1568 1569 tnew = t(i, k) + d_t(i, k) * dtime 1570 !jyg< 1571 ! Formule pour q : 1572 ! d_q = (1/tau) [rh_targ*qsat(T_new) - q] 1573 1574 ! Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new) 1575 ! qui n'etait pas correcte. 1576 1577 IF (tnew<RTT) THEN 1578 zx_qs = qsats(tnew) / (pplay(i, k)) 1579 ELSE 1580 zx_qs = qsatl(tnew) / (pplay(i, k)) 1581 ENDIF 1582 !! d_q(i,k) = d_q(i,k) + d_rh*zx_qs 1583 d_q(i, k) = d_q(i, k) + (1. / tau) * (rh_targ(i, k) * zx_qs - q(i, k)) 1584 rhnew = (q(i, k) + d_q(i, k) * dtime) / zx_qs 1585 1586 print *, ' k,d_t,rh,d_rh,rhnew,d_q ', & 1587 k, d_t(i, k), rh, d_rh, rhnew, d_q(i, k) 1634 1588 ENDIF 1635 !! d_q(i,k) = d_q(i,k) + d_rh*zx_qs 1636 d_q(i, k) = d_q(i, k) + (1. / tau) * (rh_targ(i, k) * zx_qs - q(i, k)) 1637 rhnew = (q(i, k) + d_q(i, k) * dtime) / zx_qs 1638 1639 print *, ' k,d_t,rh,d_rh,rhnew,d_q ', & 1640 k, d_t(i, k), rh, d_rh, rhnew, d_q(i, k) 1641 ENDIF 1642 1589 1590 ENDDO 1643 1591 ENDDO 1644 ENDDO 1645 1646 RETURN1647 END 1648 1649 Subroutine Nudge_UV (dtime, paprs, pplay, u_targ, v_targ, u, v, & 1650 & d_u, d_v)1651 ! ========================================================1652 USE dimphy 1653 1654 implicit none 1655 1656 ! ========================================================1657 REAL dtime1658 REAL paprs(klon, klevp1)1659 REAL pplay(klon, klev) 1660 1661 ! Variables d'etat1662 REAL u(klon, klev)1663 REAL v(klon, klev) 1664 1665 ! Tendances1666 REAL d_u(klon, klev)1667 REAL d_v(klon, klev) 1668 1669 ! Profiles cible1670 REAL u_targ(klon, klev)1671 REAL v_targ(klon, klev) 1672 1673 ! Temps de relaxation1674 REAL tau1675 !c DATA tau /3600./1676 ! DATA tau /5400./1677 DATA tau /43200./ 1678 1679 INTEGER k, i 1680 1681 !print *,'dtime, tau ',dtime,tau1682 !print *, 'u_targ',u_targ1683 !print *, 'v_targ',v_targ1684 !print *,'zonal velocity ',u1685 !print *,'meridional velocity ',v1686 DO k = 1, klev1687 DO i = 1, klon1688 !CR: nudging everywhere1689 ! IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN 1690 1691 d_u(i, k) = d_u(i, k) + 1. / tau * (u_targ(i, k) - u(i, k))1692 d_v(i, k) = d_v(i, k) + 1. / tau * (v_targ(i, k) - v(i, k)) 1693 1694 ! print *,' k,u,d_u,v,d_v ', &1695 ! k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)1696 ! ENDIF 1697 1592 1593 RETURN 1594 END 1595 1596 Subroutine Nudge_UV (dtime, paprs, pplay, u_targ, v_targ, u, v, & 1597 & d_u, d_v) 1598 ! ======================================================== 1599 USE dimphy 1600 1601 implicit none 1602 1603 ! ======================================================== 1604 REAL dtime 1605 REAL paprs(klon, klevp1) 1606 REAL pplay(klon, klev) 1607 1608 ! Variables d'etat 1609 REAL u(klon, klev) 1610 REAL v(klon, klev) 1611 1612 ! Tendances 1613 REAL d_u(klon, klev) 1614 REAL d_v(klon, klev) 1615 1616 ! Profiles cible 1617 REAL u_targ(klon, klev) 1618 REAL v_targ(klon, klev) 1619 1620 ! Temps de relaxation 1621 REAL tau 1622 !c DATA tau /3600./ 1623 ! DATA tau /5400./ 1624 DATA tau /43200./ 1625 1626 INTEGER k, i 1627 1628 !print *,'dtime, tau ',dtime,tau 1629 !print *, 'u_targ',u_targ 1630 !print *, 'v_targ',v_targ 1631 !print *,'zonal velocity ',u 1632 !print *,'meridional velocity ',v 1633 DO k = 1, klev 1634 DO i = 1, klon 1635 !CR: nudging everywhere 1636 ! IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN 1637 1638 d_u(i, k) = d_u(i, k) + 1. / tau * (u_targ(i, k) - u(i, k)) 1639 d_v(i, k) = d_v(i, k) + 1. / tau * (v_targ(i, k) - v(i, k)) 1640 1641 ! print *,' k,u,d_u,v,d_v ', & 1642 ! k,u(i,k),d_u(i,k),v(i,k),d_v(i,k) 1643 ! ENDIF 1644 1645 ENDDO 1698 1646 ENDDO 1699 ENDDO 1700 1701 RETURN 1702 END 1703 1704 !===================================================================== 1705 SUBROUTINE interp2_case_vertical(play, nlev_cas, plev_prof_cas & 1706 &, t_prof_cas, th_prof_cas, thv_prof_cas, thl_prof_cas & 1707 &, qv_prof_cas, ql_prof_cas, qi_prof_cas, u_prof_cas, v_prof_cas & 1708 &, ug_prof_cas, vg_prof_cas, vitw_prof_cas, omega_prof_cas & 1709 &, du_prof_cas, hu_prof_cas, vu_prof_cas, dv_prof_cas, hv_prof_cas, vv_prof_cas & 1710 &, dt_prof_cas, ht_prof_cas, vt_prof_cas, dtrad_prof_cas, dq_prof_cas, hq_prof_cas, vq_prof_cas & 1711 &, dth_prof_cas, hth_prof_cas, vth_prof_cas & 1712 1713 &, t_mod_cas, theta_mod_cas, thv_mod_cas, thl_mod_cas & 1714 &, qv_mod_cas, ql_mod_cas, qi_mod_cas, u_mod_cas, v_mod_cas & 1715 &, ug_mod_cas, vg_mod_cas, w_mod_cas, omega_mod_cas & 1716 &, du_mod_cas, hu_mod_cas, vu_mod_cas, dv_mod_cas, hv_mod_cas, vv_mod_cas & 1717 &, dt_mod_cas, ht_mod_cas, vt_mod_cas, dtrad_mod_cas, dq_mod_cas, hq_mod_cas, vq_mod_cas & 1718 &, dth_mod_cas, hth_mod_cas, vth_mod_cas, mxcalc) 1719 1720 implicit none 1721 1722 include "YOMCST.h" 1723 include "dimensions.h" 1724 1725 !------------------------------------------------------------------------- 1726 ! Vertical interpolation of generic case forcing data onto mod_casel levels 1727 !------------------------------------------------------------------------- 1728 1729 integer nlevmax 1730 parameter (nlevmax = 41) 1731 integer nlev_cas, mxcalc 1732 ! real play(llm), plev_prof(nlevmax) 1733 ! real t_prof(nlevmax),q_prof(nlevmax) 1734 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax) 1735 ! real ht_prof(nlevmax),vt_prof(nlevmax) 1736 ! real hq_prof(nlevmax),vq_prof(nlevmax) 1737 1738 real play(llm), plev_prof_cas(nlev_cas) 1739 real t_prof_cas(nlev_cas), th_prof_cas(nlev_cas), thv_prof_cas(nlev_cas), thl_prof_cas(nlev_cas) 1740 real qv_prof_cas(nlev_cas), ql_prof_cas(nlev_cas), qi_prof_cas(nlev_cas) 1741 real u_prof_cas(nlev_cas), v_prof_cas(nlev_cas) 1742 real ug_prof_cas(nlev_cas), vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas), omega_prof_cas(nlev_cas) 1743 real du_prof_cas(nlev_cas), hu_prof_cas(nlev_cas), vu_prof_cas(nlev_cas) 1744 real dv_prof_cas(nlev_cas), hv_prof_cas(nlev_cas), vv_prof_cas(nlev_cas) 1745 real dt_prof_cas(nlev_cas), ht_prof_cas(nlev_cas), vt_prof_cas(nlev_cas), dtrad_prof_cas(nlev_cas) 1746 real dth_prof_cas(nlev_cas), hth_prof_cas(nlev_cas), vth_prof_cas(nlev_cas) 1747 real dq_prof_cas(nlev_cas), hq_prof_cas(nlev_cas), vq_prof_cas(nlev_cas) 1748 1749 real t_mod_cas(llm), theta_mod_cas(llm), thv_mod_cas(llm), thl_mod_cas(llm) 1750 real qv_mod_cas(llm), ql_mod_cas(llm), qi_mod_cas(llm) 1751 real u_mod_cas(llm), v_mod_cas(llm) 1752 real ug_mod_cas(llm), vg_mod_cas(llm), w_mod_cas(llm), omega_mod_cas(llm) 1753 real du_mod_cas(llm), hu_mod_cas(llm), vu_mod_cas(llm) 1754 real dv_mod_cas(llm), hv_mod_cas(llm), vv_mod_cas(llm) 1755 real dt_mod_cas(llm), ht_mod_cas(llm), vt_mod_cas(llm), dtrad_mod_cas(llm) 1756 real dth_mod_cas(llm), hth_mod_cas(llm), vth_mod_cas(llm) 1757 real dq_mod_cas(llm), hq_mod_cas(llm), vq_mod_cas(llm) 1758 1759 integer l, k, k1, k2 1760 real frac, frac1, frac2, fact 1761 1762 ! do l = 1, llm 1763 ! print *,'debut interp2, play=',l,play(l) 1764 ! enddo 1765 ! do l = 1, nlev_cas 1766 ! print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l) 1767 ! enddo 1768 1769 do l = 1, llm 1770 1771 if (play(l)>=plev_prof_cas(nlev_cas)) then 1772 1773 mxcalc = l 1774 ! print *,'debut interp2, mxcalc=',mxcalc 1775 k1 = 0 1776 k2 = 0 1777 1778 if (play(l)<=plev_prof_cas(1)) then 1779 1780 do k = 1, nlev_cas - 1 1781 if (play(l)<=plev_prof_cas(k).and. play(l)>plev_prof_cas(k + 1)) then 1782 k1 = k 1783 k2 = k + 1 1647 1648 RETURN 1649 END 1650 1651 SUBROUTINE interp2_case_vertical(play, nlev_cas, plev_prof_cas & 1652 &, t_prof_cas, th_prof_cas, thv_prof_cas, thl_prof_cas & 1653 &, qv_prof_cas, ql_prof_cas, qi_prof_cas, u_prof_cas, v_prof_cas & 1654 &, ug_prof_cas, vg_prof_cas, vitw_prof_cas, omega_prof_cas & 1655 &, du_prof_cas, hu_prof_cas, vu_prof_cas, dv_prof_cas, hv_prof_cas, vv_prof_cas & 1656 &, dt_prof_cas, ht_prof_cas, vt_prof_cas, dtrad_prof_cas, dq_prof_cas, hq_prof_cas, vq_prof_cas & 1657 &, dth_prof_cas, hth_prof_cas, vth_prof_cas & 1658 1659 &, t_mod_cas, theta_mod_cas, thv_mod_cas, thl_mod_cas & 1660 &, qv_mod_cas, ql_mod_cas, qi_mod_cas, u_mod_cas, v_mod_cas & 1661 &, ug_mod_cas, vg_mod_cas, w_mod_cas, omega_mod_cas & 1662 &, du_mod_cas, hu_mod_cas, vu_mod_cas, dv_mod_cas, hv_mod_cas, vv_mod_cas & 1663 &, dt_mod_cas, ht_mod_cas, vt_mod_cas, dtrad_mod_cas, dq_mod_cas, hq_mod_cas, vq_mod_cas & 1664 &, dth_mod_cas, hth_mod_cas, vth_mod_cas, mxcalc) 1665 1666 implicit none 1667 1668 include "YOMCST.h" 1669 include "dimensions.h" 1670 1671 !------------------------------------------------------------------------- 1672 ! Vertical interpolation of generic case forcing data onto mod_casel levels 1673 !------------------------------------------------------------------------- 1674 1675 integer nlevmax 1676 parameter (nlevmax = 41) 1677 integer nlev_cas, mxcalc 1678 ! real play(llm), plev_prof(nlevmax) 1679 ! real t_prof(nlevmax),q_prof(nlevmax) 1680 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax) 1681 ! real ht_prof(nlevmax),vt_prof(nlevmax) 1682 ! real hq_prof(nlevmax),vq_prof(nlevmax) 1683 1684 real play(llm), plev_prof_cas(nlev_cas) 1685 real t_prof_cas(nlev_cas), th_prof_cas(nlev_cas), thv_prof_cas(nlev_cas), thl_prof_cas(nlev_cas) 1686 real qv_prof_cas(nlev_cas), ql_prof_cas(nlev_cas), qi_prof_cas(nlev_cas) 1687 real u_prof_cas(nlev_cas), v_prof_cas(nlev_cas) 1688 real ug_prof_cas(nlev_cas), vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas), omega_prof_cas(nlev_cas) 1689 real du_prof_cas(nlev_cas), hu_prof_cas(nlev_cas), vu_prof_cas(nlev_cas) 1690 real dv_prof_cas(nlev_cas), hv_prof_cas(nlev_cas), vv_prof_cas(nlev_cas) 1691 real dt_prof_cas(nlev_cas), ht_prof_cas(nlev_cas), vt_prof_cas(nlev_cas), dtrad_prof_cas(nlev_cas) 1692 real dth_prof_cas(nlev_cas), hth_prof_cas(nlev_cas), vth_prof_cas(nlev_cas) 1693 real dq_prof_cas(nlev_cas), hq_prof_cas(nlev_cas), vq_prof_cas(nlev_cas) 1694 1695 real t_mod_cas(llm), theta_mod_cas(llm), thv_mod_cas(llm), thl_mod_cas(llm) 1696 real qv_mod_cas(llm), ql_mod_cas(llm), qi_mod_cas(llm) 1697 real u_mod_cas(llm), v_mod_cas(llm) 1698 real ug_mod_cas(llm), vg_mod_cas(llm), w_mod_cas(llm), omega_mod_cas(llm) 1699 real du_mod_cas(llm), hu_mod_cas(llm), vu_mod_cas(llm) 1700 real dv_mod_cas(llm), hv_mod_cas(llm), vv_mod_cas(llm) 1701 real dt_mod_cas(llm), ht_mod_cas(llm), vt_mod_cas(llm), dtrad_mod_cas(llm) 1702 real dth_mod_cas(llm), hth_mod_cas(llm), vth_mod_cas(llm) 1703 real dq_mod_cas(llm), hq_mod_cas(llm), vq_mod_cas(llm) 1704 1705 integer l, k, k1, k2 1706 real frac, frac1, frac2, fact 1707 1708 ! do l = 1, llm 1709 ! print *,'debut interp2, play=',l,play(l) 1710 ! enddo 1711 ! do l = 1, nlev_cas 1712 ! print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l) 1713 ! enddo 1714 1715 do l = 1, llm 1716 1717 if (play(l)>=plev_prof_cas(nlev_cas)) then 1718 1719 mxcalc = l 1720 ! print *,'debut interp2, mxcalc=',mxcalc 1721 k1 = 0 1722 k2 = 0 1723 1724 if (play(l)<=plev_prof_cas(1)) then 1725 1726 do k = 1, nlev_cas - 1 1727 if (play(l)<=plev_prof_cas(k).and. play(l)>plev_prof_cas(k + 1)) then 1728 k1 = k 1729 k2 = k + 1 1730 endif 1731 enddo 1732 1733 if (k1==0 .or. k2==0) then 1734 write(*, *) 'PB! k1, k2 = ', k1, k2 1735 write(*, *) 'l,play(l) = ', l, play(l) / 100 1736 do k = 1, nlev_cas - 1 1737 write(*, *) 'k,plev_prof_cas(k) = ', k, plev_prof_cas(k) / 100 1738 enddo 1784 1739 endif 1785 enddo 1786 1787 if (k1==0 .or. k2==0) then 1788 write(*, *) 'PB! k1, k2 = ', k1, k2 1789 write(*, *) 'l,play(l) = ', l, play(l) / 100 1790 do k = 1, nlev_cas - 1 1791 write(*, *) 'k,plev_prof_cas(k) = ', k, plev_prof_cas(k) / 100 1792 enddo 1793 endif 1794 1795 frac = (plev_prof_cas(k2) - play(l)) / (plev_prof_cas(k2) - plev_prof_cas(k1)) 1796 t_mod_cas(l) = t_prof_cas(k2) - frac * (t_prof_cas(k2) - t_prof_cas(k1)) 1797 theta_mod_cas(l) = th_prof_cas(k2) - frac * (th_prof_cas(k2) - th_prof_cas(k1)) 1798 if(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD) 1799 thv_mod_cas(l) = thv_prof_cas(k2) - frac * (thv_prof_cas(k2) - thv_prof_cas(k1)) 1800 thl_mod_cas(l) = thl_prof_cas(k2) - frac * (thl_prof_cas(k2) - thl_prof_cas(k1)) 1801 qv_mod_cas(l) = qv_prof_cas(k2) - frac * (qv_prof_cas(k2) - qv_prof_cas(k1)) 1802 ql_mod_cas(l) = ql_prof_cas(k2) - frac * (ql_prof_cas(k2) - ql_prof_cas(k1)) 1803 qi_mod_cas(l) = qi_prof_cas(k2) - frac * (qi_prof_cas(k2) - qi_prof_cas(k1)) 1804 u_mod_cas(l) = u_prof_cas(k2) - frac * (u_prof_cas(k2) - u_prof_cas(k1)) 1805 v_mod_cas(l) = v_prof_cas(k2) - frac * (v_prof_cas(k2) - v_prof_cas(k1)) 1806 ug_mod_cas(l) = ug_prof_cas(k2) - frac * (ug_prof_cas(k2) - ug_prof_cas(k1)) 1807 vg_mod_cas(l) = vg_prof_cas(k2) - frac * (vg_prof_cas(k2) - vg_prof_cas(k1)) 1808 w_mod_cas(l) = vitw_prof_cas(k2) - frac * (vitw_prof_cas(k2) - vitw_prof_cas(k1)) 1809 omega_mod_cas(l) = omega_prof_cas(k2) - frac * (omega_prof_cas(k2) - omega_prof_cas(k1)) 1810 du_mod_cas(l) = du_prof_cas(k2) - frac * (du_prof_cas(k2) - du_prof_cas(k1)) 1811 hu_mod_cas(l) = hu_prof_cas(k2) - frac * (hu_prof_cas(k2) - hu_prof_cas(k1)) 1812 vu_mod_cas(l) = vu_prof_cas(k2) - frac * (vu_prof_cas(k2) - vu_prof_cas(k1)) 1813 dv_mod_cas(l) = dv_prof_cas(k2) - frac * (dv_prof_cas(k2) - dv_prof_cas(k1)) 1814 hv_mod_cas(l) = hv_prof_cas(k2) - frac * (hv_prof_cas(k2) - hv_prof_cas(k1)) 1815 vv_mod_cas(l) = vv_prof_cas(k2) - frac * (vv_prof_cas(k2) - vv_prof_cas(k1)) 1816 dt_mod_cas(l) = dt_prof_cas(k2) - frac * (dt_prof_cas(k2) - dt_prof_cas(k1)) 1817 ht_mod_cas(l) = ht_prof_cas(k2) - frac * (ht_prof_cas(k2) - ht_prof_cas(k1)) 1818 vt_mod_cas(l) = vt_prof_cas(k2) - frac * (vt_prof_cas(k2) - vt_prof_cas(k1)) 1819 dth_mod_cas(l) = dth_prof_cas(k2) - frac * (dth_prof_cas(k2) - dth_prof_cas(k1)) 1820 hth_mod_cas(l) = hth_prof_cas(k2) - frac * (hth_prof_cas(k2) - hth_prof_cas(k1)) 1821 vth_mod_cas(l) = vth_prof_cas(k2) - frac * (vth_prof_cas(k2) - vth_prof_cas(k1)) 1822 dq_mod_cas(l) = dq_prof_cas(k2) - frac * (dq_prof_cas(k2) - dq_prof_cas(k1)) 1823 hq_mod_cas(l) = hq_prof_cas(k2) - frac * (hq_prof_cas(k2) - hq_prof_cas(k1)) 1824 vq_mod_cas(l) = vq_prof_cas(k2) - frac * (vq_prof_cas(k2) - vq_prof_cas(k1)) 1825 dtrad_mod_cas(l) = dtrad_prof_cas(k2) - frac * (dtrad_prof_cas(k2) - dtrad_prof_cas(k1)) 1826 1827 else !play>plev_prof_cas(1) 1828 1829 k1 = 1 1830 k2 = 2 1831 print *, 'interp2_vert, k1,k2=', plev_prof_cas(k1), plev_prof_cas(k2) 1832 frac1 = (play(l) - plev_prof_cas(k2)) / (plev_prof_cas(k1) - plev_prof_cas(k2)) 1833 frac2 = (play(l) - plev_prof_cas(k1)) / (plev_prof_cas(k1) - plev_prof_cas(k2)) 1834 t_mod_cas(l) = frac1 * t_prof_cas(k1) - frac2 * t_prof_cas(k2) 1835 theta_mod_cas(l) = frac1 * th_prof_cas(k1) - frac2 * th_prof_cas(k2) 1836 if(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD) 1837 thv_mod_cas(l) = frac1 * thv_prof_cas(k1) - frac2 * thv_prof_cas(k2) 1838 thl_mod_cas(l) = frac1 * thl_prof_cas(k1) - frac2 * thl_prof_cas(k2) 1839 qv_mod_cas(l) = frac1 * qv_prof_cas(k1) - frac2 * qv_prof_cas(k2) 1840 ql_mod_cas(l) = frac1 * ql_prof_cas(k1) - frac2 * ql_prof_cas(k2) 1841 qi_mod_cas(l) = frac1 * qi_prof_cas(k1) - frac2 * qi_prof_cas(k2) 1842 u_mod_cas(l) = frac1 * u_prof_cas(k1) - frac2 * u_prof_cas(k2) 1843 v_mod_cas(l) = frac1 * v_prof_cas(k1) - frac2 * v_prof_cas(k2) 1844 ug_mod_cas(l) = frac1 * ug_prof_cas(k1) - frac2 * ug_prof_cas(k2) 1845 vg_mod_cas(l) = frac1 * vg_prof_cas(k1) - frac2 * vg_prof_cas(k2) 1846 w_mod_cas(l) = frac1 * vitw_prof_cas(k1) - frac2 * vitw_prof_cas(k2) 1847 omega_mod_cas(l) = frac1 * omega_prof_cas(k1) - frac2 * omega_prof_cas(k2) 1848 du_mod_cas(l) = frac1 * du_prof_cas(k1) - frac2 * du_prof_cas(k2) 1849 hu_mod_cas(l) = frac1 * hu_prof_cas(k1) - frac2 * hu_prof_cas(k2) 1850 vu_mod_cas(l) = frac1 * vu_prof_cas(k1) - frac2 * vu_prof_cas(k2) 1851 dv_mod_cas(l) = frac1 * dv_prof_cas(k1) - frac2 * dv_prof_cas(k2) 1852 hv_mod_cas(l) = frac1 * hv_prof_cas(k1) - frac2 * hv_prof_cas(k2) 1853 vv_mod_cas(l) = frac1 * vv_prof_cas(k1) - frac2 * vv_prof_cas(k2) 1854 dt_mod_cas(l) = frac1 * dt_prof_cas(k1) - frac2 * dt_prof_cas(k2) 1855 ht_mod_cas(l) = frac1 * ht_prof_cas(k1) - frac2 * ht_prof_cas(k2) 1856 vt_mod_cas(l) = frac1 * vt_prof_cas(k1) - frac2 * vt_prof_cas(k2) 1857 dth_mod_cas(l) = frac1 * dth_prof_cas(k1) - frac2 * dth_prof_cas(k2) 1858 hth_mod_cas(l) = frac1 * hth_prof_cas(k1) - frac2 * hth_prof_cas(k2) 1859 vth_mod_cas(l) = frac1 * vth_prof_cas(k1) - frac2 * vth_prof_cas(k2) 1860 dq_mod_cas(l) = frac1 * dq_prof_cas(k1) - frac2 * dq_prof_cas(k2) 1861 hq_mod_cas(l) = frac1 * hq_prof_cas(k1) - frac2 * hq_prof_cas(k2) 1862 vq_mod_cas(l) = frac1 * vq_prof_cas(k1) - frac2 * vq_prof_cas(k2) 1863 dtrad_mod_cas(l) = frac1 * dtrad_prof_cas(k1) - frac2 * dtrad_prof_cas(k2) 1864 1865 endif ! play.le.plev_prof_cas(1) 1866 1867 else ! above max altitude of forcing file 1868 1869 !jyg 1870 fact = 20. * (plev_prof_cas(nlev_cas) - play(l)) / plev_prof_cas(nlev_cas) !jyg 1871 fact = max(fact, 0.) !jyg 1872 fact = exp(-fact) !jyg 1873 t_mod_cas(l) = t_prof_cas(nlev_cas) !jyg 1874 theta_mod_cas(l) = th_prof_cas(nlev_cas) !jyg 1875 thv_mod_cas(l) = thv_prof_cas(nlev_cas) !jyg 1876 thl_mod_cas(l) = thl_prof_cas(nlev_cas) !jyg 1877 qv_mod_cas(l) = qv_prof_cas(nlev_cas) * fact !jyg 1878 ql_mod_cas(l) = ql_prof_cas(nlev_cas) * fact !jyg 1879 qi_mod_cas(l) = qi_prof_cas(nlev_cas) * fact !jyg 1880 u_mod_cas(l) = u_prof_cas(nlev_cas) * fact !jyg 1881 v_mod_cas(l) = v_prof_cas(nlev_cas) * fact !jyg 1882 ug_mod_cas(l) = ug_prof_cas(nlev_cas) * fact !jyg 1883 vg_mod_cas(l) = vg_prof_cas(nlev_cas) * fact !jyg 1884 w_mod_cas(l) = 0.0 !jyg 1885 omega_mod_cas(l) = 0.0 !jyg 1886 du_mod_cas(l) = du_prof_cas(nlev_cas) * fact 1887 hu_mod_cas(l) = hu_prof_cas(nlev_cas) * fact !jyg 1888 vu_mod_cas(l) = vu_prof_cas(nlev_cas) * fact !jyg 1889 dv_mod_cas(l) = dv_prof_cas(nlev_cas) * fact 1890 hv_mod_cas(l) = hv_prof_cas(nlev_cas) * fact !jyg 1891 vv_mod_cas(l) = vv_prof_cas(nlev_cas) * fact !jyg 1892 dt_mod_cas(l) = dt_prof_cas(nlev_cas) 1893 ht_mod_cas(l) = ht_prof_cas(nlev_cas) !jyg 1894 vt_mod_cas(l) = vt_prof_cas(nlev_cas) !jyg 1895 dth_mod_cas(l) = dth_prof_cas(nlev_cas) 1896 hth_mod_cas(l) = hth_prof_cas(nlev_cas) !jyg 1897 vth_mod_cas(l) = vth_prof_cas(nlev_cas) !jyg 1898 dq_mod_cas(l) = dq_prof_cas(nlev_cas) * fact 1899 hq_mod_cas(l) = hq_prof_cas(nlev_cas) * fact !jyg 1900 vq_mod_cas(l) = vq_prof_cas(nlev_cas) * fact !jyg 1901 dtrad_mod_cas(l) = dtrad_prof_cas(nlev_cas) * fact !jyg 1902 1903 endif ! play 1904 1905 enddo ! l 1906 1907 return 1908 end 1909 !***************************************************************************** 1910 1911 1912 1913 1740 1741 frac = (plev_prof_cas(k2) - play(l)) / (plev_prof_cas(k2) - plev_prof_cas(k1)) 1742 t_mod_cas(l) = t_prof_cas(k2) - frac * (t_prof_cas(k2) - t_prof_cas(k1)) 1743 theta_mod_cas(l) = th_prof_cas(k2) - frac * (th_prof_cas(k2) - th_prof_cas(k1)) 1744 if(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD) 1745 thv_mod_cas(l) = thv_prof_cas(k2) - frac * (thv_prof_cas(k2) - thv_prof_cas(k1)) 1746 thl_mod_cas(l) = thl_prof_cas(k2) - frac * (thl_prof_cas(k2) - thl_prof_cas(k1)) 1747 qv_mod_cas(l) = qv_prof_cas(k2) - frac * (qv_prof_cas(k2) - qv_prof_cas(k1)) 1748 ql_mod_cas(l) = ql_prof_cas(k2) - frac * (ql_prof_cas(k2) - ql_prof_cas(k1)) 1749 qi_mod_cas(l) = qi_prof_cas(k2) - frac * (qi_prof_cas(k2) - qi_prof_cas(k1)) 1750 u_mod_cas(l) = u_prof_cas(k2) - frac * (u_prof_cas(k2) - u_prof_cas(k1)) 1751 v_mod_cas(l) = v_prof_cas(k2) - frac * (v_prof_cas(k2) - v_prof_cas(k1)) 1752 ug_mod_cas(l) = ug_prof_cas(k2) - frac * (ug_prof_cas(k2) - ug_prof_cas(k1)) 1753 vg_mod_cas(l) = vg_prof_cas(k2) - frac * (vg_prof_cas(k2) - vg_prof_cas(k1)) 1754 w_mod_cas(l) = vitw_prof_cas(k2) - frac * (vitw_prof_cas(k2) - vitw_prof_cas(k1)) 1755 omega_mod_cas(l) = omega_prof_cas(k2) - frac * (omega_prof_cas(k2) - omega_prof_cas(k1)) 1756 du_mod_cas(l) = du_prof_cas(k2) - frac * (du_prof_cas(k2) - du_prof_cas(k1)) 1757 hu_mod_cas(l) = hu_prof_cas(k2) - frac * (hu_prof_cas(k2) - hu_prof_cas(k1)) 1758 vu_mod_cas(l) = vu_prof_cas(k2) - frac * (vu_prof_cas(k2) - vu_prof_cas(k1)) 1759 dv_mod_cas(l) = dv_prof_cas(k2) - frac * (dv_prof_cas(k2) - dv_prof_cas(k1)) 1760 hv_mod_cas(l) = hv_prof_cas(k2) - frac * (hv_prof_cas(k2) - hv_prof_cas(k1)) 1761 vv_mod_cas(l) = vv_prof_cas(k2) - frac * (vv_prof_cas(k2) - vv_prof_cas(k1)) 1762 dt_mod_cas(l) = dt_prof_cas(k2) - frac * (dt_prof_cas(k2) - dt_prof_cas(k1)) 1763 ht_mod_cas(l) = ht_prof_cas(k2) - frac * (ht_prof_cas(k2) - ht_prof_cas(k1)) 1764 vt_mod_cas(l) = vt_prof_cas(k2) - frac * (vt_prof_cas(k2) - vt_prof_cas(k1)) 1765 dth_mod_cas(l) = dth_prof_cas(k2) - frac * (dth_prof_cas(k2) - dth_prof_cas(k1)) 1766 hth_mod_cas(l) = hth_prof_cas(k2) - frac * (hth_prof_cas(k2) - hth_prof_cas(k1)) 1767 vth_mod_cas(l) = vth_prof_cas(k2) - frac * (vth_prof_cas(k2) - vth_prof_cas(k1)) 1768 dq_mod_cas(l) = dq_prof_cas(k2) - frac * (dq_prof_cas(k2) - dq_prof_cas(k1)) 1769 hq_mod_cas(l) = hq_prof_cas(k2) - frac * (hq_prof_cas(k2) - hq_prof_cas(k1)) 1770 vq_mod_cas(l) = vq_prof_cas(k2) - frac * (vq_prof_cas(k2) - vq_prof_cas(k1)) 1771 dtrad_mod_cas(l) = dtrad_prof_cas(k2) - frac * (dtrad_prof_cas(k2) - dtrad_prof_cas(k1)) 1772 1773 else !play>plev_prof_cas(1) 1774 1775 k1 = 1 1776 k2 = 2 1777 print *, 'interp2_vert, k1,k2=', plev_prof_cas(k1), plev_prof_cas(k2) 1778 frac1 = (play(l) - plev_prof_cas(k2)) / (plev_prof_cas(k1) - plev_prof_cas(k2)) 1779 frac2 = (play(l) - plev_prof_cas(k1)) / (plev_prof_cas(k1) - plev_prof_cas(k2)) 1780 t_mod_cas(l) = frac1 * t_prof_cas(k1) - frac2 * t_prof_cas(k2) 1781 theta_mod_cas(l) = frac1 * th_prof_cas(k1) - frac2 * th_prof_cas(k2) 1782 if(theta_mod_cas(l)/=0) t_mod_cas(l) = theta_mod_cas(l) * (play(l) / 100000.)**(RD / RCPD) 1783 thv_mod_cas(l) = frac1 * thv_prof_cas(k1) - frac2 * thv_prof_cas(k2) 1784 thl_mod_cas(l) = frac1 * thl_prof_cas(k1) - frac2 * thl_prof_cas(k2) 1785 qv_mod_cas(l) = frac1 * qv_prof_cas(k1) - frac2 * qv_prof_cas(k2) 1786 ql_mod_cas(l) = frac1 * ql_prof_cas(k1) - frac2 * ql_prof_cas(k2) 1787 qi_mod_cas(l) = frac1 * qi_prof_cas(k1) - frac2 * qi_prof_cas(k2) 1788 u_mod_cas(l) = frac1 * u_prof_cas(k1) - frac2 * u_prof_cas(k2) 1789 v_mod_cas(l) = frac1 * v_prof_cas(k1) - frac2 * v_prof_cas(k2) 1790 ug_mod_cas(l) = frac1 * ug_prof_cas(k1) - frac2 * ug_prof_cas(k2) 1791 vg_mod_cas(l) = frac1 * vg_prof_cas(k1) - frac2 * vg_prof_cas(k2) 1792 w_mod_cas(l) = frac1 * vitw_prof_cas(k1) - frac2 * vitw_prof_cas(k2) 1793 omega_mod_cas(l) = frac1 * omega_prof_cas(k1) - frac2 * omega_prof_cas(k2) 1794 du_mod_cas(l) = frac1 * du_prof_cas(k1) - frac2 * du_prof_cas(k2) 1795 hu_mod_cas(l) = frac1 * hu_prof_cas(k1) - frac2 * hu_prof_cas(k2) 1796 vu_mod_cas(l) = frac1 * vu_prof_cas(k1) - frac2 * vu_prof_cas(k2) 1797 dv_mod_cas(l) = frac1 * dv_prof_cas(k1) - frac2 * dv_prof_cas(k2) 1798 hv_mod_cas(l) = frac1 * hv_prof_cas(k1) - frac2 * hv_prof_cas(k2) 1799 vv_mod_cas(l) = frac1 * vv_prof_cas(k1) - frac2 * vv_prof_cas(k2) 1800 dt_mod_cas(l) = frac1 * dt_prof_cas(k1) - frac2 * dt_prof_cas(k2) 1801 ht_mod_cas(l) = frac1 * ht_prof_cas(k1) - frac2 * ht_prof_cas(k2) 1802 vt_mod_cas(l) = frac1 * vt_prof_cas(k1) - frac2 * vt_prof_cas(k2) 1803 dth_mod_cas(l) = frac1 * dth_prof_cas(k1) - frac2 * dth_prof_cas(k2) 1804 hth_mod_cas(l) = frac1 * hth_prof_cas(k1) - frac2 * hth_prof_cas(k2) 1805 vth_mod_cas(l) = frac1 * vth_prof_cas(k1) - frac2 * vth_prof_cas(k2) 1806 dq_mod_cas(l) = frac1 * dq_prof_cas(k1) - frac2 * dq_prof_cas(k2) 1807 hq_mod_cas(l) = frac1 * hq_prof_cas(k1) - frac2 * hq_prof_cas(k2) 1808 vq_mod_cas(l) = frac1 * vq_prof_cas(k1) - frac2 * vq_prof_cas(k2) 1809 dtrad_mod_cas(l) = frac1 * dtrad_prof_cas(k1) - frac2 * dtrad_prof_cas(k2) 1810 1811 endif ! play.le.plev_prof_cas(1) 1812 1813 else ! above max altitude of forcing file 1814 1815 !jyg 1816 fact = 20. * (plev_prof_cas(nlev_cas) - play(l)) / plev_prof_cas(nlev_cas) !jyg 1817 fact = max(fact, 0.) !jyg 1818 fact = exp(-fact) !jyg 1819 t_mod_cas(l) = t_prof_cas(nlev_cas) !jyg 1820 theta_mod_cas(l) = th_prof_cas(nlev_cas) !jyg 1821 thv_mod_cas(l) = thv_prof_cas(nlev_cas) !jyg 1822 thl_mod_cas(l) = thl_prof_cas(nlev_cas) !jyg 1823 qv_mod_cas(l) = qv_prof_cas(nlev_cas) * fact !jyg 1824 ql_mod_cas(l) = ql_prof_cas(nlev_cas) * fact !jyg 1825 qi_mod_cas(l) = qi_prof_cas(nlev_cas) * fact !jyg 1826 u_mod_cas(l) = u_prof_cas(nlev_cas) * fact !jyg 1827 v_mod_cas(l) = v_prof_cas(nlev_cas) * fact !jyg 1828 ug_mod_cas(l) = ug_prof_cas(nlev_cas) * fact !jyg 1829 vg_mod_cas(l) = vg_prof_cas(nlev_cas) * fact !jyg 1830 w_mod_cas(l) = 0.0 !jyg 1831 omega_mod_cas(l) = 0.0 !jyg 1832 du_mod_cas(l) = du_prof_cas(nlev_cas) * fact 1833 hu_mod_cas(l) = hu_prof_cas(nlev_cas) * fact !jyg 1834 vu_mod_cas(l) = vu_prof_cas(nlev_cas) * fact !jyg 1835 dv_mod_cas(l) = dv_prof_cas(nlev_cas) * fact 1836 hv_mod_cas(l) = hv_prof_cas(nlev_cas) * fact !jyg 1837 vv_mod_cas(l) = vv_prof_cas(nlev_cas) * fact !jyg 1838 dt_mod_cas(l) = dt_prof_cas(nlev_cas) 1839 ht_mod_cas(l) = ht_prof_cas(nlev_cas) !jyg 1840 vt_mod_cas(l) = vt_prof_cas(nlev_cas) !jyg 1841 dth_mod_cas(l) = dth_prof_cas(nlev_cas) 1842 hth_mod_cas(l) = hth_prof_cas(nlev_cas) !jyg 1843 vth_mod_cas(l) = vth_prof_cas(nlev_cas) !jyg 1844 dq_mod_cas(l) = dq_prof_cas(nlev_cas) * fact 1845 hq_mod_cas(l) = hq_prof_cas(nlev_cas) * fact !jyg 1846 vq_mod_cas(l) = vq_prof_cas(nlev_cas) * fact !jyg 1847 dtrad_mod_cas(l) = dtrad_prof_cas(nlev_cas) * fact !jyg 1848 1849 endif ! play 1850 1851 enddo ! l 1852 1853 return 1854 end 1855 1856 END MODULE lmdz_1dutils
Note: See TracChangeset
for help on using the changeset viewer.
