| 1 | SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad, & |
|---|
| 2 | dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc,NLEV) |
|---|
| 3 | |
|---|
| 4 | !-------------------------------------------------------------------- |
|---|
| 5 | ! |
|---|
| 6 | ! Introduction d une routine chimique |
|---|
| 7 | ! |
|---|
| 8 | ! Auteurs: S. Lebonnois, 01/2000 | 09/2003 |
|---|
| 9 | ! adaptation pour Titan 3D: 02/2009 |
|---|
| 10 | ! adaptation pour // : 04/2013 |
|---|
| 11 | ! extension chimie jusqu a 1300km : 10/2013 |
|---|
| 12 | ! |
|---|
| 13 | ! J. Vatant d'Ollone, 02/2017 |
|---|
| 14 | ! + adaptation pour le nouveau titan issu du generic |
|---|
| 15 | ! |
|---|
| 16 | !--------------------------------------------------------------------- |
|---|
| 17 | ! |
|---|
| 18 | |
|---|
| 19 | use dimphy |
|---|
| 20 | |
|---|
| 21 | use datafile_mod, only: datadir |
|---|
| 22 | |
|---|
| 23 | use comcstfi_mod, only: rad, pi, kbol |
|---|
| 24 | use moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat |
|---|
| 25 | use mod_grid_phy_lmdz, only: nbp_lat |
|---|
| 26 | |
|---|
| 27 | implicit none |
|---|
| 28 | |
|---|
| 29 | ! Parameters ( dans common_mod in old titan) |
|---|
| 30 | ! + doivent etre en accord avec titan.h |
|---|
| 31 | ! ------------------------------------------------- |
|---|
| 32 | |
|---|
| 33 | INTEGER,PARAMETER :: NC = 44 ! Nb de composes dans la chimie |
|---|
| 34 | INTEGER,PARAMETER :: ND = 54 ! Nb de photodissociations |
|---|
| 35 | INTEGER,PARAMETER :: NR = 377 ! Nb de reactions dans la chimie |
|---|
| 36 | INTEGER,PARAMETER :: NLRT = 650 ! Pour l'UV 650 niveaux de 2 km |
|---|
| 37 | |
|---|
| 38 | ! Dummy parameters without microphysics |
|---|
| 39 | ! + Just here to keep the whole stuff running without modif C sources |
|---|
| 40 | ! --------------------------------------------------------------------- |
|---|
| 41 | |
|---|
| 42 | INTEGER :: utilaer(16) |
|---|
| 43 | INTEGER :: aerprod = 0 |
|---|
| 44 | INTEGER :: htoh2 = 0 |
|---|
| 45 | |
|---|
| 46 | ! Arguments |
|---|
| 47 | ! --------- |
|---|
| 48 | |
|---|
| 49 | INTEGER nlon ! nb of horiz points |
|---|
| 50 | INTEGER ny ! nb de composes (nqmax-nmicro) |
|---|
| 51 | REAL*8 qy_c(nlon,klev,NC) ! Especes chimiques apres adv.+diss. |
|---|
| 52 | character*10 nomqy_c(NC+1) ! Noms des especes chimiques |
|---|
| 53 | REAL*8 declin_rad,ls_rad ! declinaison et long solaire en radians |
|---|
| 54 | REAL*8 dtchim ! pas de temps chimie |
|---|
| 55 | REAL*8 ctemp(nlon,klev) ! Temperature |
|---|
| 56 | REAL*8 cplay(nlon,klev) ! pression (Pa) |
|---|
| 57 | REAL*8 cplev(nlon,klev+1) ! pression intercouches (Pa) |
|---|
| 58 | REAL*8 czlay(nlon,klev) ! altitude (m) |
|---|
| 59 | REAL*8 czlev(nlon,klev+1) ! altitude intercouches (m) |
|---|
| 60 | |
|---|
| 61 | REAL*8 dqyc(nlon,klev,NC) ! Tendances especes chimiques |
|---|
| 62 | |
|---|
| 63 | INTEGER NLEV ! nbp_lev+70 - a changer si =/ 55 layers ? |
|---|
| 64 | |
|---|
| 65 | ! Local variables : |
|---|
| 66 | ! ----------------- |
|---|
| 67 | |
|---|
| 68 | integer i,j,l,ic,jm1 |
|---|
| 69 | |
|---|
| 70 | ! variables envoyees dans la chimie: double precision |
|---|
| 71 | |
|---|
| 72 | REAL*8 temp_c(NLEV) |
|---|
| 73 | REAL*8 press_c(NLEV) ! T,p(mbar) a 1 lat donnee |
|---|
| 74 | REAL*8 cqy(NLEV,NC) ! frac mol qui seront modifiees |
|---|
| 75 | REAL*8 cqy0(NLEV,NC) ! frac mol avant modif |
|---|
| 76 | REAL*8 surfhaze(NLEV) |
|---|
| 77 | REAL*8 cprodaer(NLEV,4),cmaer(NLEV,4) |
|---|
| 78 | REAL*8 ccsn(NLEV,4),ccsh(NLEV,4) |
|---|
| 79 | ! rmil: milieu de couche, grille pour K,D,p,ct,T,y |
|---|
| 80 | ! rinter: intercouche (grille RA dans la chimie) |
|---|
| 81 | REAL*8 rmil(NLEV),rinter(NLEV),nb(NLEV) |
|---|
| 82 | REAL*8,save,allocatable :: kedd(:) |
|---|
| 83 | |
|---|
| 84 | REAL*8 a,b |
|---|
| 85 | character str1*1,str2*2 |
|---|
| 86 | REAL*8 latit |
|---|
| 87 | character*20 formt1,formt2,formt3 |
|---|
| 88 | |
|---|
| 89 | ! variables locales initialisees au premier appel |
|---|
| 90 | |
|---|
| 91 | LOGICAL firstcall |
|---|
| 92 | DATA firstcall/.true./ |
|---|
| 93 | SAVE firstcall |
|---|
| 94 | |
|---|
| 95 | integer dec,declinint,ialt |
|---|
| 96 | REAL*8 declin_c ! declinaison en degres |
|---|
| 97 | real*8 factalt,factdec,krpddec,krpddecp1,krpddecm1 |
|---|
| 98 | real*8 duree |
|---|
| 99 | REAL*8,save :: mass(NC) |
|---|
| 100 | REAL*8,save,allocatable :: md(:,:) |
|---|
| 101 | REAL*8,save :: botCH4 |
|---|
| 102 | DATA botCH4/0.05/ |
|---|
| 103 | real*8,save :: r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver |
|---|
| 104 | REAL*8,save,allocatable :: krpd(:,:,:,:),krate(:,:) |
|---|
| 105 | integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC) |
|---|
| 106 | integer,save :: prod(200,NC),perte(2,200,NC) |
|---|
| 107 | character dummy*30,fich*7,ficdec*15,curdec*15,name*10 |
|---|
| 108 | real*8 ficalt,oldq,newq,zalt |
|---|
| 109 | logical okfic |
|---|
| 110 | |
|---|
| 111 | |
|---|
| 112 | !----------------------------------------------------------------------- |
|---|
| 113 | !*********************************************************************** |
|---|
| 114 | ! |
|---|
| 115 | ! Initialisations : |
|---|
| 116 | ! ---------------- |
|---|
| 117 | |
|---|
| 118 | duree = dtchim ! passage en real*4 pour appel a routines C |
|---|
| 119 | |
|---|
| 120 | IF (firstcall) THEN |
|---|
| 121 | |
|---|
| 122 | print*,'CHIMIE, premier appel' |
|---|
| 123 | |
|---|
| 124 | ! ************************************ |
|---|
| 125 | ! Au premier appel, initialisation de toutes les variables |
|---|
| 126 | ! pour les routines de la chimie. |
|---|
| 127 | ! ************************************ |
|---|
| 128 | |
|---|
| 129 | allocate(kedd(NLEV)) |
|---|
| 130 | |
|---|
| 131 | allocate(krpd(15,ND+1,NLRT,nbp_lat),krate(NLEV,NR),md(NLEV,NC)) |
|---|
| 132 | |
|---|
| 133 | ! Verification du nombre de composes: coherence common_mod et nqmax-nmicro |
|---|
| 134 | ! ---------------------------------- |
|---|
| 135 | |
|---|
| 136 | if (ny.ne.NC) then |
|---|
| 137 | print*,'PROBLEME de coherence dans le nombre de composes:',ny,NC |
|---|
| 138 | STOP |
|---|
| 139 | endif |
|---|
| 140 | |
|---|
| 141 | ! calcul de temp_c, densites et press_c en moyenne planetaire : |
|---|
| 142 | ! -------------------------------------------------------------- |
|---|
| 143 | |
|---|
| 144 | print*,'pression, densites et temp (init chimie):' |
|---|
| 145 | print*,'level, press_c, nb, temp_c' |
|---|
| 146 | DO l=1,klev |
|---|
| 147 | rinter(l) = (zlevmoy(l)+rad)/1000. |
|---|
| 148 | rmil(l) = (zlaymoy(l)+rad)/1000. |
|---|
| 149 | ! temp_c (K): |
|---|
| 150 | temp_c(l) = tmoy(l) |
|---|
| 151 | ! press_c (mbar): |
|---|
| 152 | press_c(l) = playmoy(l)/100. |
|---|
| 153 | ! nb (cm-3): |
|---|
| 154 | nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l)) |
|---|
| 155 | print*,l,rmil(l)-rad/1000.,press_c(l),nb(l),temp_c(l) |
|---|
| 156 | ENDDO |
|---|
| 157 | rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000. |
|---|
| 158 | |
|---|
| 159 | ! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km. |
|---|
| 160 | do l=klev+2,NLEV |
|---|
| 161 | rinter(l) = rinter(klev+1) & |
|---|
| 162 | + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1) |
|---|
| 163 | rmil(l-1) = (rinter(l-1)+rinter(l))/2. |
|---|
| 164 | enddo |
|---|
| 165 | rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2. |
|---|
| 166 | ! rmil et rinter ne servent que pour la diffusion (> plafond-100km) donc ok en moyenne planetaire |
|---|
| 167 | |
|---|
| 168 | ! lecture de tcp.ver, une seule fois |
|---|
| 169 | ! remplissage r1d,t1d,ct1d,p1d |
|---|
| 170 | open(11,file=TRIM(datadir)//'/tcp.ver',status='old') |
|---|
| 171 | read(11,*) dummy |
|---|
| 172 | do i=1,131 |
|---|
| 173 | read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i) |
|---|
| 174 | !print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i) |
|---|
| 175 | enddo |
|---|
| 176 | close(11) |
|---|
| 177 | |
|---|
| 178 | ! extension pour klev+1 a NLEV avec tcp.ver |
|---|
| 179 | ! la jonction klev/klev+1 est brutale... Tant pis ? |
|---|
| 180 | ialt=1 |
|---|
| 181 | do l=klev+1,NLEV |
|---|
| 182 | zalt=rmil(l)-rad/1000. |
|---|
| 183 | do i=ialt,130 |
|---|
| 184 | if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then |
|---|
| 185 | ialt=i |
|---|
| 186 | endif |
|---|
| 187 | enddo |
|---|
| 188 | factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt)) |
|---|
| 189 | press_c(l) = exp( log(p1d(ialt)) * (1-factalt) & |
|---|
| 190 | + log(p1d(ialt+1)) * factalt ) |
|---|
| 191 | nb(l) = exp( log(ct1d(ialt)) * (1-factalt) & |
|---|
| 192 | + log(ct1d(ialt+1)) * factalt ) |
|---|
| 193 | temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt |
|---|
| 194 | ! print*,l,zalt,press_c(l),nb(l),temp_c(l) |
|---|
| 195 | enddo |
|---|
| 196 | |
|---|
| 197 | ! caracteristiques des composes: |
|---|
| 198 | ! ----------------------------- |
|---|
| 199 | mass(:) = 0.0 |
|---|
| 200 | call comp(nomqy_c,nb,temp_c,mass,md) |
|---|
| 201 | print*,' Mass' |
|---|
| 202 | do ic=1,NC |
|---|
| 203 | print*,nomqy_c(ic),mass(ic) |
|---|
| 204 | ! if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic) |
|---|
| 205 | enddo |
|---|
| 206 | |
|---|
| 207 | ! taux de photodissociations: |
|---|
| 208 | ! -------------------------- |
|---|
| 209 | call disso(krpd,nbp_lat) |
|---|
| 210 | |
|---|
| 211 | ! reactions chimiques: |
|---|
| 212 | ! ------------------- |
|---|
| 213 | call chimie(nomqy_c,nb,temp_c,krate,reactif, & |
|---|
| 214 | nom_perte,nom_prod,perte,prod) |
|---|
| 215 | ! print*,'nom_prod, nom_perte:' |
|---|
| 216 | ! do ic=1,NC |
|---|
| 217 | ! print*,nom_prod(ic),nom_perte(ic) |
|---|
| 218 | ! enddo |
|---|
| 219 | ! print*,'premieres prod, perte(1:reaction,2:compagnon):' |
|---|
| 220 | ! do ic=1,NC |
|---|
| 221 | ! print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic) |
|---|
| 222 | ! enddo |
|---|
| 223 | |
|---|
| 224 | ! l = klev-3 |
|---|
| 225 | ! print*,'krate a p=',press_c(l),' reactifs et produits:' |
|---|
| 226 | ! do ic=1,ND+1 |
|---|
| 227 | ! print*,ic,krpd(7,ic,l,4)*nb(l)," ", |
|---|
| 228 | ! . nomqy_c(reactif(1,ic)+1), |
|---|
| 229 | ! . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1), |
|---|
| 230 | ! . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1) |
|---|
| 231 | ! enddo |
|---|
| 232 | ! do ic=ND+2,NR |
|---|
| 233 | ! print*,ic,krate(l,ic)," ", |
|---|
| 234 | ! . nomqy_c(reactif(1,ic)+1), |
|---|
| 235 | ! . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1), |
|---|
| 236 | ! . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1) |
|---|
| 237 | ! enddo |
|---|
| 238 | |
|---|
| 239 | |
|---|
| 240 | ! init kedd |
|---|
| 241 | ! --------- |
|---|
| 242 | ! kedd stays constant with time and space |
|---|
| 243 | ! => linked to pressure rather than altitude... |
|---|
| 244 | |
|---|
| 245 | kedd(:) = 5e2 ! valeur mise par defaut |
|---|
| 246 | ! UNITE ? doit etre ok pour gptitan |
|---|
| 247 | do l=1,NLEV |
|---|
| 248 | zalt=rmil(l)-rad/1000. ! z en km |
|---|
| 249 | if ((zalt.ge.250.).and.(zalt.le.400.)) then |
|---|
| 250 | kedd(l) = 10.**(3.+(zalt-250.)/50.) |
|---|
| 251 | ! 1E3 at 250 km |
|---|
| 252 | ! 1E6 at 400 km |
|---|
| 253 | elseif ((zalt.gt.400.).and.(zalt.le.600.)) then |
|---|
| 254 | kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.) |
|---|
| 255 | ! 2E7 at 600 km |
|---|
| 256 | elseif ((zalt.gt.600.).and.(zalt.le.900.)) then |
|---|
| 257 | kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.) |
|---|
| 258 | ! 1E8 above 900 km |
|---|
| 259 | elseif ( zalt.gt.900. ) then |
|---|
| 260 | kedd(l) = 1.e8 |
|---|
| 261 | endif |
|---|
| 262 | enddo |
|---|
| 263 | |
|---|
| 264 | ENDIF ! premier appel |
|---|
| 265 | |
|---|
| 266 | !*********************************************************************** |
|---|
| 267 | !----------------------------------------------------------------------- |
|---|
| 268 | |
|---|
| 269 | ! calcul declin_c (en degres) |
|---|
| 270 | ! --------------------------- |
|---|
| 271 | |
|---|
| 272 | declin_c = declin_rad*180./pi |
|---|
| 273 | ! print*,'declinaison en degre=',declin_c |
|---|
| 274 | |
|---|
| 275 | !*********************************************************************** |
|---|
| 276 | !*********************************************************************** |
|---|
| 277 | ! |
|---|
| 278 | ! BOUCLE SUR LES LATITUDES |
|---|
| 279 | ! |
|---|
| 280 | ! * Permet de faire le calcul une seule fois par lat |
|---|
| 281 | ! |
|---|
| 282 | DO j=1,nlon |
|---|
| 283 | |
|---|
| 284 | if (j.eq.1) then |
|---|
| 285 | jm1=1 |
|---|
| 286 | else |
|---|
| 287 | jm1=j-1 |
|---|
| 288 | endif |
|---|
| 289 | |
|---|
| 290 | if((j.eq.1).or.(klat(j).ne.klat(jm1))) then |
|---|
| 291 | |
|---|
| 292 | !*********************************************************************** |
|---|
| 293 | !*********************************************************************** |
|---|
| 294 | |
|---|
| 295 | !----------------------------------------------------------------------- |
|---|
| 296 | ! |
|---|
| 297 | ! Distance radiale (intercouches, en km) |
|---|
| 298 | ! ---------------------------------------- |
|---|
| 299 | |
|---|
| 300 | do l=1,klev |
|---|
| 301 | rinter(l) = (rad+czlev(j,l))/1000. |
|---|
| 302 | rmil(l) = (rad+czlay(j,l))/1000. |
|---|
| 303 | ! print*,rinter(l) |
|---|
| 304 | enddo |
|---|
| 305 | rinter(klev+1)=(rad+czlev(j,klev+1))/1000. |
|---|
| 306 | |
|---|
| 307 | ! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km. |
|---|
| 308 | do l=klev+2,NLEV |
|---|
| 309 | rinter(l) = rinter(klev+1) & |
|---|
| 310 | + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1) |
|---|
| 311 | rmil(l-1) = (rinter(l-1)+rinter(l))/2. |
|---|
| 312 | enddo |
|---|
| 313 | rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2. |
|---|
| 314 | |
|---|
| 315 | !----------------------------------------------------------------------- |
|---|
| 316 | ! |
|---|
| 317 | ! Temperature, pression (mbar), densite (cm-3) |
|---|
| 318 | ! ------------------------------------------- |
|---|
| 319 | |
|---|
| 320 | DO l=1,klev |
|---|
| 321 | ! temp_c (K): |
|---|
| 322 | temp_c(l) = ctemp(j,l) |
|---|
| 323 | ! press_c (mbar): |
|---|
| 324 | press_c(l) = cplay(j,l)/100. |
|---|
| 325 | ! nb (cm-3): |
|---|
| 326 | nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l)) |
|---|
| 327 | ENDDO |
|---|
| 328 | ! extension pour klev+1 a NLEV avec tcp.ver |
|---|
| 329 | ialt=1 |
|---|
| 330 | do l=klev+1,NLEV |
|---|
| 331 | zalt=rmil(l)-rad/1000. |
|---|
| 332 | do i=ialt,130 |
|---|
| 333 | if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then |
|---|
| 334 | ialt=i |
|---|
| 335 | endif |
|---|
| 336 | enddo |
|---|
| 337 | factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt)) |
|---|
| 338 | press_c(l) = exp( log(p1d(ialt)) * (1-factalt) & |
|---|
| 339 | + log(p1d(ialt+1)) * factalt ) |
|---|
| 340 | nb(l) = exp( log(ct1d(ialt)) * (1-factalt) & |
|---|
| 341 | + log(ct1d(ialt+1)) * factalt ) |
|---|
| 342 | temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt |
|---|
| 343 | enddo |
|---|
| 344 | |
|---|
| 345 | !----------------------------------------------------------------------- |
|---|
| 346 | ! |
|---|
| 347 | ! passage krpd => krate |
|---|
| 348 | ! --------------------- |
|---|
| 349 | ! initialisation krate pour dissociations |
|---|
| 350 | |
|---|
| 351 | if ((declin_c*10+267).lt.14.) then |
|---|
| 352 | declinint = 0 |
|---|
| 353 | dec = 0 |
|---|
| 354 | else |
|---|
| 355 | if ((declin_c*10+267).gt.520.) then |
|---|
| 356 | declinint = 14 |
|---|
| 357 | dec = 534 |
|---|
| 358 | else |
|---|
| 359 | declinint = 1 |
|---|
| 360 | dec = 27 |
|---|
| 361 | do while( (declin_c*10+267).ge.real(dec+20) ) |
|---|
| 362 | dec = dec+40 |
|---|
| 363 | declinint = declinint+1 |
|---|
| 364 | enddo |
|---|
| 365 | endif |
|---|
| 366 | endif |
|---|
| 367 | if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then |
|---|
| 368 | factdec = ( declin_c - (dec-267)/10. ) / 4. |
|---|
| 369 | else |
|---|
| 370 | factdec = ( declin_c - (dec-267)/10. ) / 2.7 |
|---|
| 371 | endif |
|---|
| 372 | |
|---|
| 373 | do l=1,NLEV |
|---|
| 374 | |
|---|
| 375 | ! INTERPOL EN ALT POUR k (krpd tous les deux km) |
|---|
| 376 | ialt = int((rmil(l)-rad/1000.)/2.)+1 |
|---|
| 377 | factalt = (rmil(l)-rad/1000.)/2.-(ialt-1) |
|---|
| 378 | |
|---|
| 379 | do i=1,ND+1 |
|---|
| 380 | krpddecm1 = krpd(declinint ,i,ialt ,klat(j)) * (1-factalt) & |
|---|
| 381 | + krpd(declinint ,i,ialt+1,klat(j)) * factalt |
|---|
| 382 | krpddec = krpd(declinint+1,i,ialt ,klat(j)) * (1-factalt) & |
|---|
| 383 | + krpd(declinint+1,i,ialt+1,klat(j)) * factalt |
|---|
| 384 | krpddecp1 = krpd(declinint+2,i,ialt ,klat(j)) * (1-factalt) & |
|---|
| 385 | + krpd(declinint+2,i,ialt+1,klat(j)) * factalt |
|---|
| 386 | |
|---|
| 387 | ! ND+1 correspond a la dissociation de N2 par les GCR |
|---|
| 388 | if ( factdec.lt.0. ) then |
|---|
| 389 | krate(l,i) = krpddecm1 * abs(factdec) & |
|---|
| 390 | + krpddec * ( 1 + factdec) |
|---|
| 391 | endif |
|---|
| 392 | if ( factdec.gt.0. ) then |
|---|
| 393 | krate(l,i) = krpddecp1 * factdec & |
|---|
| 394 | + krpddec * ( 1 - factdec) |
|---|
| 395 | endif |
|---|
| 396 | if ( factdec.eq.0. ) krate(l,i) = krpddec |
|---|
| 397 | enddo |
|---|
| 398 | enddo |
|---|
| 399 | |
|---|
| 400 | !----------------------------------------------------------------------- |
|---|
| 401 | ! |
|---|
| 402 | ! composition |
|---|
| 403 | ! ------------ |
|---|
| 404 | |
|---|
| 405 | do ic=1,NC |
|---|
| 406 | do l=1,klev |
|---|
| 407 | cqy(l,ic) = qy_c(j,l,ic) |
|---|
| 408 | cqy0(l,ic) = cqy(l,ic) |
|---|
| 409 | enddo |
|---|
| 410 | enddo |
|---|
| 411 | |
|---|
| 412 | ! lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV |
|---|
| 413 | |
|---|
| 414 | WRITE(str2,'(i2.2)') klat(j) |
|---|
| 415 | fich = "comp_"//str2//".dat" |
|---|
| 416 | inquire (file=fich,exist=okfic) |
|---|
| 417 | if (okfic) then |
|---|
| 418 | open(11,file=fich,status='old') |
|---|
| 419 | ! premiere ligne=declin |
|---|
| 420 | read(11,'(A15)') ficdec |
|---|
| 421 | write(curdec,'(E15.5)') declin_c |
|---|
| 422 | ! si la declin est la meme: ce fichier a deja ete reecrit |
|---|
| 423 | ! on lit la colonne t-1 au lieu de la colonne t |
|---|
| 424 | ! (cas d une bande de latitude partagee par 2 procs) |
|---|
| 425 | do ic=1,NC |
|---|
| 426 | read(11,'(A10)') name |
|---|
| 427 | if (name.ne.nomqy_c(ic)) then |
|---|
| 428 | print*,"probleme lecture ",fich |
|---|
| 429 | print*,name,nomqy_c(ic) |
|---|
| 430 | stop |
|---|
| 431 | endif |
|---|
| 432 | if (ficdec.eq.curdec) then |
|---|
| 433 | do l=klev+1,NLEV |
|---|
| 434 | read(11,*) ficalt,cqy(l,ic),newq |
|---|
| 435 | enddo |
|---|
| 436 | else |
|---|
| 437 | do l=klev+1,NLEV |
|---|
| 438 | read(11,*) ficalt,oldq,cqy(l,ic) |
|---|
| 439 | enddo |
|---|
| 440 | endif |
|---|
| 441 | enddo |
|---|
| 442 | close(11) |
|---|
| 443 | else ! le fichier n'est pas la |
|---|
| 444 | do ic=1,NC |
|---|
| 445 | do l=klev+1,NLEV |
|---|
| 446 | cqy(l,ic)=cqy(klev,ic) |
|---|
| 447 | enddo |
|---|
| 448 | enddo |
|---|
| 449 | endif |
|---|
| 450 | cqy0 = cqy |
|---|
| 451 | |
|---|
| 452 | !----------------------------------------------------------------------- |
|---|
| 453 | ! |
|---|
| 454 | ! Appel de chimietitan |
|---|
| 455 | ! -------------------- |
|---|
| 456 | |
|---|
| 457 | call gptitan(rinter,temp_c,nb, & |
|---|
| 458 | nomqy_c,cqy, & |
|---|
| 459 | duree,(klat(j)-1),mass,md, & |
|---|
| 460 | kedd,botCH4,krate,reactif, & |
|---|
| 461 | nom_prod,nom_perte,prod,perte, & |
|---|
| 462 | aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh, & |
|---|
| 463 | htoh2,surfhaze) |
|---|
| 464 | |
|---|
| 465 | ! Tendances composition |
|---|
| 466 | ! --------------------- |
|---|
| 467 | |
|---|
| 468 | do ic=1,NC |
|---|
| 469 | do l=1,klev |
|---|
| 470 | dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! en /s |
|---|
| 471 | enddo |
|---|
| 472 | enddo |
|---|
| 473 | |
|---|
| 474 | |
|---|
| 475 | !----------------------------------------------------------------------- |
|---|
| 476 | ! |
|---|
| 477 | ! sauvegarde compo sur NLEV |
|---|
| 478 | ! ------------------------- |
|---|
| 479 | |
|---|
| 480 | ! dans fichier compo_klat(j) (01 à 49) |
|---|
| 481 | |
|---|
| 482 | open(11,file=fich,status='replace') |
|---|
| 483 | ! premiere ligne=declin |
|---|
| 484 | write(11,'(E15.5)') declin_c |
|---|
| 485 | do ic=1,NC |
|---|
| 486 | write(11,'(A10)') nomqy_c(ic) |
|---|
| 487 | do l=klev+1,NLEV |
|---|
| 488 | write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-rad/1000.,cqy0(l,ic),cqy(l,ic) |
|---|
| 489 | enddo |
|---|
| 490 | enddo |
|---|
| 491 | close(11) |
|---|
| 492 | |
|---|
| 493 | !*********************************************************************** |
|---|
| 494 | !*********************************************************************** |
|---|
| 495 | ! FIN: BOUCLE SUR LES LATITUDES |
|---|
| 496 | |
|---|
| 497 | else ! same latitude, we don't do calculations again, only adjust longitudinal variations |
|---|
| 498 | dqyc(j,:,:) = dqyc(jm1,:,:)/qy_c(jm1,:,:)*qy_c(j,:,:) |
|---|
| 499 | endif |
|---|
| 500 | |
|---|
| 501 | ENDDO |
|---|
| 502 | |
|---|
| 503 | !*********************************************************************** |
|---|
| 504 | !*********************************************************************** |
|---|
| 505 | |
|---|
| 506 | firstcall = .false. |
|---|
| 507 | |
|---|
| 508 | end subroutine calchim |
|---|
