| 1 | PROGRAM gcm1d |
|---|
| 2 | IMPLICIT NONE |
|---|
| 3 | |
|---|
| 4 | #include "dimensions.h" |
|---|
| 5 | #include "dimphy.h" |
|---|
| 6 | #include "YOMCST.h" |
|---|
| 7 | #include "comg1d.h" |
|---|
| 8 | #include "clesphys.h" |
|---|
| 9 | #include "control.h" |
|---|
| 10 | |
|---|
| 11 | c Arguments : |
|---|
| 12 | c ----------- |
|---|
| 13 | |
|---|
| 14 | integer ngrid,nlayer,longcles,nqmax |
|---|
| 15 | parameter (longcles=20,nqmax=3) |
|---|
| 16 | cfleur on ne se sert pas de nqmax mais de nqmx defini dnas dimensions.h depend de la |
|---|
| 17 | c compilation |
|---|
| 18 | integer radpas |
|---|
| 19 | integer nday,it,iflag_con,unit,i,itap,n_cooling |
|---|
| 20 | integer iphys_ver,iadv_tvl,i_cvg,i_hum |
|---|
| 21 | |
|---|
| 22 | real fnday |
|---|
| 23 | real h,kappa,z1,z2,sbid,sigbid |
|---|
| 24 | real plev(klev+1),play(klev), psol |
|---|
| 25 | real temp(klev),u(klev),v(klev),tsurf,q(klev,nqmx),w(klev+1) |
|---|
| 26 | real ug(klev),vg(klev) |
|---|
| 27 | real du(klev),dv(klev),dt(klev),dpsrf,dq(klev,nqmx) |
|---|
| 28 | real du_dyn(klev),dv_dyn(klev) |
|---|
| 29 | : ,dt_dyn(klev),dq_dyn(klev,nqmax) |
|---|
| 30 | real phis,presnivs(klev),clesphy0(longcles) |
|---|
| 31 | real time,timestep,ecritphy,day,tho |
|---|
| 32 | real co2_ppm,solaire |
|---|
| 33 | real rlat,rlon,tsol,radsol,psol_f,tsol_f,qsol_f,qsol |
|---|
| 34 | real rugmer,rugsrel,snow,agesno,deltat,zmea,zstd,zsig,sn |
|---|
| 35 | real zgam,zthe,zpic,zval |
|---|
| 36 | real ema_sig(klev),ema_w(klev),ncst_cbmf,ema_cbmfz,ema_pcb |
|---|
| 37 | : ,zz_f(klev),vu_f(klev),vv_f(klev),t_f(klev),q_f(klev,nqmax) |
|---|
| 38 | |
|---|
| 39 | real temp0(klev),q0(klev,3) |
|---|
| 40 | |
|---|
| 41 | real dt_cooling(klev),dq_cooling(klev) |
|---|
| 42 | real d_t_cool(klev),d_q_cool(klev) |
|---|
| 43 | real d_t_adv(klev),d_q_adv(klev,nqmax) |
|---|
| 44 | real d_t_cvg(klev),d_q_cvg(klev) |
|---|
| 45 | real ht(100),hq(100),hw(100) |
|---|
| 46 | |
|---|
| 47 | real phy_nat(360) |
|---|
| 48 | real phy_alb(360) |
|---|
| 49 | real phy_sst(360) |
|---|
| 50 | real phy_bil(360) |
|---|
| 51 | real phy_rug(360) |
|---|
| 52 | real phy_ice(360) |
|---|
| 53 | |
|---|
| 54 | logical cycle_diurne,soil_model,new_oliq,ok_orodr |
|---|
| 55 | : ,ok_orolf,ok_limitvrai |
|---|
| 56 | |
|---|
| 57 | logical firstcall,lastcall,flag_cool |
|---|
| 58 | logical itsourcecont |
|---|
| 59 | c itsourcecont permet de choisir entre les sources pour plus de 5 traceurs |
|---|
| 60 | character*80 ans,file_forctl, file_fordat, file_start,file |
|---|
| 61 | character*2 cnbl |
|---|
| 62 | |
|---|
| 63 | c----------------------------------------------------------------------- |
|---|
| 64 | |
|---|
| 65 | COMMON/comvert/ |
|---|
| 66 | * s(llm),sig(llm+1),ds(llm),dsig(llm),dsig1(llm),sdsig(llm) |
|---|
| 67 | . ,sig_s(llm) |
|---|
| 68 | |
|---|
| 69 | REAL s,sig,ds,dsig,dsig1,sdsig,sig_s |
|---|
| 70 | |
|---|
| 71 | c----------------------------------------------------------------------- |
|---|
| 72 | c INCLUDE 'temps.h': |
|---|
| 73 | |
|---|
| 74 | COMMON/temps/itaufin,dtd, |
|---|
| 75 | s day_ini,day_end,anne_ini |
|---|
| 76 | |
|---|
| 77 | INTEGER itaufin |
|---|
| 78 | INTEGER*4 day_ini,day_end,anne_ini |
|---|
| 79 | REAL dtd |
|---|
| 80 | |
|---|
| 81 | c----------------------------------------------------------------------- |
|---|
| 82 | c dynamical tendencies |
|---|
| 83 | |
|---|
| 84 | INTEGER l,ierr,aslun,nlevel,iq,ll |
|---|
| 85 | |
|---|
| 86 | REAL longitude,latitude |
|---|
| 87 | REAL zlay(klev),phi(klev) |
|---|
| 88 | REAL paire |
|---|
| 89 | |
|---|
| 90 | DATA latitude,longitude/0.,0./ |
|---|
| 91 | |
|---|
| 92 | c----------------------------------------------------------------------- |
|---|
| 93 | c Initialisations des constantes |
|---|
| 94 | c ------------------------------- |
|---|
| 95 | |
|---|
| 96 | c constantes |
|---|
| 97 | c ---------- |
|---|
| 98 | |
|---|
| 99 | time=0. |
|---|
| 100 | tho=3600. |
|---|
| 101 | it=0 |
|---|
| 102 | |
|---|
| 103 | call suphec |
|---|
| 104 | |
|---|
| 105 | c parametres lus dans execution_1D: |
|---|
| 106 | c --------------------------------- |
|---|
| 107 | |
|---|
| 108 | read(*,*) fnday |
|---|
| 109 | print*, 'fnday',fnday |
|---|
| 110 | read(*,*) ecritphy |
|---|
| 111 | print*, 'ecritphy',ecritphy |
|---|
| 112 | read(*,*) timestep |
|---|
| 113 | print*, 'timestep',timestep |
|---|
| 114 | read(*,*) iflag_con |
|---|
| 115 | print*, 'iflag_con',iflag_con |
|---|
| 116 | read(*,*) flag_cool |
|---|
| 117 | print*, 'flag_cool',flag_cool |
|---|
| 118 | read(*,'(a)') cnbl |
|---|
| 119 | print*, 'cnbl',cnbl |
|---|
| 120 | read(*,*) radpas |
|---|
| 121 | print*, 'radpas',radpas |
|---|
| 122 | read(*,*) anne_ini ! sans importance mais il faut qqchose |
|---|
| 123 | print*, 'anne_ini',anne_ini |
|---|
| 124 | |
|---|
| 125 | c si ecritphy=0: on ecrit la physique a chaque pas de temps: |
|---|
| 126 | |
|---|
| 127 | if (ecritphy.eq.0.) then |
|---|
| 128 | ecritphy = timestep/rday |
|---|
| 129 | endif |
|---|
| 130 | write(*,*) 'ECRITPHY = ',ecritphy |
|---|
| 131 | |
|---|
| 132 | firstcall=.true. |
|---|
| 133 | lastcall=.false. |
|---|
| 134 | day=1. |
|---|
| 135 | itsourcecont=.true. |
|---|
| 136 | |
|---|
| 137 | open(78,file='transil.dat',form='unformatted', |
|---|
| 138 | s access='direct',recl=4*llm*llm) |
|---|
| 139 | |
|---|
| 140 | ngrid=1 |
|---|
| 141 | IF (ngrid.NE.klon) THEN |
|---|
| 142 | PRINT*,'STOP in inifis' |
|---|
| 143 | PRINT*,'Probleme de dimenesions :' |
|---|
| 144 | PRINT*,'ngrid = ',ngrid |
|---|
| 145 | PRINT*,'klon = ',klon |
|---|
| 146 | STOP |
|---|
| 147 | ENDIF |
|---|
| 148 | nlayer=klev |
|---|
| 149 | nlevel=nlayer+1 |
|---|
| 150 | |
|---|
| 151 | c ------------------------------------------------------------------- |
|---|
| 152 | c Initialisation de la discretisation verticale |
|---|
| 153 | c --------------------------------------------- |
|---|
| 154 | c sb CALL disvert(llm,rkappa,sig,dsig,s,ds,dsig1,sdsig) |
|---|
| 155 | |
|---|
| 156 | c ------------------------------------------------------------------- |
|---|
| 157 | c Profils initiaux. |
|---|
| 158 | c ----------------- |
|---|
| 159 | |
|---|
| 160 | unit = 80 |
|---|
| 161 | open(unit,file='start'//cnbl//'.data',form='formatted') |
|---|
| 162 | read(unit,*) |
|---|
| 163 | read(unit,*) |
|---|
| 164 | read(unit,*) (play(l),l=1,nlayer) |
|---|
| 165 | read(unit,*) (plev(l),l=1,nlevel) |
|---|
| 166 | read(unit,*) (temp(l),l=1,nlayer) |
|---|
| 167 | read(unit,*) (q(l,1),l=1,nlayer) |
|---|
| 168 | read(unit,*) (q(l,2),l=1,nlayer) |
|---|
| 169 | read(unit,*) (ema_sig(l),l=1,nlayer) |
|---|
| 170 | read(unit,*) (ema_w(l),l=1,nlayer) ! reinitialise + loin??? |
|---|
| 171 | read(unit,*) ncst_cbmf, ema_cbmfz, ema_pcb |
|---|
| 172 | read(unit,*) (zz_f(l),l=1,nlayer) |
|---|
| 173 | read(unit,*) (vu_f(l),l=1,nlayer) |
|---|
| 174 | read(unit,*) (vv_f(l),l=1,nlayer) |
|---|
| 175 | close(unit) |
|---|
| 176 | write(*,*) 'Lecture fichier start.data ok' |
|---|
| 177 | |
|---|
| 178 | DO l = 1, nlayer |
|---|
| 179 | u(l) = vu_f(l) |
|---|
| 180 | v(l) = vv_f(l) |
|---|
| 181 | c################## |
|---|
| 182 | cfleur ATTENTION q(3) ce nest plus la glace mais le radon si on active phytrac |
|---|
| 183 | c##################3 |
|---|
| 184 | q(l,3) = 0. ! on initialise la glace a zero |
|---|
| 185 | PRINT*,plev(l),play(l) |
|---|
| 186 | c on recopie le profil initial dans les champs de forcing: |
|---|
| 187 | t_f(l) = temp(l) |
|---|
| 188 | q_f(l,1) = q(l,1) |
|---|
| 189 | q_f(l,2) = q(l,2) |
|---|
| 190 | q_f(l,3) = q(l,3) |
|---|
| 191 | ENDDO |
|---|
| 192 | |
|---|
| 193 | |
|---|
| 194 | c Lecture/creation des conditions aux limites: |
|---|
| 195 | c -------------------------------------------- |
|---|
| 196 | |
|---|
| 197 | unit = 81 |
|---|
| 198 | open(unit,file='startphy.data',form='formatted') |
|---|
| 199 | read(unit,*) |
|---|
| 200 | read(unit,*) co2_ppm |
|---|
| 201 | read(unit,*) solaire |
|---|
| 202 | read(unit,*) rlat |
|---|
| 203 | read(unit,*) tsol |
|---|
| 204 | read(unit,*) radsol |
|---|
| 205 | close(unit) |
|---|
| 206 | write(*,*) 'Lecture fichier startphy.data ok' |
|---|
| 207 | |
|---|
| 208 | unit = 82 |
|---|
| 209 | open(unit,file='condsurf.data',form='formatted') |
|---|
| 210 | read(unit,*) psol_f |
|---|
| 211 | read(unit,*) tsol_f |
|---|
| 212 | read(unit,*) qsol_f |
|---|
| 213 | close(unit) |
|---|
| 214 | write(*,*) 'Lecture fichier condsurf.data ok' |
|---|
| 215 | |
|---|
| 216 | if (play(1).lt.10000.) then |
|---|
| 217 | do l = 1, nlayer |
|---|
| 218 | play(l) = play(l)*100. |
|---|
| 219 | enddo |
|---|
| 220 | endif |
|---|
| 221 | if (plev(1).lt.10000.) then |
|---|
| 222 | do l = 1, nlevel |
|---|
| 223 | plev(l) = plev(l)*100. |
|---|
| 224 | enddo |
|---|
| 225 | endif |
|---|
| 226 | |
|---|
| 227 | if (abs(psol_f-plev(1)) .gt. 1.) then |
|---|
| 228 | print *,' Incompatibilite entre psol et profil' |
|---|
| 229 | : ,' de pression' |
|---|
| 230 | print *,' psol = ',psol_f,' plev(1) = ',plev(1) |
|---|
| 231 | stop |
|---|
| 232 | endif |
|---|
| 233 | |
|---|
| 234 | tsol = tsol_f ! temp au sol prise dans condsurf.data |
|---|
| 235 | psol = psol_f ! pression au sol prise dans condsurf.data |
|---|
| 236 | qsol = qsol_f |
|---|
| 237 | rugmer = 0.0001 ! valeur de cchlim.data |
|---|
| 238 | rugsrel = 0.0 ! (rugsrel = rugoro) |
|---|
| 239 | snow = 0.0 |
|---|
| 240 | sn=0. |
|---|
| 241 | agesno = 50.0 |
|---|
| 242 | rlon = 0.0 |
|---|
| 243 | deltat = 0.0 ! ne sert que pour les slab_ocean |
|---|
| 244 | phis = 0. |
|---|
| 245 | zmea = 0. |
|---|
| 246 | zstd = 0. |
|---|
| 247 | zsig = 0. |
|---|
| 248 | zgam = 0. |
|---|
| 249 | zthe = 0. |
|---|
| 250 | zpic = 0. |
|---|
| 251 | zval = 0. |
|---|
| 252 | |
|---|
| 253 | do i=1,360 |
|---|
| 254 | phy_sst(i) = tsol |
|---|
| 255 | phy_nat(i) = 0.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise |
|---|
| 256 | phy_alb(i) = 0.15 ! albedo land only (old value condsurf_jyg=0.3) |
|---|
| 257 | phy_bil(i) = 1.0 ! ne sert que pour les slab_ocean |
|---|
| 258 | phy_rug(i) = 0.1 ! longueur rugosite utilisee sur land only |
|---|
| 259 | phy_ice(i) = 0.0 ! fraction de glace (?) |
|---|
| 260 | enddo |
|---|
| 261 | |
|---|
| 262 | call writelim |
|---|
| 263 | s (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice) |
|---|
| 264 | |
|---|
| 265 | |
|---|
| 266 | |
|---|
| 267 | c controles du run (en 3D: lus dans run.def) : |
|---|
| 268 | |
|---|
| 269 | cycle_diurne = .FALSE. |
|---|
| 270 | soil_model = .FALSE. |
|---|
| 271 | new_oliq = .FALSE. |
|---|
| 272 | ok_orodr = .FALSE. |
|---|
| 273 | ok_orolf = .FALSE. |
|---|
| 274 | ok_limitvrai = .FALSE. |
|---|
| 275 | |
|---|
| 276 | do i = 1, longcles |
|---|
| 277 | clesphy0(i) = 0. |
|---|
| 278 | enddo |
|---|
| 279 | |
|---|
| 280 | nbapp_rad = NINT(86400./radpas/timestep) |
|---|
| 281 | |
|---|
| 282 | clesphy0(1) = FLOAT( iflag_con ) |
|---|
| 283 | clesphy0(2) = FLOAT( nbapp_rad ) |
|---|
| 284 | if (cycle_diurne ) clesphy0(3) = 1. |
|---|
| 285 | if ( soil_model ) clesphy0(4) = 1. |
|---|
| 286 | if ( new_oliq ) clesphy0(5) = 1. |
|---|
| 287 | if ( ok_orodr ) clesphy0(6) = 1. |
|---|
| 288 | if ( ok_orolf ) clesphy0(7) = 1. |
|---|
| 289 | if ( ok_limitvrai) clesphy0(8) = 1. |
|---|
| 290 | |
|---|
| 291 | |
|---|
| 292 | c ------------------------------------------------------------------- |
|---|
| 293 | c Discretisation verticale: |
|---|
| 294 | c ------------------------- |
|---|
| 295 | |
|---|
| 296 | c calcul sig,dsig,s,ds,dsig1,sdsig a partir des play, plev lus: |
|---|
| 297 | |
|---|
| 298 | do l = 1, nlevel |
|---|
| 299 | sig(l)= plev(l)/plev(1) |
|---|
| 300 | enddo |
|---|
| 301 | |
|---|
| 302 | do l = 1, nlayer |
|---|
| 303 | s(l)= sig(l)**RKAPPA |
|---|
| 304 | enddo |
|---|
| 305 | |
|---|
| 306 | do l = 2, nlayer |
|---|
| 307 | ds(l) = s(l-1) - s(l) |
|---|
| 308 | enddo |
|---|
| 309 | ds(1) = 1. - s(1) |
|---|
| 310 | |
|---|
| 311 | do l = 1, nlayer |
|---|
| 312 | dsig(l) = sig(l)-sig(l+1) |
|---|
| 313 | sdsig(l) = s(l) * dsig(l) |
|---|
| 314 | dsig1(l) = 1./dsig(l) |
|---|
| 315 | presnivs(l) = play(l)/100./100. |
|---|
| 316 | enddo |
|---|
| 317 | |
|---|
| 318 | print*,'Diagnostique de la discretisation verticale' |
|---|
| 319 | |
|---|
| 320 | print* |
|---|
| 321 | h=7. |
|---|
| 322 | kappa = RKAPPA |
|---|
| 323 | print*,'comparaison de sig(l) et (s(l)+s(l+1))/2)**(1/K)' |
|---|
| 324 | do 14 l=1,llm-1 |
|---|
| 325 | sigbid=(0.5*(s(l)+s(l+1)))**(1./kappa) |
|---|
| 326 | print*,'sig(',l+1,') = ',sig(l+1), |
|---|
| 327 | : ' valeur approchee :',sigbid,' ',dsig(l) |
|---|
| 328 | 14 continue |
|---|
| 329 | print* |
|---|
| 330 | print*,'comparaison de s(l) et (sig(l)+sig(l+1))/2)**K' |
|---|
| 331 | do 15 l=1,llm |
|---|
| 332 | sbid=(0.5*(sig(l+1)+sig(l)))**kappa |
|---|
| 333 | print*,' s(',l,') = ',s(l), |
|---|
| 334 | : ' valeur approchee :',sbid |
|---|
| 335 | 15 continue |
|---|
| 336 | |
|---|
| 337 | print*,'Altitude approchee z,dz' |
|---|
| 338 | print* |
|---|
| 339 | z1=0. |
|---|
| 340 | print*,' l Z DZ Ztop dsig' |
|---|
| 341 | DO 18 l=1,llm-1 |
|---|
| 342 | z2=-h*log(sig(l+1)) |
|---|
| 343 | write(*,'(i5,3x,4f8.4)') l,-h*log(s(l))/kappa,z2-z1,z2 |
|---|
| 344 | : ,dsig(l) |
|---|
| 345 | write(14,'(3x,i5,1f10.4)') l,-h*log(s(l))/kappa |
|---|
| 346 | z1=z2 |
|---|
| 347 | 18 CONTINUE |
|---|
| 348 | write(*,'(i5,3x,3f8.4)') l,-h*log(s(llm))/kappa |
|---|
| 349 | write(14,'(3x,i5,1f10.4)') l,-h*log(s(llm))/kappa |
|---|
| 350 | |
|---|
| 351 | c print*,'Correspondance approx pression-altitude:' |
|---|
| 352 | c z1 = 0. |
|---|
| 353 | c do l = 1, llm |
|---|
| 354 | c z1 = z1 + 287.*temp(l)/9.81/play(l)*(plev(l)-plev(l+1)) |
|---|
| 355 | c write(*,*) l,play(l),z1 |
|---|
| 356 | c enddo |
|---|
| 357 | |
|---|
| 358 | c ------------------------------------------------------------------- |
|---|
| 359 | c Ecriture de l'etat initial: |
|---|
| 360 | c --------------------------- |
|---|
| 361 | |
|---|
| 362 | c call ini_fis(timestep,radpas,iflag_con,anne_ini |
|---|
| 363 | c : ,co2_ppm,solaire,rlat,rlon,tsol,deltat |
|---|
| 364 | c : ,qsol,snow,radsol,rugmer,agesno,zmea,zstd,zsig,zgam |
|---|
| 365 | c : ,zthe,zpic,zval,rugsrel |
|---|
| 366 | c : ,phy_sst,phy_nat,phy_alb,phy_bil,phy_rug,phy_ice) |
|---|
| 367 | |
|---|
| 368 | call physdem(rlon, rlat, timestep,radpas,co2_ppm, |
|---|
| 369 | . solaire,tsol, qsol, |
|---|
| 370 | . sn, radsol, deltat, rugmer, |
|---|
| 371 | . agesno, zmea, zstd, zsig, |
|---|
| 372 | . zgam, zthe, zpic, zval, |
|---|
| 373 | . rugsrel) |
|---|
| 374 | |
|---|
| 375 | |
|---|
| 376 | |
|---|
| 377 | c ------------------------------------------------------------------- |
|---|
| 378 | c Options de la simulation: |
|---|
| 379 | c ------------------------- |
|---|
| 380 | |
|---|
| 381 | unit = 83 |
|---|
| 382 | open(unit,file='version.data',form='formatted') |
|---|
| 383 | read(unit,*) iphys_ver,iadv_tvl,i_cvg,i_hum |
|---|
| 384 | close(unit) |
|---|
| 385 | write(*,*) 'Lecture fichier version.data ok' |
|---|
| 386 | |
|---|
| 387 | write(*,*) 'physiq # ',iphys_ver,' (0=std; 1=forced glob cis; |
|---|
| 388 | : 2=forced loc cis)' |
|---|
| 389 | write(*,*) 'advection switch : ',iadv_tvl,' (0=no adv ; |
|---|
| 390 | : 1=forced adv)' |
|---|
| 391 | write(*,*) 'convergence switch : ',i_cvg,' (0=no conv; |
|---|
| 392 | : 1=forced conv)' |
|---|
| 393 | |
|---|
| 394 | |
|---|
| 395 | |
|---|
| 396 | c Eventuellement, preparation de la "bulle froide": |
|---|
| 397 | c ------------------------------------------------- |
|---|
| 398 | |
|---|
| 399 | IF (flag_cool) THEN |
|---|
| 400 | unit = 84 |
|---|
| 401 | open(unit,file='cool_buble.data',form='formatted') |
|---|
| 402 | read(unit,*) n_cooling |
|---|
| 403 | do l = 1, nlayer |
|---|
| 404 | read(unit,*) dt_cooling(l),dq_cooling(l) |
|---|
| 405 | enddo |
|---|
| 406 | close(unit) |
|---|
| 407 | write(*,*) 'Lecture fichier cool_buble.data ok' |
|---|
| 408 | ENDIF |
|---|
| 409 | |
|---|
| 410 | |
|---|
| 411 | c Eventuellement, preparation du forcage par la convergence: |
|---|
| 412 | c ---------------------------------------------------------- |
|---|
| 413 | |
|---|
| 414 | IF (i_cvg .EQ. 1) THEN |
|---|
| 415 | file_forctl = 'forcing.ctl' |
|---|
| 416 | file_fordat = 'forcing.dat' |
|---|
| 417 | file_start = 'start'//cnbl//'.data' |
|---|
| 418 | |
|---|
| 419 | call copie(klev,play,psol,file_forctl) |
|---|
| 420 | call get_uvd2(itap,file_forctl,file_fordat,file_start |
|---|
| 421 | : ,ht,hq,hw) |
|---|
| 422 | ENDIF |
|---|
| 423 | |
|---|
| 424 | |
|---|
| 425 | c----------------------------------------------------------------------- |
|---|
| 426 | c initialisation pour GRADS-1D |
|---|
| 427 | c ---------------------------- |
|---|
| 428 | |
|---|
| 429 | g1d_nlayer=nlayer |
|---|
| 430 | g1d_nomfich='grads1d.dat' |
|---|
| 431 | g1d_unitfich=30 |
|---|
| 432 | g1d_nomctl='grads1d.ctl' |
|---|
| 433 | g1d_unitctl=31 |
|---|
| 434 | g1d_premier=.true. |
|---|
| 435 | |
|---|
| 436 | file='sort' |
|---|
| 437 | call inigrads(1,1 |
|---|
| 438 | s ,0.,1.,-2.,2.,1,0.,-2.,2.,1. |
|---|
| 439 | s ,llm,presnivs,1000. |
|---|
| 440 | s ,1800.,file,'Diagconvect') |
|---|
| 441 | print*,'Fin de Initialisation de wrgras' |
|---|
| 442 | |
|---|
| 443 | |
|---|
| 444 | c======================================================================= |
|---|
| 445 | c DEBUT DE L'INTEGRATION TEMPORELLE: |
|---|
| 446 | c ================================== |
|---|
| 447 | |
|---|
| 448 | itap = 1 |
|---|
| 449 | |
|---|
| 450 | 1 continue |
|---|
| 451 | |
|---|
| 452 | c calcul du geopotentiel: |
|---|
| 453 | c ----------------------- |
|---|
| 454 | |
|---|
| 455 | phi(1)=RD*temp(1)* |
|---|
| 456 | :(plev(1)-play(1))/(.5*(plev(1)+play(1))) |
|---|
| 457 | do l = 1, nlayer-1 |
|---|
| 458 | phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))* |
|---|
| 459 | : (play(l)-play(l+1))/(play(l)+play(l+1)) |
|---|
| 460 | enddo |
|---|
| 461 | |
|---|
| 462 | c PRINT*,'altitude km et T (K)' |
|---|
| 463 | c WRITE(6,'(2f10.2)') (.001*phi(l)/g,temp(l),l=1,nlayer) |
|---|
| 464 | |
|---|
| 465 | |
|---|
| 466 | c pour etre coherent avec LMDZ.3, on passe 2 arguments |
|---|
| 467 | c en plus a physiq: phis et aire (pas utilises en 1D): |
|---|
| 468 | paire = 1. ! aire de la maille |
|---|
| 469 | |
|---|
| 470 | |
|---|
| 471 | c eventuellement, induction de la convection par chauffage du |
|---|
| 472 | c bas et refroidissement du haut de la couche limite: |
|---|
| 473 | c ------------------------------------------------------------ |
|---|
| 474 | |
|---|
| 475 | IF (flag_cool) THEN |
|---|
| 476 | |
|---|
| 477 | call cool_pool(itap |
|---|
| 478 | e ,n_cooling,dt_cooling,dq_cooling |
|---|
| 479 | s ,d_t_cool,d_q_cool) |
|---|
| 480 | |
|---|
| 481 | do l = 1, nlayer |
|---|
| 482 | temp(l) = temp(l) + d_t_cool(l) |
|---|
| 483 | q(l,1) = q(l,1) + d_q_cool(l) |
|---|
| 484 | enddo |
|---|
| 485 | |
|---|
| 486 | ELSE |
|---|
| 487 | |
|---|
| 488 | do l = 1, nlayer |
|---|
| 489 | d_t_cool(l) = 0. |
|---|
| 490 | d_q_cool(l) = 0. |
|---|
| 491 | enddo |
|---|
| 492 | |
|---|
| 493 | ENDIF |
|---|
| 494 | |
|---|
| 495 | |
|---|
| 496 | c eventuellement, ajouter une tendance relative a une |
|---|
| 497 | c translation du domaine (ex: pour suivre une ligne de grains): |
|---|
| 498 | c ------------------------------------------------------------- |
|---|
| 499 | |
|---|
| 500 | IF (iadv_tvl .EQ. 1) THEN |
|---|
| 501 | |
|---|
| 502 | call advect_tvl(timestep,temp,q,vu_f,vv_f,t_f,q_f |
|---|
| 503 | : ,d_t_adv,d_q_adv) |
|---|
| 504 | |
|---|
| 505 | do l = 1, nlayer |
|---|
| 506 | temp(l) = temp(l) + d_t_adv(l) |
|---|
| 507 | q(l,1) = q(l,1) + d_q_adv(l,1) |
|---|
| 508 | q(l,2) = q(l,2) + d_q_adv(l,2) |
|---|
| 509 | q(l,3) = q(l,3) + d_q_adv(l,3) |
|---|
| 510 | enddo |
|---|
| 511 | |
|---|
| 512 | ELSE |
|---|
| 513 | |
|---|
| 514 | do l = 1, nlayer |
|---|
| 515 | d_t_adv(l) = 0. |
|---|
| 516 | d_q_adv(l,1) = 0. |
|---|
| 517 | d_q_adv(l,2) = 0. |
|---|
| 518 | d_q_adv(l,3) = 0. |
|---|
| 519 | enddo |
|---|
| 520 | |
|---|
| 521 | ENDIF |
|---|
| 522 | |
|---|
| 523 | |
|---|
| 524 | c eventuellement, ajouter une tendance relative a la |
|---|
| 525 | c convergence de grande echelle: |
|---|
| 526 | c ------------------------------------------------------------- |
|---|
| 527 | |
|---|
| 528 | IF (i_cvg .EQ. 1) THEN |
|---|
| 529 | |
|---|
| 530 | call get_uvd(itap,timestep,tsol,qsol,file_fordat,ht,hq,hw) |
|---|
| 531 | |
|---|
| 532 | do l = 1, nlayer |
|---|
| 533 | |
|---|
| 534 | temp0(l) = temp(l) ! memoire de "l'avant dynamique" |
|---|
| 535 | q0(l,1) = q(l,1) |
|---|
| 536 | q0(l,2) = q(l,2) |
|---|
| 537 | q0(l,3) = q(l,3) |
|---|
| 538 | |
|---|
| 539 | d_t_cvg(l) = ht(l) * timestep |
|---|
| 540 | d_q_cvg(l) = hq(l) * timestep |
|---|
| 541 | temp(l) = temp(l) + d_t_cvg(l) |
|---|
| 542 | q(l,1) = q(l,1) + d_q_cvg(l) |
|---|
| 543 | c write(*,*)'d_t_cvg,d_q_cvg:',d_t_cvg(l),d_q_cvg(l) |
|---|
| 544 | if (q(l,1).lt.0.) then ! sb |
|---|
| 545 | print*,'OVAP negative dans main!' |
|---|
| 546 | print*,'itap,l,q,d_q_cvg:',itap,l,q(l,1),d_q_cvg(l),hq(l) |
|---|
| 547 | q(l,1) = MAX(q(l,1),1.e-10) ! evite les humidites negatives |
|---|
| 548 | endif |
|---|
| 549 | enddo |
|---|
| 550 | |
|---|
| 551 | if (itap.ge.30) then |
|---|
| 552 | print*,'hq(l),ht(l),dq,dt,q:' |
|---|
| 553 | : ,hq(16),ht(16),d_q_cvg(16),d_t_cvg(16),q(16,1) |
|---|
| 554 | endif |
|---|
| 555 | |
|---|
| 556 | ELSE |
|---|
| 557 | |
|---|
| 558 | do l = 1, nlayer |
|---|
| 559 | d_t_cvg(l) = 0. |
|---|
| 560 | d_q_cvg(l) = 0. |
|---|
| 561 | enddo |
|---|
| 562 | |
|---|
| 563 | ENDIF |
|---|
| 564 | |
|---|
| 565 | do l = 1, nlayer |
|---|
| 566 | dt_dyn(l) = d_t_cvg(l) / timestep |
|---|
| 567 | dq_dyn(l,1) = d_q_cvg(l) / timestep |
|---|
| 568 | dq_dyn(l,2) = 0.0 |
|---|
| 569 | dq_dyn(l,3) = 0.0 |
|---|
| 570 | enddo |
|---|
| 571 | |
|---|
| 572 | |
|---|
| 573 | c calcul des tendances physiques: |
|---|
| 574 | c ------------------------------- |
|---|
| 575 | |
|---|
| 576 | if (itsourcecont) then |
|---|
| 577 | do iq=5,nqmx |
|---|
| 578 | ll=min(iq-4,18) |
|---|
| 579 | if(firstcall) then |
|---|
| 580 | c q(ll,iq)=5/((plev(ll)-plev(ll+1))/rg) |
|---|
| 581 | c q(ll,iq)=5 |
|---|
| 582 | c q(ll,iq)=5*play(ll)/((plev(ll)-plev(ll+1))/rg)/(RD*temp(ll)) |
|---|
| 583 | q(ll,iq)=5*play(ll)/(RD*temp(ll)) |
|---|
| 584 | endif |
|---|
| 585 | print*,'q ll iq',q(ll,iq),ll,iq |
|---|
| 586 | enddo |
|---|
| 587 | else |
|---|
| 588 | do iq=5,nqmx |
|---|
| 589 | ll=min(iq-4,18) |
|---|
| 590 | q(ll,iq)=q(ll,iq)+1 |
|---|
| 591 | do l=1,llm |
|---|
| 592 | q(l,iq)=q(l,iq)*(1-timestep/tho) |
|---|
| 593 | print*,'q ll iq',q(ll,iq),ll,iq |
|---|
| 594 | enddo |
|---|
| 595 | enddo |
|---|
| 596 | |
|---|
| 597 | endif |
|---|
| 598 | |
|---|
| 599 | c CALL physiq(ecritphy,ngrid,nlayer,nqmx, |
|---|
| 600 | c : firstcall,lastcall, |
|---|
| 601 | c : day,day,time,timestep, |
|---|
| 602 | c : plev,play,phi,phis,paire, |
|---|
| 603 | c : u,v,temp,q, |
|---|
| 604 | c : du_dyn, dv_dyn, dt_dyn, dq_dyn, |
|---|
| 605 | c : w, |
|---|
| 606 | c : du,dv,dt,dq,dpsrf) |
|---|
| 607 | |
|---|
| 608 | c CALL physiq (ecritphy,ngrid,nlayer,nqmx,firstcall,lastcall, |
|---|
| 609 | c , day,day,time,timestep,plev,play,phi,phis,paire, |
|---|
| 610 | c , presnivs,clesphy0,u,v,temp,q, |
|---|
| 611 | c pour calculer proprement la tendance de cape liee a la |
|---|
| 612 | c dynamique, il faut entrer aussi temp0 et q0: |
|---|
| 613 | c , presnivs,clesphy0,u,v,temp,q,temp0,q0, |
|---|
| 614 | c , du_dyn,dv_dyn,dt_dyn,dq_dyn,w, |
|---|
| 615 | c - sorties |
|---|
| 616 | c s du,dv,dt,dq,dpsrf) |
|---|
| 617 | |
|---|
| 618 | print*,'PAS DE TEMPS ',timestep |
|---|
| 619 | print*,'ATTENTION!!! Il faudra passer temp0 et q0' |
|---|
| 620 | CALL physiq(ngrid,nlayer,nqmx, |
|---|
| 621 | : firstcall,lastcall, |
|---|
| 622 | : day,day,time,timestep, |
|---|
| 623 | : plev,play,phi,phis,paire,presnivs,clesphy0, |
|---|
| 624 | : u,v,temp,q, |
|---|
| 625 | : w, |
|---|
| 626 | : du,dv,dt,dq,dpsrf) |
|---|
| 627 | |
|---|
| 628 | |
|---|
| 629 | |
|---|
| 630 | |
|---|
| 631 | firstcall=.false. |
|---|
| 632 | |
|---|
| 633 | |
|---|
| 634 | c Ajout des tendances |
|---|
| 635 | c ------------------- |
|---|
| 636 | DO l=1,nlayer |
|---|
| 637 | u(l)=u(l)+timestep*du(l) |
|---|
| 638 | v(l)=v(l)+timestep*dv(l) |
|---|
| 639 | |
|---|
| 640 | temp(l)=temp(l)+timestep*dt(l) |
|---|
| 641 | ENDDO |
|---|
| 642 | |
|---|
| 643 | do iq=1,nqmx |
|---|
| 644 | do l=1,nlayer |
|---|
| 645 | q(l,iq)=q(l,iq)+timestep*dq(l,iq) |
|---|
| 646 | enddo |
|---|
| 647 | enddo |
|---|
| 648 | c write(78,rec=itap) ((q(l,iq),l=1,llm),iq=5,22) |
|---|
| 649 | write(78,rec=itap) ((q(l,iq),iq=5,22),l=1,llm) |
|---|
| 650 | itap = itap + 1 |
|---|
| 651 | |
|---|
| 652 | |
|---|
| 653 | time=time+timestep/rday |
|---|
| 654 | it=it+1 |
|---|
| 655 | if(mod(it,2000).eq.0) print*,'TIME=',time |
|---|
| 656 | if(time.gt.fnday) then |
|---|
| 657 | CALL endg1d(1,nlayer,play,int(time/ecritphy),timestep) |
|---|
| 658 | stop |
|---|
| 659 | endif |
|---|
| 660 | |
|---|
| 661 | GOTO 1 |
|---|
| 662 | |
|---|
| 663 | END |
|---|
| 664 | |
|---|
| 665 | c======================================================================= |
|---|
| 666 | c======================================================================= |
|---|
| 667 | c FIN DU PROGRAMMES |
|---|
| 668 | c CI-DESSOUS, QUELQUES SOUS-PROGRAMMES UTILS |
|---|
| 669 | c======================================================================= |
|---|
| 670 | c======================================================================= |
|---|
| 671 | |
|---|
| 672 | #include "1DUTILS.h" |
|---|
| 673 | #include "1Dconv.h" |
|---|
