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