source: trunk/mars/libf/phymars/meso_testphys1d.F @ 82

Last change on this file since 82 was 47, checked in by aslmd, 14 years ago

mars: ajout des fichiers mesoscale dans la physique du GCM martien

File size: 22.5 KB
Line 
1
2      PROGRAM meso_testphys1d
3! to use  'getin'
4      USE ioipsl_getincom
5      IMPLICIT NONE
6
7c=======================================================================
8c   subject:
9c   --------
10c   PROGRAM useful to run physical part of the martian GCM in a 1D column
11c       
12c Can be compiled with a command like (e.g. for 25 layers)
13c  "makegcm -p mars -d 25 testphys1d"
14c It requires the files "testphys1d.def" "callphys.def"
15c   and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
16c      and a file describing the sigma layers (e.g. "z2sig.def")
17c
18c   author: Frederic Hourdin, R.Fournier,F.Forget
19c   -------
20c   
21c   update: 12/06/2003 including chemistry (S. Lebonnois)
22c                            and water ice (F. Montmessin)
23c
24c=======================================================================
25
26#include "dimensions.h"
27#include "dimphys.h"
28#include "dimradmars.h"
29#include "comgeomfi.h"
30#include "surfdat.h"
31#include "comsoil.h"
32#include "comdiurn.h"
33#include "callkeys.h"
34#include "comcstfi.h"
35#include "planete.h"
36#include "comsaison.h"
37#include "yomaer.h"
38#include "control.h"
39#include "comvert.h"
40#include "netcdf.inc"
41#include "comg1d.h"
42#include "watercap.h"
43#include "logic.h"
44#include "advtrac.h"
45
46c --------------------------------------------------------------
47c  Declarations
48c --------------------------------------------------------------
49c
50      INTEGER unitstart      ! unite d'ecriture de "startfi"
51      INTEGER nlayer,nlevel,nsoil,ndt
52      INTEGER ilayer,ilevel,isoil,idt,iq
53      LOGICAl firstcall,lastcall
54c
55      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
56      REAL day           ! date durant le run
57      REAL time             ! time (0<time<1 ; time=0.5 a midi)
58      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
59      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
60      REAL psurf,tsurf     
61      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
62      REAL gru,grv   ! prescribed "geostrophic" background wind
63      REAL temp(nlayermx)   ! temperature at the middle of the layers
64      REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
65      REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
66      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
67      REAL co2ice           ! co2ice layer (kg.m-2)
68      REAL emis             ! surface layer
69      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
70      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
71
72c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
73      REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
74      REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
75      REAL dpsurf   
76      REAL dq(nlayermx,nqmx)
77      REAL dqdyn(nlayermx,nqmx)
78
79c   Various intermediate variables
80      INTEGER thermo
81      REAL zls
82      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
83      REAL pks, ptif, w(nlayermx)
84      REAL qtotinit, mqtot(nqmx),qtot
85      INTEGER ierr, aslun
86      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
87      Logical  tracerdyn
88      integer :: nq=1 ! number of tracers
89
90      character*2 str2
91      character (len=7) :: str7
92      character(len=44) :: txt
93
94c=======================================================================
95
96c=======================================================================
97c INITIALISATION
98c=======================================================================
99
100c ------------------------------------------------------
101c  Constantes prescrites ICI
102c ------------------------------------------------------
103
104      pi=2.E+0*asin(1.E+0)
105
106c     Constante de la Planete Mars
107c     ----------------------------
108      rad=3397200.               ! rayon de mars (m)  ~3397200 m
109      daysec=88775.              ! duree du sol (s)  ~88775 s
110      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
111      g=3.72                     ! gravite (m.s-2) ~3.72 
112      mugaz=43.49                ! Masse molaire de l'atm (g.mol-1) ~43.49
113      rcp=.256793                ! = r/cp  ~0.256793
114      r= 8.314511E+0 *1000.E+0/mugaz
115      cpp= r/rcp
116      year_day = 669           ! duree de l'annee (sols) ~668.6
117      periheli = 206.66          ! dist.min. soleil-mars (Mkm) ~206.66
118      aphelie = 249.22           ! dist.max. soleil-mars (Mkm) ~249.22
119      peri_day =  485.           ! date du perihelie (sols depuis printemps)
120      obliquit = 25.2           ! Obliquite de la planete (deg) ~25.2         
121 
122c     Parametres Couche limite et Turbulence
123c     --------------------------------------
124      z0 =  1.e-2                ! surface roughness (m) ~0.01
125      emin_turb = 1.e-6          ! energie minimale ~1.e-8
126      lmixmin = 30               ! longueur de melange ~100
127 
128c     propriete optiques des calottes et emissivite du sol
129c     ----------------------------------------------------
130      emissiv= 0.95              ! Emissivite du sol martien ~.95
131      emisice(1)=0.95            ! Emissivite calotte nord
132      emisice(2)=0.95            ! Emissivite calotte sud 
133      albedice(1)=0.5           ! Albedo calotte nord
134      albedice(2)=0.5            ! Albedo calotte sud
135      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
136      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
137      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
138      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
139      hybrid=.false.
140
141c ------------------------------------------------------
142c  Lecture des parametres dans "run.def"
143c ------------------------------------------------------
144
145
146! check if 'run.def' file is around (otherwise reading parameters
147! from callphys.def via getin() routine won't work.
148      open(99,file='run.def',status='old',form='formatted',
149     &     iostat=ierr)
150      if (ierr.ne.0) then
151        write(*,*) 'Cannot find required file "run.def"'
152        write(*,*) '  (which should contain some input parameters'
153        write(*,*) '   along with the following line:'
154        write(*,*) '   INCLUDEDEF=callphys.def'
155        write(*,*) '   )'
156        write(*,*) ' ... might as well stop here ...'
157        stop
158      else
159        close(99)
160      endif
161
162! check if we are going to run with or without tracers
163      write(*,*) "Run with or without tracer transport ?"
164      tracer=.false. ! default value
165      call getin("tracer",tracer)
166      write(*,*) " tracer = ",tracer
167
168! while we're at it, check if there is a 'traceur.def' file
169! and preocess it, if necessary. Otherwise initialize tracer names
170      if (tracer) then
171      ! load tracer names from file 'traceur.def'
172        open(90,file='traceur.def',status='old',form='formatted',
173     &       iostat=ierr)
174        if (ierr.ne.0) then
175          write(*,*) 'Cannot find required file "traceur.def"'
176          write(*,*) ' If you want to run with tracers, I need it'
177          write(*,*) ' ... might as well stop here ...'
178          stop
179        else
180          write(*,*) "testphys1d: Reading file traceur.def"
181          ! read number of tracers:
182          read(90,*,iostat=ierr) nq
183          if (ierr.ne.0) then
184            write(*,*) "testphys1d: error reading number of tracers"
185            write(*,*) "   (first line of traceur.def) "
186            stop
187          else
188            ! check that the number of tracers is indeed nqmx
189            if (nq.ne.nqmx) then
190              write(*,*) "testphys1d: error, wrong number of tracers:"
191              write(*,*) "nq=",nq," whereas nqmx=",nqmx
192              stop
193            endif
194          endif
195        endif
196        ! read tracer names from file traceur.def
197        do iq=1,nqmx
198          read(90,*,iostat=ierr) tnom(iq)
199          if (ierr.ne.0) then
200            write(*,*) 'testphys1d: error reading tracer names...'
201            stop
202          endif
203        enddo
204        close(90)
205
206        ! initialize tracers here:
207        write(*,*) "testphys1d: initializing tracers"
208        q(:,:)=0 ! default, set everything to zero
209        qsurf(:)=0
210        ! "smarter" initialization of some tracers
211        ! (get values from "profile_*" files, if these are available)
212        do iq=1,nqmx
213          txt=""
214          write(txt,"(a)") tnom(iq)
215          write(*,*)"  tracer:",trim(txt)
216          ! CO2
217          if (txt.eq."co2") then
218            q(:,iq)=0.95   ! kg /kg of atmosphere
219            qsurf(iq)=0. ! kg/m2 (not used for CO2)
220            ! even better, look for a "profile_co2" input file
221            open(91,file='profile_co2',status='old',
222     &       form='formatted',iostat=ierr)
223            if (ierr.eq.0) then
224              read(91,*) qsurf(iq)
225              do ilayer=1,nlayermx
226                read(91,*) q(ilayer,iq)
227              enddo
228            endif
229            close(91)
230          endif ! of if (txt.eq."co2")
231          ! WATER VAPOUR
232          if (txt.eq."h2o_vap") then
233            !look for a "profile_h2o_vap" input file
234            open(91,file='profile_h2o_vap',status='old',
235     &       form='formatted',iostat=ierr)
236            if (ierr.eq.0) then
237              read(91,*) qsurf(iq)
238              do ilayer=1,nlayermx
239                read(91,*) q(ilayer,iq)
240              enddo
241            else
242              write(*,*) "No profile_h2o_vap file!"
243            endif
244            close(91)
245          endif ! of if (txt.eq."h2o_ice")
246          ! WATER ICE
247          if (txt.eq."h2o_ice") then
248            !look for a "profile_h2o_vap" input file
249            open(91,file='profile_h2o_ice',status='old',
250     &       form='formatted',iostat=ierr)
251            if (ierr.eq.0) then
252              read(91,*) qsurf(iq)
253              do ilayer=1,nlayermx
254                read(91,*) q(ilayer,iq)
255              enddo
256            else
257              write(*,*) "No profile_h2o_ice file!"
258            endif
259            close(91)
260          endif ! of if (txt.eq."h2o_ice")
261          ! DUST
262          !if (txt(1:4).eq."dust") then
263          !  q(:,iq)=0.4    ! kg/kg of atmosphere
264          !  qsurf(iq)=100 ! kg/m2
265          !endif
266          ! DUST MMR
267          if (txt.eq."dust_mass") then
268            !look for a "profile_dust_mass" input file
269            open(91,file='profile_dust_mass',status='old',
270     &       form='formatted',iostat=ierr)
271            if (ierr.eq.0) then
272              read(91,*) qsurf(iq)
273              do ilayer=1,nlayermx
274                read(91,*) q(ilayer,iq)
275!                write(*,*) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
276              enddo
277            else
278              write(*,*) "No profile_dust_mass file!"
279            endif
280            close(91)
281          endif ! of if (txt.eq."dust_mass")
282          ! DUST NUMBER
283          if (txt.eq."dust_number") then
284            !look for a "profile_dust_number" input file
285            open(91,file='profile_dust_number',status='old',
286     &       form='formatted',iostat=ierr)
287            if (ierr.eq.0) then
288              read(91,*) qsurf(iq)
289              do ilayer=1,nlayermx
290                read(91,*) q(ilayer,iq)
291              enddo
292            else
293              write(*,*) "No profile_dust_number file!"
294            endif
295            close(91)
296          endif ! of if (txt.eq."dust_number")
297        enddo ! of do iq=1,nqmx
298        ! NB: some more initializations (chemistry) is done later
299
300      else
301      ! we still need to set (dummy) tracer names for physdem1
302        nq=nqmx
303        do iq=1,nq
304          write(str7,'(a1,i2.2)')'q',iq
305          tnom(iq)=str7
306        enddo
307      endif ! of if (tracer)
308
309c  Date et heure locale du debut du run
310c  ------------------------------------
311c    Date (en sols depuis le solstice de printemps) du debut du run
312      day0 = 0 ! default value for day0
313      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
314      call getin("day0",day0)
315      day=float(day0)
316      write(*,*) " day0 = ",day0
317c  Heure de demarrage
318      time=0 ! default value for time
319      write(*,*)'Initial local time (in hours, between 0 and 24)?'
320      call getin("time",time)
321      write(*,*)" time = ",time
322      time=time/24.E+0 ! convert time (hours) to fraction of sol
323
324c  Discretisation (Definition de la grille et des pas de temps)
325c  --------------
326c
327      nlayer=nlayermx
328      nlevel=nlayer+1
329      nsoil=nsoilmx
330
331      day_step=48 ! default value for day_step
332      PRINT *,'Number of time steps per sol ?'
333      call getin("day_step",day_step)
334      write(*,*) " day_step = ",day_step
335
336      ndt=10 ! default value for ndt
337      PRINT *,'Number of sols to run ?'
338      call getin("ndt",ndt)
339      write(*,*) " ndt = ",ndt
340
341      ndt=ndt*day_step     
342      dtphys=daysec/day_step 
343c Pression de surface sur la planete
344c ------------------------------------
345c
346      psurf=610. ! default value for psurf
347      PRINT *,'Surface pressure (Pa) ?'
348      call getin("psurf",psurf)
349      write(*,*) " psurf = ",psurf
350c Pression de reference
351      pa=20.
352      preff=610.     
353 
354c Proprietes des poussiere aerosol
355c --------------------------------
356      tauvis=0.2 ! default value for tauvis
357      print *,'Reference dust opacity at 700 Pa ?'
358      call getin("tauvis",tauvis)
359      write(*,*) " tauvis = ",tauvis
360 
361c  latitude/longitude
362c  ------------------
363      lati(1)=0 ! default value for lati(1)
364      PRINT *,'latitude (in degrees) ?'
365      call getin("latitude",lati(1))
366      write(*,*) " latitude = ",lati(1)
367      lati(1)=lati(1)*pi/180.E+0
368      long(1)=0.E+0
369      long(1)=long(1)*pi/180.E+0
370
371c  Initialisation albedo / inertie du sol
372c  --------------------------------------
373c
374      albedodat(1)=0.2 ! default value for albedodat
375      PRINT *,'Albedo of bare ground ?'
376      call getin("albedo",albedodat(1))
377      write(*,*) " albedo = ",albedodat(1)
378
379      inertiedat(1,1)=400 ! default value for inertiedat
380      PRINT *,'Soil thermal inertia (SI) ?'
381      call getin("inertia",inertiedat(1,1))
382      write(*,*) " inertia = ",inertiedat(1,1)
383c
384c  pour le schema d'ondes de gravite
385c  ---------------------------------
386c
387      zmea(1)=0.E+0
388      zstd(1)=0.E+0
389      zsig(1)=0.E+0
390      zgam(1)=0.E+0
391      zthe(1)=0.E+0
392
393
394
395c   Initialisation speciales "physiq"
396c   ---------------------------------
397c   la surface de chaque maille est inutile en 1D --->
398      area(1)=1.E+0
399
400c   le geopotentiel au sol est inutile en 1D car tout est controle
401c   par la pression de surface --->
402      phisfi(1)=0.E+0
403
404c  "inifis" reproduit un certain nombre d'initialisations deja faites
405c  + lecture des clefs de callphys.def
406c  + calcul de la frequence d'appel au rayonnement
407c  + calcul des sinus et cosinus des longitude latitude
408
409!Mars possible matter with dtphys in input and include!!!
410      CALL meso_inifis(1,llm,day0,daysec,dtphys,
411     .            lati,long,area,rad,g,r,cpp)
412c   Initialisation pour prendre en compte les vents en 1-D
413c   ------------------------------------------------------
414      ptif=2.E+0*omeg*sinlat(1)
415 
416c    vent geostrophique
417      gru=10. ! default value for gru
418      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
419      call getin("u",gru)
420      write(*,*) " u = ",gru
421      grv=0. !default value for grv
422      PRINT *,'meridional northward component of the geostrophic',
423     &' wind (m/s) ?'
424      call getin("v",grv)
425      write(*,*) " v = ",grv
426
427c     Initialisation des vents  au premier pas de temps
428      DO ilayer=1,nlayer
429         u(ilayer)=gru
430         v(ilayer)=grv
431      ENDDO
432
433c     energie cinetique turbulente
434      DO ilevel=1,nlevel
435         q2(ilevel)=0.E+0
436      ENDDO
437
438c  glace de CO2 au sol
439c  -------------------
440      co2ice=0.E+0 ! default value for co2ice
441      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
442      call getin("co2ice",co2ice)
443      write(*,*) " co2ice = ",co2ice
444
445c
446c  emissivite
447c  ----------
448      emis=emissiv
449      IF (co2ice.eq.1.E+0) THEN
450         emis=emisice(1) ! northern hemisphere
451         IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
452      ENDIF
453
454 
455
456c  calcul des pressions et altitudes en utilisant les niveaux sigma
457c  ----------------------------------------------------------------
458
459c    Vertical Coordinates
460c    """"""""""""""""""""
461      hybrid=.true.
462      PRINT *,'Hybrid coordinates ?'
463      call getin("hybrid",hybrid)
464      write(*,*) " hybrid = ", hybrid
465
466      CALL  disvert
467
468      DO ilevel=1,nlevel
469        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
470      ENDDO
471
472      DO ilayer=1,nlayer
473        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
474      ENDDO
475
476      DO ilayer=1,nlayer
477        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))
478     &   /g
479      ENDDO
480
481
482c  profil de temperature au premier appel
483c  --------------------------------------
484      pks=psurf**rcp
485
486c altitude en km dans profile: on divise zlay par 1000
487      tmp1(0)=0.E+0
488      DO ilayer=1,nlayer
489        tmp1(ilayer)=zlay(ilayer)/1000.E+0
490      ENDDO
491
492      call profile(nlayer+1,tmp1,tmp2)
493
494      tsurf=tmp2(0)
495      DO ilayer=1,nlayer
496        temp(ilayer)=tmp2(ilayer)
497      ENDDO
498     
499
500
501! Initialize soil properties and temperature
502! ------------------------------------------
503      volcapa=1.e6 ! volumetric heat capacity
504      DO isoil=1,nsoil
505         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
506         tsoil(isoil)=tsurf  ! soil temperature
507      ENDDO
508
509! Initialize depths
510! -----------------
511      do isoil=0,nsoil-1
512        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
513      enddo
514      do isoil=1,nsoil
515        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
516      enddo
517
518c    Initialisation des traceurs
519c    ---------------------------
520
521      if (photochem.or.callthermos) then
522         write(*,*) 'Especes chimiques initialisees'
523         ! thermo=0: initialisation pour toutes les couches
524         thermo=0
525         call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
526     $   thermo,qsurf)
527      endif
528      watercaptag(ngridmx)=.false.
529     
530
531c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
532c    ----------------------------------------------------------------
533c    (effectuee avec "writeg1d" a partir de "physiq.F"
534c    ou tout autre subroutine physique, y compris celle ci).
535
536        g1d_nlayer=nlayer
537        g1d_nomfich='g1d.dat'
538        g1d_unitfich=40
539        g1d_nomctl='g1d.ctl'
540        g1d_unitctl=41
541        g1d_premier=.true.
542        g2d_premier=.true.
543
544c  Ecriture de "startfi"
545c  --------------------
546c  (Ce fichier sera aussitot relu au premier
547c   appel de "physiq", mais il est necessaire pour passer
548c   les variables purement physiques a "physiq"...
549
550      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
551     .              dtphys,float(day0),
552     .              time,tsurf,tsoil,co2ice,emis,q2,qsurf,
553     .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
554c=======================================================================
555c  BOUCLE TEMPORELLE DU MODELE 1D
556c=======================================================================
557c
558      firstcall=.true.
559      lastcall=.false.
560
561      DO idt=1,ndt
562c        IF (idt.eq.ndt) lastcall=.true.
563        IF (idt.eq.ndt-day_step-1) then       !test
564         lastcall=.true.
565         call solarlong(day*1.0,zls)
566         write(103,*) 'Ls=',zls*180./pi
567         write(103,*) 'Lat=', lati(1)*180./pi
568         write(103,*) 'Tau=', tauvis/700*psurf
569         write(103,*) 'RunEnd - Atmos. Temp. File'
570         write(103,*) 'RunEnd - Atmos. Temp. File'
571         write(104,*) 'Ls=',zls*180./pi
572         write(104,*) 'Lat=', lati(1)
573         write(104,*) 'Tau=', tauvis/700*psurf
574         write(104,*) 'RunEnd - Atmos. Temp. File'
575        ENDIF
576
577c    calcul du geopotentiel
578c     ~~~~~~~~~~~~~~~~~~~~~
579      DO ilayer=1,nlayer
580        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
581        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
582      ENDDO
583      phi(1)=pks*h(1)*(1.E+0-s(1))
584      DO ilayer=2,nlayer
585         phi(ilayer)=phi(ilayer-1)+
586     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
587     &                  *(s(ilayer-1)-s(ilayer))
588
589      ENDDO
590
591c       appel de la physique
592c       --------------------
593
594      CALL meso_physiq (1,llm,nqmx,
595     ,     firstcall,lastcall,
596     ,     day,time,dtphys,
597     ,     plev,play,phi,
598     ,     u, v,temp, q, 
599     ,     w,
600C - sorties
601     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
602
603c       evolution du vent : modele 1D
604c       -----------------------------
605 
606c       la physique calcule les derivees temporelles de u et v.
607c       on y rajoute betement un effet Coriolis.
608c
609c       DO ilayer=1,nlayer
610c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
611c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
612c       ENDDO
613
614c       Pour certain test : pas de coriolis a l'equateur
615c       if(lati(1).eq.0.) then
616          DO ilayer=1,nlayer
617             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
618             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
619          ENDDO
620c       end if
621c     
622c
623c       Calcul du temps au pas de temps suivant
624c       ---------------------------------------
625        firstcall=.false.
626        time=time+dtphys/daysec
627        IF (time.gt.1.E+0) then
628            time=time-1.E+0
629            day=day+1
630        ENDIF
631
632c       calcul des vitesses et temperature au pas de temps suivant
633c       ----------------------------------------------------------
634
635        DO ilayer=1,nlayer
636           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
637           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
638           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
639        ENDDO
640
641c       calcul des pressions au pas de temps suivant
642c       ----------------------------------------------------------
643
644           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
645           DO ilevel=1,nlevel
646             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
647           ENDDO
648           DO ilayer=1,nlayer
649             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
650           ENDDO
651
652!       increment tracer values
653        DO iq = 1, nqmx
654          DO ilayer=1,nlayer
655             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
656          ENDDO
657        ENDDO
658
659      ENDDO   ! fin de la boucle temporelle
660
661c    ========================================================
662c    GESTION DES SORTIE
663c    ========================================================
664
665c    fermeture pour conclure les sorties format grads dans "g1d.dat"
666c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
667c    ou tout autre subroutine physique, y compris celle ci).
668
669c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
670        CALL endg1d(1,nlayer,zlay/1000.,ndt)
671
672c    ========================================================
673      END
674 
675c***********************************************************************
676c***********************************************************************
677c     Subroutines Bidons utilise seulement en 3D, mais
678c     necessaire a la compilation de testphys1d en 1-D
679
680      subroutine gr_fi_dyn
681      RETURN
682      END
683 
684c***********************************************************************
685c***********************************************************************
686
687#include "../dyn3d/disvert.F"
Note: See TracBrowser for help on using the repository browser.