[887] | 1 | PROGRAM rcm1d |
---|
| 2 | |
---|
| 3 | USE infotrac |
---|
[3540] | 4 | use control_mod, only: planet_type,day_step,cpofT |
---|
[901] | 5 | USE phys_state_var_mod |
---|
[1305] | 6 | use chemparam_mod |
---|
[1517] | 7 | USE comconst_mod, ONLY: cpp,t0_venus,nu_venus |
---|
[1017] | 8 | use cpdet_mod, only: ini_cpdet |
---|
[1301] | 9 | use moyzon_mod, only: tmoy |
---|
[1442] | 10 | USE comvert_mod, ONLY: ap,bp,presnivs,pa,preff,nivsigs,nivsig, |
---|
| 11 | . aps,bps,scaleheight,pseudoalt, |
---|
| 12 | . disvert_type,pressure_exner |
---|
[2684] | 13 | use conc, only: rho,mmean |
---|
[1523] | 14 | USE iniphysiq_mod, ONLY: iniphysiq |
---|
[1543] | 15 | USE mod_const_mpi, ONLY: comm_lmdz |
---|
[1549] | 16 | USE physiq_mod, ONLY: physiq |
---|
[1692] | 17 | USE logic_mod, ONLY: iflag_trac |
---|
[3255] | 18 | ! For XIOS outputs: |
---|
| 19 | USE mod_const_mpi, ONLY: init_const_mpi |
---|
| 20 | USE parallel_lmdz, ONLY: init_parallel |
---|
[887] | 21 | IMPLICIT NONE |
---|
| 22 | |
---|
| 23 | c======================================================================= |
---|
| 24 | c subject: |
---|
| 25 | c -------- |
---|
| 26 | c PROGRAM useful to run physical part of the venusian GCM in a 1D column |
---|
| 27 | c |
---|
| 28 | c Can be compiled with a command like (e.g. for 50 layers) |
---|
[3255] | 29 | c "makelmdz_lmdz -p venus -d 50 rcm1d" |
---|
| 30 | c If you want XIOS outputs, then you'll need to compile with MPI and xios: |
---|
| 31 | c "makelmdz_lmdz -p venus -parallel mpi -io xios -d 50 rcm1d" |
---|
| 32 | c but the model should then be run using a single core, i.e. without mpirun |
---|
[887] | 33 | |
---|
| 34 | c It requires the files "rcm1d.def" "physiq.def" |
---|
| 35 | c and a file describing the sigma layers (e.g. "z2sig.def") |
---|
| 36 | c |
---|
| 37 | c author: Frederic Hourdin, R.Fournier,F.Forget (original Mars version) |
---|
| 38 | c ------- Sebastien Lebonnois (Venus version) |
---|
| 39 | c |
---|
| 40 | c======================================================================= |
---|
| 41 | |
---|
| 42 | #include "dimensions.h" |
---|
| 43 | #include "dimsoil.h" |
---|
| 44 | #include "comcstfi.h" |
---|
| 45 | #include "netcdf.inc" |
---|
| 46 | #include "clesphys.h" |
---|
| 47 | #include "iniprint.h" |
---|
| 48 | #include "tabcontrol.h" |
---|
| 49 | |
---|
| 50 | c -------------------------------------------------------------- |
---|
| 51 | c Declarations |
---|
| 52 | c -------------------------------------------------------------- |
---|
| 53 | c |
---|
| 54 | INTEGER unit ! unite de lecture de "rcm1d.def" |
---|
| 55 | INTEGER unitstart ! unite d'ecriture de "startphy.nc" |
---|
| 56 | INTEGER nlayer,nlevel,nsoil,ndt |
---|
[2684] | 57 | INTEGER ilayer,ilevel,isoil,idt,iq,i |
---|
[887] | 58 | LOGICAl firstcall,lastcall |
---|
[3254] | 59 | REAL :: nb_days ! number of Vdays (and/or fraction thererof) to run |
---|
[887] | 60 | c |
---|
| 61 | INTEGER day0 ! date initial (sol ; =0 a Ls=0) |
---|
| 62 | REAL day ! date durant le run |
---|
| 63 | REAL time ! time (0<time<1 ; time=0.5 a midi) |
---|
| 64 | REAL play(llm) ! Pressure at the middle of the layers (Pa) |
---|
| 65 | REAL plev(llm+1) ! intermediate pressure levels (pa) |
---|
[901] | 66 | REAL psurf |
---|
[887] | 67 | REAL u(llm),v(llm) ! zonal, meridional wind |
---|
| 68 | REAL gru,grv ! prescribed "geostrophic" background wind |
---|
| 69 | REAL temp(llm) ! temperature at the middle of the layers |
---|
[894] | 70 | REAL,allocatable :: q(:,:) ! tracer mixing ratio (e.g. kg/kg) |
---|
[887] | 71 | REAL zlay(llm) ! altitude estimee dans les couches (km) |
---|
| 72 | REAL long(1),lati(1),area(1) |
---|
| 73 | REAL cufi(1),cvfi(1) |
---|
[901] | 74 | REAL phisfi(1) |
---|
[887] | 75 | |
---|
| 76 | c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s) |
---|
| 77 | REAL du(llm),dv(llm),dtemp(llm) |
---|
| 78 | REAL dudyn(llm),dvdyn(llm),dtempdyn(llm) |
---|
[1549] | 79 | REAL dpsurf(1) |
---|
[894] | 80 | REAL,allocatable :: dq(:,:) |
---|
[887] | 81 | |
---|
| 82 | c Various intermediate variables |
---|
| 83 | REAL zls |
---|
| 84 | REAL phi(llm),s(llm) |
---|
| 85 | REAL pk(llm),pks, w(llm) |
---|
| 86 | INTEGER l, ierr, aslun |
---|
| 87 | REAL tmp1(0:llm),tmp2(0:llm) |
---|
| 88 | |
---|
| 89 | character*2 str2 |
---|
| 90 | |
---|
[894] | 91 | real pi |
---|
[887] | 92 | |
---|
[2684] | 93 | ! initialisation des traceurs |
---|
[3254] | 94 | logical :: file_is_present |
---|
[2684] | 95 | integer :: idummy |
---|
| 96 | real :: dummy |
---|
| 97 | |
---|
[887] | 98 | c======================================================================= |
---|
| 99 | c INITIALISATION |
---|
| 100 | |
---|
| 101 | lunout = 6 |
---|
| 102 | |
---|
[3255] | 103 | #ifdef CPP_XIOS |
---|
| 104 | call init_const_mpi |
---|
| 105 | call init_parallel |
---|
| 106 | #endif |
---|
| 107 | |
---|
[887] | 108 | c ------------------------------------------------------ |
---|
| 109 | c Constantes prescrites ICI |
---|
| 110 | c ------------------------------------------------------ |
---|
| 111 | |
---|
| 112 | pi=2.E+0*asin(1.E+0) |
---|
| 113 | |
---|
| 114 | c Constante de la Planete Venus |
---|
| 115 | c ----------------------------- |
---|
| 116 | planet_type = "venus" |
---|
| 117 | rad=6051300. ! rayon de Venus (m) ~6051300 m |
---|
| 118 | daysec= 1.0087e7 ! duree du sol (s) ~1.e7 s |
---|
| 119 | omeg=4.*asin(1.)/19.4141e6 ! vitesse de rotation (rad.s-1) |
---|
| 120 | g= 8.87 ! gravite (m.s-2) ~8.87 |
---|
| 121 | mugaz=43.44 ! Masse molaire de l'atm (g.mol-1) ~43.44 |
---|
| 122 | ! ADAPTATION GCM POUR CP(T) |
---|
| 123 | ! VENUS: Cp(T) = cpp*(T/T0)^nu |
---|
[3540] | 124 | ! avec cpp=1000., T0=460. et nu=0.35 |
---|
| 125 | cpofT=.true. |
---|
[887] | 126 | cpp=1.0e3 |
---|
[3540] | 127 | ! Version Cp constant |
---|
| 128 | ! cpofT=.false. |
---|
| 129 | ! cpp=9.0e2 |
---|
[887] | 130 | r= 8.314511E+0 *1000.E+0/mugaz |
---|
| 131 | rcp= r/cpp |
---|
| 132 | |
---|
[894] | 133 | c----------------------------------------------------------------------- |
---|
[3540] | 134 | c Initialisation des traceursphyven |
---|
[894] | 135 | c --------------------------- |
---|
| 136 | c Choix du nombre de traceurs et du schema pour l'advection |
---|
| 137 | c dans fichier traceur.def |
---|
| 138 | call infotrac_init |
---|
[3519] | 139 | iflag_trac=0 |
---|
[1692] | 140 | if (nqtot.gt.1) iflag_trac=1 |
---|
[894] | 141 | |
---|
| 142 | c Allocation de la tableau q : champs advectes |
---|
| 143 | allocate(q(llm,nqtot)) |
---|
| 144 | allocate(dq(llm,nqtot)) |
---|
| 145 | |
---|
[887] | 146 | c ------------------------------------------------------ |
---|
| 147 | c Lecture des parametres dans "rcm1d.def" |
---|
| 148 | c ------------------------------------------------------ |
---|
| 149 | |
---|
| 150 | c Opening parameters file "rcm1d.def" |
---|
| 151 | c --------------------------------------- |
---|
| 152 | unit =97 |
---|
| 153 | OPEN(unit,file='rcm1d.def',status='old',form='formatted' |
---|
| 154 | . ,iostat=ierr) |
---|
| 155 | |
---|
| 156 | IF(ierr.ne.0) THEN |
---|
| 157 | write(*,*) 'Problem to open "rcm1d.def' |
---|
| 158 | write(*,*) 'Is it there ?' |
---|
| 159 | stop |
---|
| 160 | END IF |
---|
| 161 | |
---|
| 162 | c Date et heure locale du debut du run |
---|
| 163 | c ------------------------------------ |
---|
| 164 | c Date (en sols depuis le solstice de printemps) du debut du run |
---|
| 165 | day0 = 0 |
---|
| 166 | PRINT *,'date de depart ?' |
---|
| 167 | READ(unit,*) day0 |
---|
| 168 | day=REAL(day0) |
---|
| 169 | PRINT *,day0 |
---|
| 170 | c Heure de demarrage |
---|
| 171 | PRINT *,'heure de debut de simulation (entre 0 et 24) ?' |
---|
| 172 | READ(unit,*) time |
---|
| 173 | time=time/24.E+0 |
---|
| 174 | |
---|
| 175 | c Discretisation (Definition de la grille et des pas de temps) |
---|
| 176 | c -------------- |
---|
| 177 | c |
---|
| 178 | nlayer=llm |
---|
| 179 | nlevel=nlayer+1 |
---|
| 180 | nsoil=nsoilmx |
---|
| 181 | PRINT *,'nombre de pas de temps par jour ?' |
---|
| 182 | READ(unit,*) day_step |
---|
| 183 | print*,day_step |
---|
| 184 | |
---|
| 185 | c PRINT *,'nombre d appel au rayonnement par jour ?' |
---|
| 186 | c READ(unit,*) nbapp_rad |
---|
| 187 | c print*,nbapp_rad |
---|
| 188 | c LU DANS PHYSIQ.DEF... |
---|
[2684] | 189 | nbapp_rad = 24000 |
---|
[887] | 190 | |
---|
| 191 | PRINT *,'nombre de jours simules ?' |
---|
[3254] | 192 | READ(unit,*) nb_days |
---|
| 193 | print*,nb_days |
---|
[887] | 194 | |
---|
[3254] | 195 | ndt=nint(nb_days*day_step) |
---|
| 196 | write(*,*) " => will run ", ndt," timesteps" |
---|
[887] | 197 | dtphys=daysec/day_step |
---|
[1121] | 198 | dtime=dtphys |
---|
[887] | 199 | |
---|
| 200 | c Pression de surface sur la planete |
---|
| 201 | c ------------------------------------ |
---|
| 202 | c |
---|
| 203 | PRINT *,'pression au sol' |
---|
| 204 | READ(unit,*) psurf |
---|
| 205 | PRINT *,psurf |
---|
| 206 | c Pression de reference ! voir dyn3d/etat0_venus |
---|
| 207 | c pa = 5.e4 |
---|
| 208 | pa = 1.e6 |
---|
| 209 | preff = 9.2e6 ! 92 bars |
---|
| 210 | c preff = psurf |
---|
| 211 | |
---|
| 212 | c latitude/longitude |
---|
| 213 | c ------------------- |
---|
| 214 | PRINT *,'latitude en degres ?' |
---|
| 215 | READ(unit,*) lati(1) |
---|
| 216 | PRINT *,lati(1) |
---|
[1691] | 217 | lati(1)=lati(1)*pi/180. ! must be in radians. |
---|
[887] | 218 | long(1)=0.E+0 |
---|
| 219 | |
---|
| 220 | c Initialisation speciales "physiq" |
---|
| 221 | c --------------------------------- |
---|
| 222 | |
---|
[1543] | 223 | ! CALL init_phys_lmdz(iim,jjm,llm,1,(/1/)) |
---|
[1517] | 224 | |
---|
[887] | 225 | c la surface de chaque maille est inutile en 1D ---> |
---|
| 226 | area(1)=1.E+0 |
---|
| 227 | c de meme ? |
---|
| 228 | cufi(1)=1.E+0 |
---|
| 229 | cvfi(1)=1.E+0 |
---|
| 230 | |
---|
[1691] | 231 | call ini_cpdet |
---|
| 232 | |
---|
[3540] | 233 | |
---|
[1523] | 234 | c Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid, |
---|
| 235 | c e.g. for cell boundaries, which are meaningless in 1D; so pad these |
---|
| 236 | c with '0.' when necessary |
---|
[1543] | 237 | CALL iniphysiq(1,1,llm, |
---|
| 238 | & 1,comm_lmdz, |
---|
| 239 | & daysec,day0,dtphys, |
---|
[1523] | 240 | & (/lati(1),0./),(/0./), |
---|
| 241 | & (/0.,0./),(/long(1),0./), |
---|
| 242 | & (/ (/area,0./),(/0.,0./) /), |
---|
| 243 | & (/cufi,0.,0.,0./), |
---|
| 244 | & (/cvfi,0./), |
---|
[2684] | 245 | & rad,g,r,cpp,1) |
---|
[1523] | 246 | |
---|
[887] | 247 | c le geopotentiel au sol est inutile en 1D car tout est controle |
---|
| 248 | c par la pression de surface ---> |
---|
| 249 | phisfi(1)=0.E+0 |
---|
| 250 | |
---|
| 251 | c Initialisation pour prendre en compte les vents en 1-D |
---|
| 252 | c ------------------------------------------------------ |
---|
| 253 | |
---|
| 254 | c vent geostrophique |
---|
| 255 | PRINT *,'composante vers l est du vent geostrophique (U) ?' |
---|
| 256 | READ(unit,*) gru |
---|
| 257 | PRINT *,'composante vers le nord du vent geostrophique (V) ?' |
---|
| 258 | READ(unit,*) grv |
---|
| 259 | |
---|
| 260 | c Initialisation des vents au premier pas de temps |
---|
| 261 | DO ilayer=1,nlayer |
---|
| 262 | u(ilayer)=gru |
---|
| 263 | v(ilayer)=grv |
---|
[3254] | 264 | w(ilayer)=0 |
---|
[887] | 265 | ENDDO |
---|
| 266 | |
---|
| 267 | c calcul des pressions et altitudes en utilisant les niveaux sigma |
---|
| 268 | c ---------------------------------------------------------------- |
---|
| 269 | |
---|
| 270 | c Vertical Coordinates (hybrids) |
---|
| 271 | c """""""""""""""""""" |
---|
| 272 | CALL disvert_noterre |
---|
| 273 | |
---|
| 274 | c Calcul au milieu des couches : Vient de la version Mars |
---|
| 275 | c WARNING : le choix de placer le milieu des couches au niveau de |
---|
| 276 | c pression intermédiaire est arbitraire et pourrait etre modifié. |
---|
| 277 | c C'est fait de la meme facon dans disvert |
---|
| 278 | |
---|
| 279 | DO l = 1, llm |
---|
| 280 | aps(l) = 0.5 *( ap(l) +ap(l+1)) |
---|
| 281 | bps(l) = 0.5 *( bp(l) +bp(l+1)) |
---|
| 282 | ENDDO |
---|
| 283 | |
---|
| 284 | DO ilevel=1,nlevel |
---|
| 285 | plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) |
---|
| 286 | ENDDO |
---|
| 287 | |
---|
| 288 | DO ilayer=1,nlayer |
---|
| 289 | play(ilayer)=aps(ilayer)+psurf*bps(ilayer) |
---|
| 290 | pk(ilayer) =cpp*(play(ilayer)/preff)**rcp |
---|
| 291 | c write(120,*) ilayer,plev(ilayer),play(ilayer) |
---|
| 292 | ENDDO |
---|
| 293 | c write(120,*) nlevel,plev(nlevel) |
---|
| 294 | c stop |
---|
| 295 | |
---|
| 296 | pks=cpp*(psurf/preff)**rcp |
---|
| 297 | |
---|
[901] | 298 | c init des variables pour phyredem |
---|
| 299 | c -------------------------------- |
---|
[2624] | 300 | call phys_state_var_init(nqtot) |
---|
[901] | 301 | |
---|
[887] | 302 | c profil de temperature et altitude au premier appel |
---|
| 303 | c -------------------------------------------------- |
---|
| 304 | |
---|
| 305 | c modif par rapport a Mars: |
---|
| 306 | c on envoie dz/T=-log(play/psurf)*r/g dans profile |
---|
| 307 | tmp1(0)=0.0 |
---|
| 308 | tmp1(1)= -log(play(1)/psurf)*r/g |
---|
| 309 | DO ilayer=2,nlayer |
---|
| 310 | tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g |
---|
| 311 | ENDDO |
---|
| 312 | call profile(unit,nlayer+1,tmp1,tmp2) |
---|
| 313 | CLOSE(unit) |
---|
| 314 | |
---|
| 315 | print*," Pression Altitude Temperature" |
---|
| 316 | ilayer=1 |
---|
[901] | 317 | ftsol(1)=tmp2(0) |
---|
[887] | 318 | temp(1)=tmp2(1) |
---|
| 319 | zlay(1)=tmp2(1)*tmp1(1) |
---|
[901] | 320 | print*," 0",ftsol(1) |
---|
[887] | 321 | print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer) |
---|
| 322 | DO ilayer=2,nlayer |
---|
| 323 | temp(ilayer)=tmp2(ilayer) |
---|
| 324 | zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer) |
---|
| 325 | print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer) |
---|
| 326 | ENDDO |
---|
[1301] | 327 | |
---|
| 328 | allocate(tmoy(llm)) |
---|
| 329 | tmoy(:)=temp(:) |
---|
[887] | 330 | |
---|
| 331 | c temperature du sous-sol |
---|
| 332 | c ~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 333 | DO isoil=1,nsoil |
---|
[901] | 334 | ftsoil(1,isoil)=ftsol(1) |
---|
[887] | 335 | ENDDO |
---|
| 336 | |
---|
| 337 | c Initialisation des traceurs |
---|
| 338 | c --------------------------- |
---|
| 339 | |
---|
| 340 | DO iq=1,nqtot |
---|
| 341 | DO ilayer=1,nlayer |
---|
| 342 | q(ilayer,iq) = 0. |
---|
| 343 | ENDDO |
---|
| 344 | ENDDO |
---|
| 345 | |
---|
[3519] | 346 | if (iflag_trac.eq.1) then |
---|
[3254] | 347 | print*,"rcm1d: Loading chemistry profiles from init_1D.txt" |
---|
| 348 | ! check if the file is indeed there |
---|
| 349 | inquire(file="init_1D.txt",exist=file_is_present) |
---|
| 350 | if (file_is_present) then |
---|
| 351 | open(21, form = 'formatted', file = 'init_1D.txt') |
---|
| 352 | read(21,*) |
---|
| 353 | do ilayer = nlayer,1,-1 |
---|
| 354 | read(21,*) idummy, dummy, dummy, (q(ilayer,iq), iq = 1,nqtot) |
---|
[2684] | 355 | ! print*, idummy, q(ilayer,1), q(ilayer,nqtot) |
---|
[3254] | 356 | end do |
---|
| 357 | close(21) |
---|
| 358 | else |
---|
| 359 | write(*,*) "Cannot find input file init_1D.txt!" |
---|
| 360 | write(*,*) "Might as well stop here" |
---|
| 361 | stop |
---|
| 362 | endif ! of if(file_is_present) |
---|
[3519] | 363 | endif ! iflag_trac |
---|
[1305] | 364 | |
---|
[887] | 365 | c Initialisation des parametres d'oro |
---|
| 366 | c ----------------------------------- |
---|
| 367 | |
---|
| 368 | zmea(1) = 0. |
---|
| 369 | zstd(1) = 0. |
---|
| 370 | zsig(1) = 0. |
---|
| 371 | zgam(1) = 0. |
---|
| 372 | zthe(1) = 0. |
---|
| 373 | zpic(1) = 0. |
---|
| 374 | zval(1) = 0. |
---|
| 375 | |
---|
[901] | 376 | c Initialisation albedo |
---|
| 377 | c ---------------------- |
---|
| 378 | |
---|
| 379 | falbe(1)=0.1 |
---|
| 380 | |
---|
[887] | 381 | c Ecriture de "startphy.nc" |
---|
| 382 | c ------------------------- |
---|
| 383 | c (Ce fichier sera aussitot relu au premier |
---|
| 384 | c appel de "physiq", mais il est necessaire pour passer |
---|
| 385 | c les variables purement physiques a "physiq"... |
---|
| 386 | |
---|
| 387 | solsw(1) = 0. |
---|
[901] | 388 | sollw(1) = 0. |
---|
[1691] | 389 | fder(1) = 0. |
---|
[3254] | 390 | dlw(1) = 0. |
---|
| 391 | sollwdown(1)= 0. |
---|
[887] | 392 | radsol(1) = 0. |
---|
[3254] | 393 | |
---|
| 394 | t_ancien(1,:)=0. |
---|
| 395 | q2(1,:)=0. |
---|
| 396 | |
---|
[887] | 397 | radpas = NINT(1.*day_step/nbapp_rad) |
---|
| 398 | soil_model = .true. |
---|
| 399 | |
---|
[901] | 400 | call phyredem("startphy.nc") |
---|
[887] | 401 | |
---|
[901] | 402 | c deallocation des variables phyredem |
---|
| 403 | c ----------------------------------- |
---|
| 404 | call phys_state_var_end |
---|
| 405 | |
---|
[887] | 406 | c======================================================================= |
---|
| 407 | c BOUCLE TEMPORELLE DU MODELE 1D |
---|
| 408 | c======================================================================= |
---|
| 409 | c |
---|
[3519] | 410 | !TEMPORAIRE |
---|
[2684] | 411 | |
---|
[887] | 412 | firstcall=.true. |
---|
| 413 | lastcall=.false. |
---|
| 414 | |
---|
[2684] | 415 | ! debut de boucle temporelle |
---|
| 416 | |
---|
[3254] | 417 | DO idt=1,ndt |
---|
[887] | 418 | IF (idt.eq.ndt) then |
---|
| 419 | lastcall=.true. |
---|
| 420 | c toujours nulle dans le cas de Venus, pour l'instant... |
---|
| 421 | zls = 0.0 |
---|
| 422 | c write(103,*) 'Ls=',zls*180./pi |
---|
| 423 | c write(103,*) 'Lat=', lati(1) |
---|
| 424 | c write(103,*) 'RunEnd - Atmos. Temp. File' |
---|
| 425 | c write(103,*) 'RunEnd - Atmos. Temp. File' |
---|
| 426 | c write(104,*) 'Ls=',zls*180./pi |
---|
| 427 | c write(104,*) 'Lat=', lati(1) |
---|
| 428 | c write(104,*) 'RunEnd - Atmos. Temp. File' |
---|
| 429 | ENDIF |
---|
| 430 | |
---|
| 431 | c calcul du geopotentiel |
---|
| 432 | c ~~~~~~~~~~~~~~~~~~~~~ |
---|
| 433 | ! ADAPTATION GCM POUR CP(T) |
---|
| 434 | DO ilayer=1,nlayer |
---|
| 435 | s(ilayer)=(play(ilayer)/psurf)**rcp |
---|
| 436 | ENDDO |
---|
| 437 | phi(1)=cpp*temp(1)*(1.E+0-s(1)) |
---|
| 438 | DO ilayer=2,nlayer |
---|
| 439 | phi(ilayer)=phi(ilayer-1)+ |
---|
| 440 | & cpp*(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))*0.5 |
---|
| 441 | & *(s(ilayer-1)-s(ilayer)) |
---|
| 442 | |
---|
| 443 | ENDDO |
---|
| 444 | |
---|
| 445 | c appel de la physique |
---|
| 446 | c -------------------- |
---|
| 447 | |
---|
| 448 | CALL physiq (1,llm,nqtot, |
---|
| 449 | , firstcall,lastcall, |
---|
| 450 | , day,time,dtphys, |
---|
| 451 | , plev,play,pk,phi,phisfi, |
---|
| 452 | , presnivs, |
---|
[1621] | 453 | , u,v,temp,q, |
---|
[887] | 454 | , w, |
---|
| 455 | C - sorties |
---|
| 456 | s du,dv,dtemp,dq,dpsurf) |
---|
| 457 | |
---|
[1517] | 458 | c calcul de rho |
---|
| 459 | rho = 0. |
---|
| 460 | c print*,rho |
---|
| 461 | |
---|
| 462 | |
---|
| 463 | c print*,"DT APRES PHYSIQ=",day,time,dtime |
---|
[887] | 464 | c print*,dtemp |
---|
| 465 | c print*,temp |
---|
| 466 | c print*," " |
---|
| 467 | c stop |
---|
| 468 | |
---|
| 469 | c evolution du vent : modele 1D |
---|
| 470 | c ----------------------------- |
---|
| 471 | |
---|
| 472 | c la physique calcule les derivees temporelles de u et v. |
---|
| 473 | c Pas de coriolis |
---|
| 474 | DO ilayer=1,nlayer |
---|
| 475 | du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4 |
---|
| 476 | dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4 |
---|
| 477 | ENDDO |
---|
| 478 | c |
---|
| 479 | c Calcul du temps au pas de temps suivant |
---|
| 480 | c --------------------------------------- |
---|
| 481 | firstcall=.false. |
---|
| 482 | time=time+dtphys/daysec |
---|
| 483 | IF (time.gt.1.E+0) then |
---|
| 484 | time=time-1.E+0 |
---|
| 485 | day=day+1 |
---|
| 486 | ENDIF |
---|
| 487 | |
---|
| 488 | c calcul des vitesses et temperature au pas de temps suivant |
---|
| 489 | c ---------------------------------------------------------- |
---|
| 490 | |
---|
| 491 | DO ilayer=1,nlayer |
---|
| 492 | u(ilayer)=u(ilayer)+dtphys*du(ilayer) |
---|
| 493 | v(ilayer)=v(ilayer)+dtphys*dv(ilayer) |
---|
| 494 | temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer) |
---|
| 495 | ENDDO |
---|
| 496 | |
---|
[1726] | 497 | c calcul des traceurs au pas de temps suivant |
---|
| 498 | c ------------------------------------------- |
---|
| 499 | if (iflag_trac.eq.1) then |
---|
| 500 | DO iq=1,nqtot |
---|
| 501 | DO ilayer=1,nlayer |
---|
| 502 | q(ilayer,iq)=q(ilayer,iq)+dq(ilayer,iq)*dtphys |
---|
| 503 | ENDDO |
---|
| 504 | ENDDO |
---|
| 505 | endif |
---|
| 506 | |
---|
[887] | 507 | c calcul des pressions au pas de temps suivant |
---|
[1726] | 508 | c -------------------------------------------- |
---|
[887] | 509 | |
---|
[1549] | 510 | psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface |
---|
[887] | 511 | DO ilevel=1,nlevel |
---|
| 512 | plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) |
---|
| 513 | ENDDO |
---|
| 514 | DO ilayer=1,nlayer |
---|
| 515 | play(ilayer)=aps(ilayer)+psurf*bps(ilayer) |
---|
[2684] | 516 | ENDDO |
---|
[3519] | 517 | |
---|
[887] | 518 | ENDDO ! fin de la boucle temporelle |
---|
| 519 | |
---|
[3254] | 520 | close(15) |
---|
[887] | 521 | c ======================================================== |
---|
| 522 | c GESTION DES SORTIE |
---|
| 523 | c ======================================================== |
---|
| 524 | |
---|
| 525 | print*,"Temperature finale:" |
---|
| 526 | print*,temp |
---|
| 527 | |
---|
| 528 | c stabilite |
---|
| 529 | DO ilayer=1,nlayer |
---|
| 530 | zlay(ilayer) = phi(ilayer)/g/1000. !en km |
---|
| 531 | ENDDO |
---|
| 532 | DO ilayer=2,nlayer |
---|
| 533 | tmp1(ilayer) = |
---|
| 534 | . (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1)) |
---|
| 535 | . + 1000.*g/cpp |
---|
| 536 | ENDDO |
---|
| 537 | |
---|
| 538 | OPEN(11,file='profile.new') |
---|
| 539 | DO ilayer=1,nlayer |
---|
| 540 | write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer) |
---|
| 541 | ENDDO |
---|
| 542 | |
---|
| 543 | c ======================================================== |
---|
| 544 | END |
---|
| 545 | |
---|
| 546 | c*********************************************************************** |
---|
| 547 | c*********************************************************************** |
---|
| 548 | |
---|
[1403] | 549 | !#include "../dyn3d_common/disvert_noterre.F" |
---|
| 550 | !#include "../dyn3d/abort_gcm.F" |
---|
[887] | 551 | |
---|