source: trunk/LMDZ.MARS/libf/phymars/testphys1d.F @ 259

Last change on this file since 259 was 246, checked in by aslmd, 14 years ago

MESOSCALE: ifort is picky when using float() on real numbers (corrected in meso_inc_inifisini.F). also changes in r234 not impacted in testphys1d.F for newphys, now corrected.
File size: 23.3 KB
RevLine 
[38]1
2      PROGRAM 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     --------------------------------------
[224]124      z0_default =  1.e-2                ! surface roughness (m) ~0.01
[38]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")
[161]231          ! Allow for an initial profile of argon
232          ! Can also be used to introduce a decaying tracer
233          ! in the 1D (TBD) to study thermals
234          if (txt.eq."ar") then
235            !look for a "profile_ar" input file
236            open(91,file='profile_ar',status='old',
237     &       form='formatted',iostat=ierr)
238            if (ierr.eq.0) then
239              read(91,*) qsurf(iq)
240              do ilayer=1,nlayermx
241                read(91,*) q(ilayer,iq)
242              enddo
243            else
244              write(*,*) "No profile_ar file!"
245            endif
246            close(91)
247          endif ! of if (txt.eq."ar")
248
[38]249          ! WATER VAPOUR
250          if (txt.eq."h2o_vap") then
251            !look for a "profile_h2o_vap" input file
252            open(91,file='profile_h2o_vap',status='old',
253     &       form='formatted',iostat=ierr)
254            if (ierr.eq.0) then
255              read(91,*) qsurf(iq)
256              do ilayer=1,nlayermx
257                read(91,*) q(ilayer,iq)
258              enddo
259            else
260              write(*,*) "No profile_h2o_vap file!"
261            endif
262            close(91)
263          endif ! of if (txt.eq."h2o_ice")
264          ! WATER ICE
265          if (txt.eq."h2o_ice") then
266            !look for a "profile_h2o_vap" input file
267            open(91,file='profile_h2o_ice',status='old',
268     &       form='formatted',iostat=ierr)
269            if (ierr.eq.0) then
270              read(91,*) qsurf(iq)
271              do ilayer=1,nlayermx
272                read(91,*) q(ilayer,iq)
273              enddo
274            else
275              write(*,*) "No profile_h2o_ice file!"
276            endif
277            close(91)
278          endif ! of if (txt.eq."h2o_ice")
279          ! DUST
280          !if (txt(1:4).eq."dust") then
281          !  q(:,iq)=0.4    ! kg/kg of atmosphere
282          !  qsurf(iq)=100 ! kg/m2
283          !endif
284          ! DUST MMR
285          if (txt.eq."dust_mass") then
286            !look for a "profile_dust_mass" input file
287            open(91,file='profile_dust_mass',status='old',
288     &       form='formatted',iostat=ierr)
289            if (ierr.eq.0) then
290              read(91,*) qsurf(iq)
291              do ilayer=1,nlayermx
292                read(91,*) q(ilayer,iq)
293!                write(*,*) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
294              enddo
295            else
296              write(*,*) "No profile_dust_mass file!"
297            endif
298            close(91)
299          endif ! of if (txt.eq."dust_mass")
300          ! DUST NUMBER
301          if (txt.eq."dust_number") then
302            !look for a "profile_dust_number" input file
303            open(91,file='profile_dust_number',status='old',
304     &       form='formatted',iostat=ierr)
305            if (ierr.eq.0) then
306              read(91,*) qsurf(iq)
307              do ilayer=1,nlayermx
308                read(91,*) q(ilayer,iq)
309              enddo
310            else
311              write(*,*) "No profile_dust_number file!"
312            endif
313            close(91)
314          endif ! of if (txt.eq."dust_number")
315        enddo ! of do iq=1,nqmx
316        ! NB: some more initializations (chemistry) is done later
317
318      else
319      ! we still need to set (dummy) tracer names for physdem1
320        nq=nqmx
321        do iq=1,nq
322          write(str7,'(a1,i2.2)')'q',iq
323          tnom(iq)=str7
324        enddo
325      endif ! of if (tracer)
326
327c  Date et heure locale du debut du run
328c  ------------------------------------
329c    Date (en sols depuis le solstice de printemps) du debut du run
330      day0 = 0 ! default value for day0
331      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
332      call getin("day0",day0)
333      day=float(day0)
334      write(*,*) " day0 = ",day0
335c  Heure de demarrage
336      time=0 ! default value for time
337      write(*,*)'Initial local time (in hours, between 0 and 24)?'
338      call getin("time",time)
339      write(*,*)" time = ",time
340      time=time/24.E+0 ! convert time (hours) to fraction of sol
341
342c  Discretisation (Definition de la grille et des pas de temps)
343c  --------------
344c
345      nlayer=nlayermx
346      nlevel=nlayer+1
347      nsoil=nsoilmx
348
349      day_step=48 ! default value for day_step
350      PRINT *,'Number of time steps per sol ?'
351      call getin("day_step",day_step)
352      write(*,*) " day_step = ",day_step
353
354      ndt=10 ! default value for ndt
355      PRINT *,'Number of sols to run ?'
356      call getin("ndt",ndt)
357      write(*,*) " ndt = ",ndt
358
359      ndt=ndt*day_step     
360      dtphys=daysec/day_step 
361c Pression de surface sur la planete
362c ------------------------------------
363c
364      psurf=610. ! default value for psurf
365      PRINT *,'Surface pressure (Pa) ?'
366      call getin("psurf",psurf)
367      write(*,*) " psurf = ",psurf
368c Pression de reference
369      pa=20.
370      preff=610.     
371 
372c Proprietes des poussiere aerosol
373c --------------------------------
374      tauvis=0.2 ! default value for tauvis
375      print *,'Reference dust opacity at 700 Pa ?'
376      call getin("tauvis",tauvis)
377      write(*,*) " tauvis = ",tauvis
378 
379c  latitude/longitude
380c  ------------------
381      lati(1)=0 ! default value for lati(1)
382      PRINT *,'latitude (in degrees) ?'
383      call getin("latitude",lati(1))
384      write(*,*) " latitude = ",lati(1)
385      lati(1)=lati(1)*pi/180.E+0
386      long(1)=0.E+0
387      long(1)=long(1)*pi/180.E+0
388
389c  Initialisation albedo / inertie du sol
390c  --------------------------------------
391c
392      albedodat(1)=0.2 ! default value for albedodat
393      PRINT *,'Albedo of bare ground ?'
394      call getin("albedo",albedodat(1))
395      write(*,*) " albedo = ",albedodat(1)
396
397      inertiedat(1,1)=400 ! default value for inertiedat
398      PRINT *,'Soil thermal inertia (SI) ?'
399      call getin("inertia",inertiedat(1,1))
400      write(*,*) " inertia = ",inertiedat(1,1)
[224]401
402      z0(1)=z0_default ! default value for roughness
403      write(*,*) 'Surface roughness length z0 (m)?'
404      call getin("z0",z0(1))
405      write(*,*) " z0 = ",z0(1)
[38]406c
407c  pour le schema d'ondes de gravite
408c  ---------------------------------
409c
410      zmea(1)=0.E+0
411      zstd(1)=0.E+0
412      zsig(1)=0.E+0
413      zgam(1)=0.E+0
414      zthe(1)=0.E+0
415
416
417
418c   Initialisation speciales "physiq"
419c   ---------------------------------
420c   la surface de chaque maille est inutile en 1D --->
421      area(1)=1.E+0
422
423c   le geopotentiel au sol est inutile en 1D car tout est controle
424c   par la pression de surface --->
425      phisfi(1)=0.E+0
426
427c  "inifis" reproduit un certain nombre d'initialisations deja faites
428c  + lecture des clefs de callphys.def
429c  + calcul de la frequence d'appel au rayonnement
430c  + calcul des sinus et cosinus des longitude latitude
431
432!Mars possible matter with dtphys in input and include!!!
433      CALL inifis(1,llm,day0,daysec,dtphys,
434     .            lati,long,area,rad,g,r,cpp)
435c   Initialisation pour prendre en compte les vents en 1-D
436c   ------------------------------------------------------
437      ptif=2.E+0*omeg*sinlat(1)
438 
439c    vent geostrophique
440      gru=10. ! default value for gru
441      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
442      call getin("u",gru)
443      write(*,*) " u = ",gru
444      grv=0. !default value for grv
445      PRINT *,'meridional northward component of the geostrophic',
446     &' wind (m/s) ?'
447      call getin("v",grv)
448      write(*,*) " v = ",grv
449
450c     Initialisation des vents  au premier pas de temps
451      DO ilayer=1,nlayer
452         u(ilayer)=gru
453         v(ilayer)=grv
454      ENDDO
455
456c     energie cinetique turbulente
457      DO ilevel=1,nlevel
458         q2(ilevel)=0.E+0
459      ENDDO
460
461c  glace de CO2 au sol
462c  -------------------
463      co2ice=0.E+0 ! default value for co2ice
464      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
465      call getin("co2ice",co2ice)
466      write(*,*) " co2ice = ",co2ice
467
468c
469c  emissivite
470c  ----------
471      emis=emissiv
472      IF (co2ice.eq.1.E+0) THEN
473         emis=emisice(1) ! northern hemisphere
474         IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
475      ENDIF
476
477 
478
479c  calcul des pressions et altitudes en utilisant les niveaux sigma
480c  ----------------------------------------------------------------
481
482c    Vertical Coordinates
483c    """"""""""""""""""""
484      hybrid=.true.
485      PRINT *,'Hybrid coordinates ?'
486      call getin("hybrid",hybrid)
487      write(*,*) " hybrid = ", hybrid
488
489      CALL  disvert
490
491      DO ilevel=1,nlevel
492        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
493      ENDDO
494
495      DO ilayer=1,nlayer
496        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
497      ENDDO
498
499      DO ilayer=1,nlayer
500        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))
501     &   /g
502      ENDDO
503
504
505c  profil de temperature au premier appel
506c  --------------------------------------
507      pks=psurf**rcp
508
509c altitude en km dans profile: on divise zlay par 1000
510      tmp1(0)=0.E+0
511      DO ilayer=1,nlayer
512        tmp1(ilayer)=zlay(ilayer)/1000.E+0
513      ENDDO
514
515      call profile(nlayer+1,tmp1,tmp2)
516
517      tsurf=tmp2(0)
518      DO ilayer=1,nlayer
519        temp(ilayer)=tmp2(ilayer)
520      ENDDO
521     
522
523
524! Initialize soil properties and temperature
525! ------------------------------------------
526      volcapa=1.e6 ! volumetric heat capacity
527      DO isoil=1,nsoil
528         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
529         tsoil(isoil)=tsurf  ! soil temperature
530      ENDDO
531
532! Initialize depths
533! -----------------
534      do isoil=0,nsoil-1
535        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
536      enddo
537      do isoil=1,nsoil
538        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
539      enddo
540
541c    Initialisation des traceurs
542c    ---------------------------
543
544      if (photochem.or.callthermos) then
545         write(*,*) 'Especes chimiques initialisees'
546         ! thermo=0: initialisation pour toutes les couches
547         thermo=0
548         call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
549     $   thermo,qsurf)
550      endif
551      watercaptag(ngridmx)=.false.
552     
553
554c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
555c    ----------------------------------------------------------------
556c    (effectuee avec "writeg1d" a partir de "physiq.F"
557c    ou tout autre subroutine physique, y compris celle ci).
558
559        g1d_nlayer=nlayer
560        g1d_nomfich='g1d.dat'
561        g1d_unitfich=40
562        g1d_nomctl='g1d.ctl'
563        g1d_unitctl=41
564        g1d_premier=.true.
565        g2d_premier=.true.
566
567c  Ecriture de "startfi"
568c  --------------------
569c  (Ce fichier sera aussitot relu au premier
570c   appel de "physiq", mais il est necessaire pour passer
571c   les variables purement physiques a "physiq"...
572
573      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
574     .              dtphys,float(day0),
575     .              time,tsurf,tsoil,co2ice,emis,q2,qsurf,
576     .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
577c=======================================================================
578c  BOUCLE TEMPORELLE DU MODELE 1D
579c=======================================================================
580c
581      firstcall=.true.
582      lastcall=.false.
583
584      DO idt=1,ndt
585c        IF (idt.eq.ndt) lastcall=.true.
586        IF (idt.eq.ndt-day_step-1) then       !test
587         lastcall=.true.
588         call solarlong(day*1.0,zls)
589         write(103,*) 'Ls=',zls*180./pi
590         write(103,*) 'Lat=', lati(1)*180./pi
591         write(103,*) 'Tau=', tauvis/700*psurf
592         write(103,*) 'RunEnd - Atmos. Temp. File'
593         write(103,*) 'RunEnd - Atmos. Temp. File'
594         write(104,*) 'Ls=',zls*180./pi
595         write(104,*) 'Lat=', lati(1)
596         write(104,*) 'Tau=', tauvis/700*psurf
597         write(104,*) 'RunEnd - Atmos. Temp. File'
598        ENDIF
599
600c    calcul du geopotentiel
601c     ~~~~~~~~~~~~~~~~~~~~~
602      DO ilayer=1,nlayer
603        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
604        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
605      ENDDO
606      phi(1)=pks*h(1)*(1.E+0-s(1))
607      DO ilayer=2,nlayer
608         phi(ilayer)=phi(ilayer-1)+
609     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
610     &                  *(s(ilayer-1)-s(ilayer))
611
612      ENDDO
613
614c       appel de la physique
615c       --------------------
616      CALL physiq (1,llm,nqmx,
617     ,     firstcall,lastcall,
618     ,     day,time,dtphys,
619     ,     plev,play,phi,
620     ,     u, v,temp, q, 
621     ,     w,
622C - sorties
623     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
624
625c       evolution du vent : modele 1D
626c       -----------------------------
627 
628c       la physique calcule les derivees temporelles de u et v.
629c       on y rajoute betement un effet Coriolis.
630c
631c       DO ilayer=1,nlayer
632c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
633c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
634c       ENDDO
635
636c       Pour certain test : pas de coriolis a l'equateur
637c       if(lati(1).eq.0.) then
638          DO ilayer=1,nlayer
639             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
640             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
641          ENDDO
642c       end if
643c     
644c
645c       Calcul du temps au pas de temps suivant
646c       ---------------------------------------
647        firstcall=.false.
648        time=time+dtphys/daysec
649        IF (time.gt.1.E+0) then
650            time=time-1.E+0
651            day=day+1
652        ENDIF
653
654c       calcul des vitesses et temperature au pas de temps suivant
655c       ----------------------------------------------------------
656
657        DO ilayer=1,nlayer
658           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
659           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
660           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
661        ENDDO
662
663c       calcul des pressions au pas de temps suivant
664c       ----------------------------------------------------------
665
666           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
667           DO ilevel=1,nlevel
668             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
669           ENDDO
670           DO ilayer=1,nlayer
671             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
672           ENDDO
673
674!       increment tracer values
675        DO iq = 1, nqmx
676          DO ilayer=1,nlayer
677             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
678          ENDDO
679        ENDDO
680
681      ENDDO   ! fin de la boucle temporelle
682
683c    ========================================================
684c    GESTION DES SORTIE
685c    ========================================================
686
687c    fermeture pour conclure les sorties format grads dans "g1d.dat"
688c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
689c    ou tout autre subroutine physique, y compris celle ci).
690
691c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
692        CALL endg1d(1,nlayer,zlay/1000.,ndt)
693
694c    ========================================================
695      END
696 
697c***********************************************************************
698c***********************************************************************
699c     Subroutines Bidons utilise seulement en 3D, mais
700c     necessaire a la compilation de testphys1d en 1-D
701
702      subroutine gr_fi_dyn
703      RETURN
704      END
705 
706c***********************************************************************
707c***********************************************************************
708
709#include "../dyn3d/disvert.F"
Note: See TracBrowser for help on using the repository browser.