Changeset 2542 for trunk/LMDZ.GENERIC/libf/phystd/dyn1d
- Timestamp:
- Jul 5, 2021, 2:44:34 PM (3 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf/phystd/dyn1d
- Files:
-
- 2 edited
-
inichim_1D.F90 (modified) (3 diffs)
-
rcm1d.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/dyn1d/inichim_1D.F90
r1882 r2542 1 subroutine inichim_1D(n q, pq, qsurf, ps, &1 subroutine inichim_1D(nlayer, nq, pq, qsurf, play, & 2 2 flagh2o,flagthermo) 3 3 4 use tracer_h 5 USE comvert_mod, ONLY: aps,bps 6 USE mod_grid_phy_lmdz, ONLY: nbp_lev 7 use callkeys_mod 8 use datafile_mod 4 use tracer_h, only: noms, mmol 5 use datafile_mod, only: datadir 9 6 10 7 implicit none … … 26 23 ! Modified 11/2011 Addition of methane Franck Lefevre 27 24 ! Rewritten 04/2012 Franck Lefevre 25 ! Rewritten 03/2021 Yassin Jaziri (Use of #Moderntrac-v1 to init thanks traceur.def) 28 26 ! 29 27 ! Arguments: 30 28 ! ---------- 31 29 ! 32 ! pq(nbp_lon+1,nbp_lat,nbp_lev,nq) Advected fields, ie chemical species here 33 ! qsurf(ngrid,nq) Amount of tracer on the surface (kg/m2) 34 ! ps(nbp_lon+1,nbp_lat) Surface pressure (Pa) 35 ! flagh2o flag for initialisation of h2o (1: yes / 0: no) 36 ! flagthermo flag for initialisation of thermosphere only (1: yes / 0: no) 30 ! nlayer Number of atmospheric layers 31 ! pq(nlayer,nq) Advected fields, ie chemical species here 32 ! qsurf(nq) Amount of tracer on the surface (kg/m2) 33 ! play(nlayer) Pressure (Pa) 34 ! flagh2o flag for initialisation of h2o (1: yes / 0: no) 35 ! flagthermo flag for initialisation of thermosphere only (1: yes / 0: no) 37 36 ! 38 37 !======================================================================= … … 41 40 ! inputs : 42 41 43 integer,intent(in) :: nq ! number of tracers 44 real,intent(in) :: ps ! surface pressure in the gcm (Pa) 45 integer,intent(in) :: flagh2o ! flag for h2o initialisation 46 integer,intent(in) :: flagthermo ! flag for thermosphere initialisation only 42 integer,intent(in) :: nlayer ! Number of atmospheric layers. 43 integer,intent(in) :: nq ! number of tracers 44 real ,intent(in) :: play(nlayer) ! Mid-layer pressure (Pa). 45 integer,intent(in) :: flagh2o ! flag for h2o initialisation 46 integer,intent(in) :: flagthermo ! flag for thermosphere initialisation only 47 47 48 48 ! outputs : 49 49 50 real,intent(out) :: pq(nbp_lev,nq) ! advected fields, ie chemical species51 real,intent(out) :: qsurf(nq)! surface values (kg/m2) of tracers50 real,intent(out) :: pq(nlayer,nq) ! tracers (kg/kg_of_air) 51 real,intent(out) :: qsurf(nq) ! surface values (kg/m2) of tracers 52 52 53 53 ! local : 54 54 55 integer :: iq, l, n, nbqchem 56 integer :: count, ierr, dummy 57 real :: mmean(nbp_lev) ! mean molecular mass (g) 58 real :: pgcm ! pressure at each layer in the gcm (Pa) 59 60 integer, parameter :: nalt = 252 ! number of altitudes in the initialization files 61 integer :: nspe ! number of species in the initialization files 62 integer, allocatable :: niq(:) ! local index of species in initialization files 63 real, dimension(nalt) :: tinit, zzfile ! temperature in initialization files 64 real, dimension(nalt) :: pinit ! pressure in initialization files 65 real, dimension(nalt) :: densinit ! total number density in initialization files 66 real, allocatable :: vmrinit(:,:) ! mixing ratios in initialization files 67 real, allocatable :: vmrint(:) ! mixing ratio interpolated onto the gcm vertical grid 68 real :: vmr 69 70 character(len=20) :: txt ! to store some text 71 logical :: flagnitro ! checks if N species present 72 73 ! 1. identify tracers by their names: (and set corresponding values of mmol) 74 75 ! 1.1 initialize tracer indexes to zero: 76 ! nqmx=nq ! initialize value of nqmx 77 78 ! do iq = 1,nq 79 ! igcm_dustbin(iq) = 0 80 ! end do 81 82 ! igcm_dust_mass = 0 83 ! igcm_dust_number = 0 84 ! igcm_ccn_mass = 0 85 ! igcm_ccn_number = 0 86 igcm_h2o_vap = 0 87 igcm_h2o_ice = 0 88 igcm_co2 = 0 89 igcm_co = 0 90 igcm_o = 0 91 igcm_o1d = 0 92 igcm_o2 = 0 93 igcm_o3 = 0 94 igcm_h = 0 95 igcm_h2 = 0 96 igcm_oh = 0 97 igcm_ho2 = 0 98 igcm_h2o2 = 0 99 igcm_ch4 = 0 100 igcm_n2 = 0 101 igcm_n = 0 102 igcm_n2d = 0 103 igcm_no = 0 104 igcm_no2 = 0 105 igcm_ar = 0 106 igcm_ar_n2 = 0 107 igcm_ch3 = 0 108 igcm_ch = 0 109 igcm_3ch2 = 0 110 igcm_1ch2 = 0 111 igcm_cho = 0 112 igcm_ch2o = 0 113 igcm_c = 0 114 igcm_c2 = 0 115 igcm_c2h = 0 116 igcm_c2h2 = 0 117 igcm_c2h3 = 0 118 igcm_c2h4 = 0 119 igcm_c2h6 = 0 120 igcm_ch2co = 0 121 igcm_ch3co = 0 122 igcm_hcaer = 0 123 124 ! 1.2 find dust tracers 125 count = 0 126 ! 127 ! if (dustbin > 0) then 128 ! do iq = 1,nq 129 ! txt = " " 130 ! write(txt,'(a4,i2.2)') 'dust', count + 1 131 ! if (noms(iq) == txt) then 132 ! count = count + 1 133 ! igcm_dustbin(count) = iq 134 ! mmol(iq) = 100. 135 ! end if 136 ! end do !do iq=1,nq 137 ! end if ! of if (dustbin.gt.0) 138 ! 139 ! if (doubleq) then 140 ! do iq = 1,nq 141 ! if (noms(iq) == "dust_mass") then 142 ! igcm_dust_mass = iq 143 ! count = count + 1 144 ! end if 145 ! if (noms(iq) == "dust_number") then 146 ! igcm_dust_number = iq 147 ! count = count + 1 148 ! end if 149 ! end do 150 ! end if ! of if (doubleq) 151 ! 152 ! if (scavenging) then 153 ! do iq = 1,nq 154 ! if (noms(iq) == "ccn_mass") then 155 ! igcm_ccn_mass = iq 156 ! count = count + 1 157 ! end if 158 ! if (noms(iq) == "ccn_number") then 159 ! igcm_ccn_number = iq 160 ! count = count + 1 161 ! end if 162 ! end do 163 ! end if ! of if (scavenging) 164 ! 165 ! if (submicron) then 166 ! do iq=1,nq 167 ! if (noms(iq) == "dust_submicron") then 168 ! igcm_dust_submicron = iq 169 ! mmol(iq) = 100. 170 ! count = count + 1 171 ! end if 172 ! end do 173 ! end if ! of if (submicron) 174 175 ! 1.3 find chemistry and water tracers 176 177 nbqchem = 0 178 179 do iq = 1,nq 180 if (noms(iq) == "co2") then 181 igcm_co2 = iq 182 mmol(igcm_co2) = 44. 183 count = count + 1 184 nbqchem = nbqchem + 1 185 end if 186 if (noms(iq) == "co") then 187 igcm_co = iq 188 mmol(igcm_co) = 28. 189 count = count + 1 190 nbqchem = nbqchem + 1 191 end if 192 if (noms(iq) == "o") then 193 igcm_o = iq 194 mmol(igcm_o) = 16. 195 count = count + 1 196 nbqchem = nbqchem + 1 197 end if 198 if (noms(iq) == "o1d") then 199 igcm_o1d = iq 200 mmol(igcm_o1d) = 16. 201 count = count + 1 202 nbqchem = nbqchem + 1 203 end if 204 if (noms(iq) == "o2") then 205 igcm_o2 = iq 206 mmol(igcm_o2) = 32. 207 count = count + 1 208 nbqchem = nbqchem + 1 209 end if 210 if (noms(iq) == "o3") then 211 igcm_o3 = iq 212 mmol(igcm_o3) = 48. 213 count = count + 1 214 nbqchem = nbqchem + 1 215 end if 216 if (noms(iq) == "h") then 217 igcm_h = iq 218 mmol(igcm_h) = 1. 219 count = count + 1 220 nbqchem = nbqchem + 1 221 end if 222 if (noms(iq) == "h2") then 223 igcm_h2 = iq 224 mmol(igcm_h2) = 2. 225 count = count + 1 226 nbqchem = nbqchem + 1 227 end if 228 if (noms(iq) == "oh") then 229 igcm_oh = iq 230 mmol(igcm_oh) = 17. 231 count = count + 1 232 nbqchem = nbqchem + 1 233 end if 234 if (noms(iq) == "ho2") then 235 igcm_ho2 = iq 236 mmol(igcm_ho2) = 33. 237 count = count + 1 238 nbqchem = nbqchem + 1 239 end if 240 if (noms(iq) == "h2o2") then 241 igcm_h2o2 = iq 242 mmol(igcm_h2o2) = 34. 243 count = count + 1 244 nbqchem = nbqchem + 1 245 end if 246 if (noms(iq) == "n2") then 247 igcm_n2 = iq 248 mmol(igcm_n2) = 28. 249 count = count + 1 250 nbqchem = nbqchem + 1 251 end if 252 if (noms(iq) == "ch4") then 253 igcm_ch4 = iq 254 mmol(igcm_ch4) = 16. 255 count = count + 1 256 nbqchem = nbqchem + 1 257 end if 258 if (noms(iq) == "ar") then 259 igcm_ar = iq 260 mmol(igcm_ar) = 40. 261 count = count + 1 262 nbqchem = nbqchem + 1 263 end if 264 if (noms(iq) == "n") then 265 igcm_n = iq 266 mmol(igcm_n) = 14. 267 count = count + 1 268 nbqchem = nbqchem + 1 269 end if 270 if (noms(iq) == "n2d") then 271 igcm_n2d = iq 272 mmol(igcm_n2d) = 14. 273 count = count + 1 274 nbqchem = nbqchem + 1 275 end if 276 if (noms(iq) == "no") then 277 igcm_no = iq 278 mmol(igcm_no) = 30. 279 count = count + 1 280 nbqchem = nbqchem + 1 281 end if 282 if (noms(iq) == "no2") then 283 igcm_no2 = iq 284 mmol(igcm_no2) = 46. 285 count = count + 1 286 nbqchem = nbqchem + 1 287 end if 288 289 if (noms(iq) == "h2o_vap") then 290 igcm_h2o_vap = iq 291 mmol(igcm_h2o_vap) = 18. 292 count = count + 1 293 nbqchem = nbqchem + 1 294 end if 295 if (noms(iq) == "h2o_ice") then 296 igcm_h2o_ice = iq 297 mmol(igcm_h2o_ice) = 18. 298 count = count + 1 299 nbqchem = nbqchem + 1 300 end if 301 302 303 if (noms(iq).eq."ch3") then 304 igcm_ch3=iq 305 mmol(igcm_ch3)=15. 306 count=count+1 307 nbqchem = nbqchem + 1 308 endif 309 if (noms(iq).eq."ch") then 310 igcm_ch=iq 311 mmol(igcm_ch)=13. 312 count=count+1 313 nbqchem = nbqchem + 1 314 endif 315 if (noms(iq).eq."3ch2") then 316 igcm_3ch2=iq 317 mmol(igcm_3ch2)=14. 318 count=count+1 319 nbqchem = nbqchem + 1 320 endif 321 if (noms(iq).eq."1ch2") then 322 igcm_1ch2=iq 323 mmol(igcm_1ch2)=14. 324 count=count+1 325 nbqchem = nbqchem + 1 326 endif 327 if (noms(iq).eq."cho") then 328 igcm_cho=iq 329 mmol(igcm_cho)=29. 330 count=count+1 331 nbqchem = nbqchem + 1 332 endif 333 if (noms(iq).eq."ch2o") then 334 igcm_ch2o=iq 335 mmol(igcm_ch2o)=30. 336 count=count+1 337 nbqchem = nbqchem + 1 338 endif 339 if (noms(iq).eq."ch3o") then 340 igcm_ch3o=iq 341 mmol(igcm_ch3o)=31. 342 count=count+1 343 nbqchem = nbqchem + 1 344 endif 345 if (noms(iq).eq."c") then 346 igcm_c=iq 347 mmol(igcm_c)=12. 348 count=count+1 349 nbqchem = nbqchem + 1 350 endif 351 if (noms(iq).eq."c2") then 352 igcm_c2=iq 353 mmol(igcm_c2)=24. 354 count=count+1 355 nbqchem = nbqchem + 1 356 endif 357 if (noms(iq).eq."c2h") then 358 igcm_c2h=iq 359 mmol(igcm_c2h)=25. 360 count=count+1 361 nbqchem = nbqchem + 1 362 endif 363 if (noms(iq).eq."c2h2") then 364 igcm_c2h2=iq 365 mmol(igcm_c2h2)=26. 366 count=count+1 367 nbqchem = nbqchem + 1 368 endif 369 if (noms(iq).eq."c2h3") then 370 igcm_c2h3=iq 371 mmol(igcm_c2h3)=27. 372 count=count+1 373 nbqchem = nbqchem + 1 374 endif 375 if (noms(iq).eq."c2h4") then 376 igcm_c2h4=iq 377 mmol(igcm_c2h4)=28. 378 count=count+1 379 nbqchem = nbqchem + 1 380 endif 381 if (noms(iq).eq."c2h6") then 382 igcm_c2h6=iq 383 mmol(igcm_c2h6)=30. 384 count=count+1 385 nbqchem = nbqchem + 1 386 endif 387 if (noms(iq).eq."ch2co") then 388 igcm_ch2co=iq 389 mmol(igcm_ch2co)=42. 390 count=count+1 391 nbqchem = nbqchem + 1 392 endif 393 if (noms(iq).eq."ch3co") then 394 igcm_ch3co=iq 395 mmol(igcm_ch3co)=43. 396 count=count+1 397 nbqchem = nbqchem + 1 398 endif 399 if (noms(iq).eq."hcaer") then 400 igcm_hcaer=iq 401 mmol(igcm_hcaer)=50. 402 count=count+1 403 nbqchem = nbqchem + 1 404 endif 405 406 407 408 409 ! 1.5 find idealized non-condensible tracer 410 411 if (noms(iq) == "Ar_N2") then 412 igcm_ar_n2 = iq 413 mmol(igcm_ar_n2) = 30. 414 count = count + 1 415 end if 416 417 end do ! of do iq=1,nq 418 419 420 ! 1.6 check that we identified all tracers: 421 422 if (count /= nq) then 423 write(*,*) "inichim_1D: found only ",count," tracers" 424 write(*,*) " expected ",nq 425 do iq = 1,count 426 write(*,*) ' ', iq, ' ', trim(noms(iq)) 427 end do 55 real, allocatable :: pf(:) ! pressure in vmr profile files set in traceur.def 56 real, allocatable :: qf(:) ! vmr in vmr profile files set in traceur.def 57 58 real :: mmean(nlayer) ! mean molecular mass (g) 59 real :: pqx(nlayer,nq) ! tracers (vmr) 60 real :: qx(nq) ! constant vmr set in traceur.def 61 integer :: iq, ilay, iline, nlines, ierr 62 63 CHARACTER(len=100) :: qxf(nq) ! vmr profile files set in traceur.def 64 CHARACTER(len=100) :: fil ! path files 65 character(len=500) :: tracline ! to store a line of text 66 67 logical :: foundback = .false. 68 69 ! 1. initialization 70 71 pq(:,:) = 0. 72 qsurf(:) = 0. 73 pqx(:,:) = 0. 74 qx(:) = 0. 75 qxf(:) = 'None' 76 77 ! 2. load in traceur.def chemistry data for initialization: 78 79 ! Skip nq 80 open(90,file='traceur.def',status='old',form='formatted',iostat=ierr) 81 if (ierr.eq.0) then 82 READ(90,'(A)') tracline 83 IF (trim(tracline).eq.'#ModernTrac-v1') THEN ! Test modern traceur.def 84 DO 85 READ(90,'(A)',iostat=ierr) tracline 86 IF (ierr.eq.0) THEN 87 IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header 88 EXIT 89 ENDIF 90 ELSE ! If pb, or if reached EOF without having found number of tracer 91 write(*,*) "calchim: error reading line of tracers" 92 write(*,*) " (first lines of traceur.def) " 93 stop 94 ENDIF 95 ENDDO 96 ENDIF ! if modern or standard traceur.def 97 else 98 write(*,*) "calchim: error opening traceur.def in inichim_1D" 428 99 stop 429 else430 write(*,*) "inichim_1D: found all expected tracers"431 do iq = 1,nq432 write(*,*) ' ', iq, ' ', trim(noms(iq))433 end do434 end if435 436 ! 1.7 check if nitrogen species are present:437 438 if(igcm_no == 0) then439 !check that no N species is in traceur.def440 if(igcm_n /= 0 .or. igcm_no2 /= 0 .or. igcm_n2d /= 0) then441 write(*,*)'inichim_1D error:'442 write(*,*)'N, NO2 and/or N2D are in traceur.def, but not NO'443 write(*,*)'stop'444 stop445 endif446 flagnitro = .false.447 nspe = 15448 else449 !check that all N species are in traceur.def450 if(igcm_n == 0 .or. igcm_no2 == 0 .or. igcm_n2d == 0) then451 write(*,*)'inichim_1D error:'452 write(*,*)'if NO is in traceur.def, N, NO2 and N2D must also be'453 write(*,*)'stop'454 stop455 endif456 flagnitro = .true.457 nspe = 18458 100 endif 459 101 460 ! 1.8 allocate arrays 461 462 allocate(niq(nspe)) 463 allocate(vmrinit(nalt,nspe)) 464 allocate(vmrint(nspe)) 465 466 ! 2. load in chemistry data for initialization: 467 468 ! order of major species in initialization file: 469 ! 470 ! 1: co2 471 ! 2: ar 472 ! 3: n2 473 ! 4: o2 474 ! 5: co 475 ! 6: o 476 ! 7: h2 477 ! 478 ! order of minor species in initialization file: 479 ! 480 ! 1: h 481 ! 2: oh 482 ! 3: ho2 483 ! 4: h2o 484 ! 5: h2o2 485 ! 6: o1d 486 ! 7: o3 487 ! 488 ! order of nitrogen species in initialization file: 489 ! 490 ! 1: n 491 ! 2: no 492 ! 3: no2 493 ! 4: n2d 494 495 ! major species: 496 497 niq(1) = igcm_co2 498 niq(2) = igcm_ar 499 niq(3) = igcm_n2 500 niq(4) = igcm_o2 501 niq(5) = igcm_co 502 niq(6) = igcm_o 503 niq(7) = igcm_h2 504 505 ! minor species: 506 507 niq(8) = igcm_h 508 niq(9) = igcm_oh 509 niq(10) = igcm_ho2 510 niq(11) = igcm_h2o_vap 511 niq(12) = igcm_h2o2 512 niq(13) = igcm_o1d 513 niq(14) = igcm_o3 514 515 ! nitrogen species: 516 517 if (flagnitro) then 518 niq(15) = igcm_n 519 niq(16) = igcm_no 520 niq(17) = igcm_no2 521 niq(18) = igcm_n2d 522 end if 523 524 ! carbon species: 525 ! niq(18) = igcm_ch4 526 ! niq(19) = igcm_ch3 527 ! niq(20) = igcm_ch 528 ! niq(21) = igcm_1ch2 529 ! niq(22) = igcm_3ch2 530 ! niq(23) = igcm_cho 531 ! niq(24) = igcm_ch2o 532 ! niq(25) = igcm_ch3o 533 ! niq(26) = igcm_c 534 ! niq(27) = igcm_c2 535 ! niq(28) = igcm_c2h 536 ! niq(29) = igcm_c2h2 537 ! niq(30) = igcm_c2h3 538 ! niq(31) = igcm_c2h4 539 ! niq(32) = igcm_c2h6 540 ! niq(33) = igcm_ch2co 541 ! niq(34) = igcm_ch3co 542 ! niq(35) = igcm_hcaer 543 544 545 546 ! 2.1 open initialization files 547 if(1.eq.1) then 548 open(210, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_may.dat') 549 if (ierr /= 0) then 550 write(*,*)'Error : cannot open file atmosfera_LMD_may.dat ' 551 write(*,*)'(in aeronomars/inichim_1D.F)' 552 write(*,*)'It should be in :', trim(datadir),'/' 553 write(*,*)'1) You can change this path in callphys.def with' 554 write(*,*)' datadir=/path/to/datafiles/' 555 write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)' 556 write(*,*)' can be obtained online on:' 557 write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' 102 ! Get data of tracers 103 do iq=1,nq 104 read(90,'(A)') tracline 105 write(*,*)"inichim_1D: iq=",iq,"noms(iq)=",trim(noms(iq)) 106 if (index(tracline,'qx=' ) /= 0) then 107 read(tracline(index(tracline,'qx=')+len('qx='):),*) qx(iq) 108 write(*,*) ' Parameter value (traceur.def) : qx=', qx(iq) 109 else 110 write(*,*) ' Parameter value (default) : qx=', qx(iq) 111 end if 112 if (index(tracline,'qxf=' ) /= 0) then 113 read(tracline(index(tracline,'qxf=')+len('qxf='):),*) qxf(iq) 114 write(*,*) ' Parameter value (traceur.def) : qxf=', qxf(iq) 115 else 116 write(*,*) ' Parameter value (default) : qxf=', qxf(iq) 117 end if 118 end do 119 120 close(90) 121 122 ! 3. initialization of tracers 123 124 ! 3.1 vertical interpolation 125 126 do iq=1,nq 127 if (qx(iq) /= 0.) then 128 pqx(:,iq) = qx(iq) 129 else if (qxf(iq) /= 'None') then 130 ! Opening file 131 fil = trim(datadir)//'/chemical_profiles/'//qxf(iq) 132 print*, 'chemical pofile '//trim(noms(iq))//': ', fil 133 open(UNIT=90,FILE=fil,STATUS='old',iostat=ierr) 134 if (ierr.eq.0) then 135 read(90,*) ! read one header line 136 do ! get number of lines 137 read(90,*,iostat=ierr) 138 if (ierr<0) exit 139 nlines = nlines + 1 140 end do 141 ! allocate reading variable 142 allocate(pf(nlines)) 143 allocate(qf(nlines)) 144 ! read file 145 rewind(90) ! restart from the beggining of the file 146 read(90,*) ! read one header line 147 do iline=1,nlines 148 read(90,*) pf(iline), qf(iline) ! pf [Pa], qf [vmr] 149 end do 150 ! interp in gcm grid 151 do ilay=1,nlayer 152 call intrplf(log(play(ilay)),pqx(ilay,iq),log(pf),qf,nlines) 153 end do 154 ! deallocate for next tracer 155 deallocate(pf) 156 deallocate(qf) 157 else 158 write(*,*) 'inichim_1D: error opening ', fil 159 stop 160 endif 161 close(90) 162 end if 163 end do 164 165 ! 3.2 background gas 166 167 do iq=1,nq 168 if (all(pqx(:,iq)==1.)) then 169 pqx(:,iq) = 0. 170 do ilay=1,nlayer 171 pqx(ilay,iq) = 1-sum(pqx(ilay,:)) 172 if (pqx(ilay,iq)<=0.) then 173 write(*,*) 'inichim_1D: vmr tot > 1 not possible' 174 stop 175 end if 176 end do 177 foundback = .true. 178 exit ! you found the background gas you can skip others 179 end if 180 end do 181 if (.not.foundback) then 182 write(*,*) 'inichim_1D: you need to set a background gas' 183 write(*,*) ' by qx=1. in traceur.def' 558 184 stop 559 185 end if 560 open(220, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_min.dat') 561 if (ierr /= 0) then 562 write(*,*)'Error : cannot open file atmosfera_LMD_min.dat ' 563 write(*,*)'(in aeronomars/inichim_1D.F)' 564 write(*,*)'It should be in :', trim(datadir),'/' 565 write(*,*)'1) You can change this path in callphys.def with' 566 write(*,*)' datadir=/path/to/datafiles/' 567 write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)' 568 write(*,*)' can be obtained online on:' 569 write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' 570 stop 571 end if 572 573 print*,'flagnitro=',flagnitro 574 if(flagnitro) then 575 open(230, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_nitr.dat') 576 if (ierr /= 0) then 577 write(*,*)'Error : cannot open file atmosfera_LMD_nitr.dat ' 578 write(*,*)'(in aeronomars/inichim_1D.F)' 579 write(*,*)'It should be in :', datadir 580 write(*,*)'1) You can change this directory address in ' 581 write(*,*)' file phymars/datafile.h' 582 write(*,*)'2) If necessary atmosfera_LMD_nitr.dat (and others)' 583 write(*,*)' can be obtained online on:' 584 write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' 585 STOP 586 endif 587 endif ! Of if(flagnitro) 588 589 590 ! 2.2 read initialization files' 591 592 ! major species 593 594 read(210,*) 595 do l = 1,nalt 596 read(210,*) dummy, tinit(l), pinit(l), densinit(l), & 597 (vmrinit(l,n), n = 1,7) 598 pinit(l) = pinit(l)*100. ! conversion in Pa 599 pinit(l) = log(pinit(l)) ! for the vertical interpolation 600 end do 601 close(210) 602 603 ! minor species 604 605 read(220,*) 606 do l = 1,nalt 607 read(220,*) dummy, (vmrinit(l,n), n = 8,14) 608 end do 609 close(220) 610 611 ! nitrogen species 612 613 if(flagnitro) then 614 read(230,*) 615 do l = 1,nalt 616 read(230,*) dummy, (vmrinit(l,n), n = 15,18) 186 187 ! 3.3 mmean 188 mmean(:) = 0. 189 do ilay=1,nlayer 190 do iq=1,nq 191 mmean(ilay) = mmean(ilay) + pqx(ilay,iq)*mmol(iq) 617 192 end do 618 close(230) 619 endif 620 endif !if(1.eq.0) 621 622 623 ! initialization for the early eath 624 if (1.eq.0) then 625 do l = 1,nalt 626 vmrinit(l,:)=0e-7 627 vmrinit(l,1)=0.1 !co2 628 vmrinit(l,2)=0.9 !n2 629 vmrinit(l,3)=0.0 !o2 630 vmrinit(l,6)=4e-3 !h2 631 vmrinit(l,10)=1e-9 !h2o 632 vmrinit(l,18)=0.0e-5 !ch4 633 vmrinit(l,10:13)=0.0e-7 !n 634 vmrinit(l,14)=0.0 !n 635 vmrinit(l,15)=0.0 !no 636 vmrinit(l,16)=0.0 !no2 637 ! pinit(l)=aps(l) + bps(l)*ps 638 ! vmrinit(l,18)=2e-3*min(pinit(l)/100.,1.) ! decrease with scale height above 1 hpa 639 vmrinit(l,2)=0.0 640 vmrinit(l,2)=1-sum(vmrinit(l,:)) !n2 641 ! vmrinit(l,4)=0.1 642 ! vmrinit(l,7)=0.001 643 end do 644 endif 645 646 647 648 ! 3. initialization of tracers 649 650 do l = 1,nbp_lev 651 652 pgcm = aps(l) + bps(l)*ps ! gcm pressure 653 pgcm = log(pgcm) ! for the vertical interpolation 654 mmean(l) = 0. 655 656 ! 3.1 vertical interpolation 657 do n = 1,nspe 658 call intrplf(pgcm,vmr,pinit,vmrinit(:,n),nalt) 659 vmrint(n) = vmr 660 iq = niq(n) 661 mmean(l) = mmean(l) + vmrint(n)*mmol(iq) 662 ! mmean(l) = mmean(l) + vmrinit(1,n)*mmol(iq) 663 end do 664 665 ! 3.2 attribute mixing ratio: - all layers or only thermosphere 666 ! - with our without h2o 667 668 if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 0.1)) then 669 do n = 1,nq 670 pq(l,iq) = 0. 671 qsurf(iq) = 0. 672 enddo 673 674 do n = 1,nspe 675 iq = niq(n) 676 ! if (iq /= igcm_h2o_vap .or. flagh2o == 1) then 677 pq(l,iq) = vmrint(n)*mmol(iq)/mmean(l) 678 ! pq(l,iq) = vmrinit(1,n)*mmol(iq)/mmean(1) 679 ! end if 680 end do 681 ! pq(l,igcm_ch4) = 2.e-3*min((aps(l) + bps(l)*ps)/100.,1.)*mmol(igcm_ch4)/mmean(l) 682 end if 683 684 end do 685 686 687 ! set surface values of chemistry tracers to zero 688 689 if (flagthermo == 0) then 690 ! NB: no problem for "surface water vapour" tracer which is always 0 691 do n = 1,nspe 692 iq = niq(n) 693 qsurf(iq) = 0. 193 end do 194 195 ! 3.4 convert vmr to mmr 196 197 do ilay=1,nlayer 198 do iq=1,nq 199 pq(ilay,iq) = pqx(ilay,iq)*mmol(iq)/mmean(ilay) 694 200 end do 695 end if 696 697 ! 3.3 initialization of tracers not contained in the initialization files 698 699 ! methane : 10 ppbv 700 701 ! if (igcm_ch4 /= 0) then 702 ! vmr = 10.e-9 703 ! do l = 1,nbp_lev 704 ! pq(l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(l) 705 ! end do 706 ! ! set surface value to zero 707 ! qsurf(igcm_ch4) = 0. 708 ! end if 709 710 711 ! deallocations 712 713 deallocate(niq) 714 deallocate(vmrinit) 715 deallocate(vmrint) 201 end do 202 203 ! 4. Hard coding 204 ! Do whatever you want here to specify pq and qsurf 205 ! Or use #ModernTrac-v1 and add another option section 2. 716 206 717 207 end -
trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F
r2436 r2542 898 898 allocate(nametmp(nq)) 899 899 nametmp(1:nq)=tname(1:nq) 900 call inichim_1D(n q, q, qsurf, psurf, 0, 0)900 call inichim_1D(nlayer, nq, q, qsurf, play, 0, 0) 901 901 tname(1:nq)=nametmp(1:nq) 902 902 noms(1:nq)=nametmp(1:nq)
Note: See TracChangeset
for help on using the changeset viewer.
