source: trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F @ 537

Last change on this file since 537 was 534, checked in by jleconte, 13 years ago

Corrects a bug in rcm1d
  • Property svn:executable set to *
File size: 24.2 KB
RevLine 
[253]1      program rcm1d
2
3      use radcommon_h, only: tauvis
4! to use  'getin'
5      use ioipsl_getincom
6
7      implicit none
8
9!==================================================================
10!     
11!     Purpose
12!     -------
13!     Run the physics package of the universal model in a 1D column.
14!     
15!     It can be compiled with a command like (e.g. for 25 layers):
16!                                  "makegcm -p std -d 25 rcm1d"
17!     It requires the files "callphys.def", "z2sig.def",
18!     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
19!
20!     Authors
21!     -------
22!     Frederic Hourdin
23!     R. Fournier
24!     F. Forget
25!     F. Montmessin (water ice added)
26!     R. Wordsworth
27!
28!==================================================================
29
30#include "dimensions.h"
[534]31#include "paramet.h"
[253]32#include "dimphys.h"
33#include "comgeomfi.h"
34#include "surfdat.h"
35#include "comsoil.h"
36#include "comdiurn.h"
37#include "callkeys.h"
38#include "comcstfi.h"
39#include "planete.h"
40#include "comsaison.h"
41#include "control.h"
42#include "comvert.h"
43#include "netcdf.inc"
44#include "watercap.h"
45#include "fisice.h"
46#include "logic.h"
47#include "advtrac.h"
[534]48#include "comgeom.h"
[253]49
50c --------------------------------------------------------------
51c  Declarations
52c --------------------------------------------------------------
53c
54      INTEGER unitstart      ! unite d'ecriture de "startfi"
55      INTEGER nlayer,nlevel,nsoil,ndt
56      INTEGER ilayer,ilevel,isoil,idt,iq
57      LOGICAl firstcall,lastcall
58c
59      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
60      REAL day           ! date durant le run
61      REAL time             ! time (0<time<1 ; time=0.5 a midi)
62      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
63      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
64      REAL psurf,tsurf     
65      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
66      REAL gru,grv   ! prescribed "geostrophic" background wind
67      REAL temp(nlayermx)   ! temperature at the middle of the layers
68      REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
69      REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
70      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
71!      REAL co2ice           ! co2ice layer (kg.m-2) !not used anymore
72      integer :: i_co2_ice=0  ! tracer index of co2 ice
73      integer :: i_h2o_ice=0  ! tracer index of co2 ice
74      REAL emis             ! surface layer
75      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
76      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
77
78c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
79      REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
80      REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
81      REAL dpsurf   
82      REAL dq(nlayermx,nqmx)
83      REAL dqdyn(nlayermx,nqmx)
84
85c   Various intermediate variables
86!      INTEGER thermo
87      REAL zls
88      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
89      REAL pks, ptif, w(nlayermx)
90      REAL qtotinit, mqtot(nqmx),qtot
91      INTEGER ierr, aslun
92      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
93      Logical  tracerdyn
94      integer :: nq=1 ! number of tracers
95 
96      character*2 str2
97      character (len=7) :: str7
98
99      logical oldcompare, earthhack,saveprofile
100
101!     added by RW for zlay computation
102      real Hscale, Hmax, rho, dz
103
104!     added by RW for autozlevs computation
105      real nu, xx, pMIN, zlev, Htop
106      real logplevs(nlayermx)
107
108!     added by RDW to allow variations in cp with location
109      REAL cpp3D(ngridmx,nlayermx)   ! specific heat capacity at const. pressure
110      REAL rcp3D(ngridmx,nlayermx)   ! R / specific heat capacity
111     
112!     added by BC
113      REAL cloudfrac(ngridmx,nlayermx)
114      REAL hice(ngridmx),totcloudfrac(ngridmx)
115
116c=======================================================================
117c INITIALISATION
118c=======================================================================
119
[526]120      saveprofile=.false.
[253]121
122c ------------------------------------------------------
123c  Default values for constants (corresponding to Mars)
124c ------------------------------------------------------
125
126      pi=2.E+0*asin(1.E+0)
127
128c     Constante de la Planete Mars
129c     ----------------------------
130      rad=3397200.               ! rayon de mars (m)  ~3397200 m
131      daysec=88775.              ! duree du sol (s)  ~88775 s
132      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
133      g=3.72                     ! gravite (m.s-2) ~3.72 
134      mugaz=43.49                ! Masse molaire de l'atm (g.mol-1) ~43.49
135      rcp=.256793                ! = r/cp  ~0.256793
136      r= 8.314511E+0 *1000.E+0/mugaz
137      cpp= r/rcp
138      year_day = 669           ! duree de l'annee (sols) ~668.6
139      periastr = 206.66          ! min. dist. star-planet (Mkm) ~206.66
140      apoastr = 249.22           ! max. dist. star-planet (Mkm) ~249.22
141      peri_day =  485.           ! date du periastron (sols depuis printemps)
142      obliquit = 25.2           ! Obliquite de la planete (deg) ~25.2         
143 
144c     Parametres Couche limite et Turbulence
145c     --------------------------------------
146      z0 =  1.e-2                ! surface roughness (m) ~0.01
147      emin_turb = 1.e-6          ! energie minimale ~1.e-8
148      lmixmin = 30               ! longueur de melange ~100
149 
150c     propriete optiques des calottes et emissivite du sol
151c     ----------------------------------------------------
152      emissiv= 0.95              ! Emissivite du sol martien ~.95
153      emisice(1)=0.95            ! Emissivite calotte nord
154      emisice(2)=0.95            ! Emissivite calotte sud 
155
156      albedice(1)=0.5            ! Albedo calotte nord
157      albedice(2)=0.5            ! Albedo calotte sud
158      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
159      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
160      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
161      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
162      hybrid=.false.
163
164c ------------------------------------------------------
165c  Load parameters from "run.def"
166c ------------------------------------------------------
167
168! also check if 'run1d.def' file is around (otherwise reading parameters
169! from callphys.def via getin() routine won't work.
170      open(90,file='run.def',status='old',form='formatted',
171     &     iostat=ierr)
172      if (ierr.ne.0) then
173        write(*,*) 'Cannot find required file "run.def"'
174        write(*,*) '  (which should contain some input parameters'
175        write(*,*) '   along with the following line:'
176        write(*,*) '   INCLUDEDEF=callphys.def'
177        write(*,*) '   )'
178        write(*,*) ' ... might as well stop here ...'
179        stop
180      else
181        close(90)
182      endif
183
184! check if we are going to run with or without tracers
185      write(*,*) "Run with or without tracer transport ?"
186      tracer=.false. ! default value
187      call getin("tracer",tracer)
188      write(*,*) " tracer = ",tracer
189
190! while we're at it, check if there is a 'traceur.def' file
191! and preocess it, if necessary. Otherwise initialize tracer names
192      if (tracer) then
193      ! load tracer names from file 'traceur.def'
194        open(90,file='traceur.def',status='old',form='formatted',
195     &       iostat=ierr)
196        if (ierr.eq.0) then
197          write(*,*) "rcm1d: Reading file traceur.def"
198          ! read number of tracers:
199          read(90,*,iostat=ierr) nq
200          if (ierr.ne.0) then
201            write(*,*) "rcm1d: error reading number of tracers"
202            write(*,*) "   (first line of traceur.def) "
203            stop
204          else
205            ! check that the number of tracers is indeed nqmx
206            if (nq.ne.nqmx) then
207              write(*,*) "rcm1d: error, wrong number of tracers:"
208              write(*,*) "nq=",nq," whereas nqmx=",nqmx
209              stop
210            endif
211          endif
212       
213          ! initialize advection schemes to Van-Leer for all tracers
214          do iq=1,nq
215            iadv(iq)=3 ! Van-Leer
216          enddo
217       
218          do iq=1,nq
219          ! minimal version, just read in the tracer names, 1 per line
220            read(90,*,iostat=ierr) tnom(iq)
221            if (ierr.ne.0) then
222              write(*,*) 'rcm1d: error reading tracer names...'
223              stop
224            endif
225          enddo !of do iq=1,nq
226! check for co2_ice / h2o_ice tracers:
227         i_co2_ice=0
228         i_h2o_ice=0
229         do iq=1,nq
230           if (tnom(iq)=="co2_ice") then
231             i_co2_ice=iq
232           elseif (tnom(iq)=="h2o_ice") then
233             i_h2o_ice=iq
234           endif
235         enddo
236         !if (i_co2_ice==0) then
237         !  write(*,*) "rcm1d: error, we need a 'co2_ice' tracer"
238         !  write(*,*) "   (add one to traceur.def)"
239         !  stop
240         !endif
241        else
242          write(*,*) 'Cannot find required file "traceur.def"'
243          write(*,*) ' If you want to run with tracers, I need it'
244          write(*,*) ' ... might as well stop here ...'
245          stop
246        endif
247        close(90)
248
249      else
250      ! we still need to set (dummy) tracer names for physdem1
251        nq=nqmx
252        do iq=1,nq
253          write(str7,'(a1,i2.2)')'q',iq
254          tnom(iq)=str7
255        enddo
256        ! actually, we'll need at least one "co2_ice" tracer
257        ! (for surface CO2 ice)
258        tnom(1)="co2_ice"
259        i_co2_ice=1
260      endif ! of if (tracer)
261
262
263c  Date et heure locale du debut du run
264c  ------------------------------------
265c    Date (en sols depuis le solstice de printemps) du debut du run
266      day0 = 0 ! default value for day0
267      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
268      call getin("day0",day0)
269      day=float(day0)
270      write(*,*) " day0 = ",day0
271c  Heure de demarrage
272      time=0 ! default value for time
273      write(*,*)'Initial local time (in hours, between 0 and 24)?'
274      call getin("time",time)
275      write(*,*)" time = ",time
276      time=time/24.E+0 ! convert time (hours) to fraction of sol
277
[526]278
[253]279c  Discretisation (Definition de la grille et des pas de temps)
280c  --------------
281c
282      nlayer=nlayermx
283      nlevel=nlayer+1
284      nsoil=nsoilmx
285
286      day_step=48 ! default value for day_step
287      PRINT *,'Number of time steps per sol ?'
288      call getin("day_step",day_step)
289      write(*,*) " day_step = ",day_step
290
[534]291       
292      ecritphy=day_step ! default value for ecritphy
293      PRINT *,'Nunber of steps between writediagfi ?'
294      call getin("ecritphy",ecritphy)
295      write(*,*) " ecritphy = ",ecritphy
296
[253]297      ndt=10 ! default value for ndt
298      PRINT *,'Number of sols to run ?'
299      call getin("ndt",ndt)
300      write(*,*) " ndt = ",ndt
301
302      ndt=ndt*day_step     
303      dtphys=daysec/day_step 
304
305
[526]306c output spectrum?
307      write(*,*)"Output spectral OLR?"
308      specOLR=.false.
309      call getin("specOLR",specOLR)
310      write(*,*)" specOLR = ",specOLR
311
312
[253]313c Pression de surface sur la planete
314c ------------------------------------
315c
316      psurf=100000. ! default value for psurf
317      PRINT *,'Surface pressure (Pa) ?'
318      call getin("psurf",psurf)
319      write(*,*) " psurf = ",psurf
320c Pression de reference
[526]321      pa=psurf/30.
322      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
[253]323
324c Gravity of planet
325c -----------------
326      g=10.0 ! default value for g
327      PRINT *,'Gravity ?'
328      call getin("g",g)
329      write(*,*) " g = ",g
330
331      mugaz=0.
332      cpp=0.
333
334      call su_gases
335      call calc_cpp_mugaz
336     
337
338c Proprietes des poussiere aerosol
339c --------------------------------
340      tauvis=0.2 ! default value for tauvis
341      print *,'Reference dust opacity at 700 Pa ?'
342      call getin("tauvis",tauvis)
343      write(*,*) " tauvis = ",tauvis
344 
345c  latitude/longitude
346c  ------------------
347      lati(1)=0 ! default value for lati(1)
348      PRINT *,'latitude (in degrees) ?'
349      call getin("latitude",lati(1))
350      write(*,*) " latitude = ",lati(1)
351      lati(1)=lati(1)*pi/180.E+0
352      long(1)=0.E+0
353      long(1)=long(1)*pi/180.E+0
354
355c  periastron/apoastron
356c  --------------------
357      periastr=227.94
358      apoastr=227.94 ! default = mean martian values
359      PRINT *,'periastron (in AU) ?'
360      call getin("periastr",periastr)
361      write(*,*) "periastron = ",periastr
362      periastr=periastr*149.598 ! AU to Mkm
363      PRINT *,'apoastron (in AU) ?'
364      call getin("apoastr",apoastr)
365      write(*,*) "apoastron = ",apoastr
366      apoastr=apoastr*149.598 ! AU to Mkm
367
368c  obliquity
369c  ---------
370      obliquit=0.0! default=zero in 1d
371      PRINT *,'obliquity (in AU) ?'
372      call getin("obliquit",obliquit)
373      write(*,*) "obliquity = ",obliquit
374
375c  Initialisation albedo / inertie du sol
376c  --------------------------------------
377      albedodat(1)=0.2 ! default value for albedodat
378      PRINT *,'Albedo of bare ground ?'
379      call getin("albedo",albedodat(1))
380      write(*,*) " albedo = ",albedodat(1)
381
382      inertiedat(1,1)=400 ! default value for inertiedat
383      PRINT *,'Soil thermal inertia (SI) ?'
384      call getin("inertia",inertiedat(1,1))
385      write(*,*) " inertia = ",inertiedat(1,1)
386c
387c  pour le schema d'ondes de gravite
388c  ---------------------------------
389c
390      zmea(1)=0.E+0
391      zstd(1)=0.E+0
392      zsig(1)=0.E+0
393      zgam(1)=0.E+0
394      zthe(1)=0.E+0
395
396c    Initialisation des traceurs
397c    ---------------------------
398
399      DO iq=1,nqmx
400        DO ilayer=1,nlayer
401           q(ilayer,iq) = 0.
402        ENDDO
403      ENDDO
404     
405      DO iq=1,nqmx
406        qsurf(iq) = 0.
407      ENDDO
408
409c   Initialisation speciales "physiq"
410c   ---------------------------------
411c   la surface de chaque maille est inutile en 1D --->
412      area(1)=1.E+0
[534]413      aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
[253]414
415c   le geopotentiel au sol est inutile en 1D car tout est controle
416c   par la pression de surface --->
417      phisfi(1)=0.E+0
418
419c  "inifis" reproduit un certain nombre d'initialisations deja faites
420c  + lecture des clefs de callphys.def
421c  + calcul de la frequence d'appel au rayonnement
422c  + calcul des sinus et cosinus des longitude latitude
423
424!Mars possible issue with dtphys in input and include!!!
425      CALL inifis(1,llm,day0,daysec,dtphys,
426     .            lati,long,area,rad,g,r,cpp)
427c   Initialisation pour prendre en compte les vents en 1-D
428c   ------------------------------------------------------
429      ptif=2.E+0*omeg*sinlat(1)
430 
431
432c    vent geostrophique
433      gru=10. ! default value for gru
434      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
435      call getin("u",gru)
436      write(*,*) " u = ",gru
437      grv=0. !default value for grv
438      PRINT *,'meridional northward component of the geostrophic',
439     &' wind (m/s) ?'
440      call getin("v",grv)
441      write(*,*) " v = ",grv
442
443c     Initialisation des vents  au premier pas de temps
444      DO ilayer=1,nlayer
445         u(ilayer)=gru
446         v(ilayer)=grv
447      ENDDO
448
449c     energie cinetique turbulente
450      DO ilevel=1,nlevel
451         q2(ilevel)=0.E+0
452      ENDDO
453
454c  emissivity / surface co2 ice ( + h2o ice??)
455c  -------------------------------------------
456      emis=emissiv ! default value for emissivity
457      PRINT *,'Emissivity of bare ground ?'
458      call getin("emis",emis)
459      write(*,*) " emis = ",emis
460      emissiv=emis ! we do this so that condense_co2 sets things to the right
461                   ! value if there is no snow
462
463      if(i_co2_ice.gt.0)then
464         qsurf(i_co2_ice)=0 ! default value for co2ice
465         print*,'Initial CO2 ice on the surface (kg.m-2)'
466         call getin("co2ice",qsurf(i_co2_ice))
467         write(*,*) " co2ice = ",qsurf(i_co2_ice)
468         IF (qsurf(i_co2_ice).ge.1.E+0) THEN
469            ! if we have some CO2 ice on the surface, change emissivity
470            if (lati(1).ge.0) then ! northern hemisphere
471              emis=emisice(1)
472            else ! southern hemisphere
473              emis=emisice(2)
474            endif
475         ENDIF
476      endif
477
478c  calcul des pressions et altitudes en utilisant les niveaux sigma
479c  ----------------------------------------------------------------
480
481c    Vertical Coordinates
482c    """"""""""""""""""""
483      hybrid=.true.
484      PRINT *,'Hybrid coordinates ?'
485      call getin("hybrid",hybrid)
486      write(*,*) " hybrid = ", hybrid
487
488
489      autozlevs=.false.
490      PRINT *,'Auto-discretise vertical levels ?'
491      call getin("autozlevs",autozlevs)
492      write(*,*) " autozlevs = ", autozlevs
493
494      pceil=100.0 ! Pascals
495      PRINT *,'Ceiling pressure (Pa) ?'
496      call getin("pceil",pceil)
497      write(*,*) " pceil = ", pceil
498
499! Test of incompatibility:
500! if autozlevs used, cannot have hybrid too
501
502      if (autozlevs.and.hybrid) then
503         print*,'Cannot use autozlevs and hybrid together!'
504         call abort
505      endif
506
507      if(autozlevs)then
508           
[305]509         open(91,file="z2sig.def",form='formatted')
510         read(91,*) Hscale
511         DO ilayer=1,nlayer-2
512            read(91,*) Hmax
513         enddo
514         close(91)
[253]515 
516           
517         print*,'Hmax = ',Hmax,' km'
518         print*,'Auto-shifting Hscale to:'
519!     Hscale = Hmax / log(psurf/100.0)
520         Hscale = Hmax / log(psurf/pceil)
521         print*,'Hscale = ',Hscale,' km'
522         
523! none of this matters if we dont care about zlay
524         
525      endif
526
527      call disvert
528
529         if(.not.autozlevs)then
530            ! we want only the scale height from z2sig, in order to compute zlay
531            open(91,file="z2sig.def",form='formatted')
532            read(91,*) Hscale
533            close(91)
534         endif
535
536!         if(autozlevs)then
537!            open(94,file="Hscale.temp",form='formatted')
538!            read(94,*) Hscale
539!            close(94)
540!         endif
541
542         DO ilevel=1,nlevel
543            plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
544         ENDDO
545         
546         DO ilayer=1,nlayer
547            play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
548         ENDDO
549         
550
551
552         DO ilayer=1,nlayer
553!     zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
554!     &   /g
555            zlay(ilayer)=-1000.0*Hscale*log(play(ilayer)/plev(1))
556         ENDDO
557
558!      endif
559
560c  profil de temperature au premier appel
561c  --------------------------------------
562      pks=psurf**rcp
563
564c altitude en km dans profile: on divise zlay par 1000
565      tmp1(0)=0.E+0
566      DO ilayer=1,nlayer
567        tmp1(ilayer)=zlay(ilayer)/1000.E+0
568      ENDDO
569      call profile(nlayer+1,tmp1,tmp2)
570
571      tsurf=tmp2(0)
572      DO ilayer=1,nlayer
573        temp(ilayer)=tmp2(ilayer)
574      ENDDO
575     
576
577
578! Initialize soil properties and temperature
579! ------------------------------------------
580      volcapa=1.e6 ! volumetric heat capacity
581      DO isoil=1,nsoil
582         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
583         tsoil(isoil)=tsurf  ! soil temperature
584      ENDDO
585
586! Initialize depths
587! -----------------
588      do isoil=0,nsoil-1
589        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
590      enddo
591      do isoil=1,nsoil
592        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
593      enddo
594
595
596c  Ecriture de "startfi"
597c  --------------------
598c  (Ce fichier sera aussitot relu au premier
599c   appel de "physiq", mais il est necessaire pour passer
600c   les variables purement physiques a "physiq"...
601
602      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
603     &              dtphys,float(day0),
604     &              time,tsurf,tsoil,emis,q2,qsurf,
605     &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
606     &              cloudfrac,totcloudfrac,hice)
607
608c=======================================================================
609c  BOUCLE TEMPORELLE DU MODELE 1D
610c=======================================================================
611
612      firstcall=.true.
613      lastcall=.false.
614
615      DO idt=1,ndt
[486]616        IF (idt.eq.ndt) then       !test
[253]617         lastcall=.true.
618         call stellarlong(day*1.0,zls)
[526]619!         write(103,*) 'Ls=',zls*180./pi
620!         write(103,*) 'Lat=', lati(1)*180./pi
621!         write(103,*) 'Tau=', tauvis/700*psurf
622!         write(103,*) 'RunEnd - Atmos. Temp. File'
623!         write(103,*) 'RunEnd - Atmos. Temp. File'
624!         write(104,*) 'Ls=',zls*180./pi
625!         write(104,*) 'Lat=', lati(1)
626!         write(104,*) 'Tau=', tauvis/700*psurf
627!         write(104,*) 'RunEnd - Atmos. Temp. File'
[253]628        ENDIF
629
630c    calcul du geopotentiel
631c     ~~~~~~~~~~~~~~~~~~~~~
632
633        if(nonideal)then
634
635           DO ilayer=1,nlayer
636              call calc_cpp3d(cpp3D(1,ilayer),
637     $             rcp3D(1,ilayer),temp(ilayer),play(ilayer))
638           ENDDO
639
640           DO ilayer=1,nlayer
641
642!              if(autozlevs)then
643!                 s(ilayer)=(play(ilayer)/psurf)**rcp3D(1,ilayer)
644!              else
645                 s(ilayer)=
646     &                (aps(ilayer)/psurf+bps(ilayer))**rcp3D(1,ilayer)
647!              endif
648
649              h(ilayer)=cpp3D(1,ilayer)*temp(ilayer)/(pks*s(ilayer))
650           ENDDO
651
652        else
653           DO ilayer=1,nlayer
654
655!              if(autozlevs)then
656!                 s(ilayer)=(play(ilayer)/psurf)**rcp
657!              else
658                 s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
659!              endif
660              !s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
661              h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
662           ENDDO
663        endif
664!      DO ilayer=1,nlayer
665!        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
666!        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
667!      ENDDO
668      phi(1)=pks*h(1)*(1.E+0-s(1))
669      DO ilayer=2,nlayer
670         phi(ilayer)=phi(ilayer-1)+
671     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
672     &                  *(s(ilayer-1)-s(ilayer))
673
674      ENDDO
675
676c       appel de la physique
677c       --------------------
678
679
680      CALL physiq (1,llm,nqmx,
681     ,     firstcall,lastcall,
682     ,     day,time,dtphys,
683     ,     plev,play,phi,
684     ,     u, v,temp, q, 
685     ,     w,
686C - sorties
687     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
688
689
690c       evolution du vent : modele 1D
691c       -----------------------------
692 
693c       la physique calcule les derivees temporelles de u et v.
694c       on y rajoute betement un effet Coriolis.
695c
696c       DO ilayer=1,nlayer
697c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
698c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
699c       ENDDO
700
701c       Pour certain test : pas de coriolis a l'equateur
702c       if(lati(1).eq.0.) then
703          DO ilayer=1,nlayer
704             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
705             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
706          ENDDO
707c       end if
708c     
709c
710c       Calcul du temps au pas de temps suivant
711c       ---------------------------------------
712        firstcall=.false.
713        time=time+dtphys/daysec
714        IF (time.gt.1.E+0) then
715            time=time-1.E+0
716            day=day+1
717        ENDIF
718
719c       calcul des vitesses et temperature au pas de temps suivant
720c       ----------------------------------------------------------
721
722        DO ilayer=1,nlayer
723           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
724           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
725           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
726        ENDDO
727
728c       calcul des pressions au pas de temps suivant
729c       ----------------------------------------------------------
730
731           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
732           DO ilevel=1,nlevel
733              plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
734           ENDDO
735           DO ilayer=1,nlayer
736                 play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
737           ENDDO
738
739c       calcul traceur au pas de temps suivant
740c       --------------------------------------
741
742        DO iq = 1, nqmx
743          DO ilayer=1,nlayer
744             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
745          ENDDO
746        END DO
747
748      ENDDO   ! fin de la boucle temporelle
749
750c    ========================================================
751c    GESTION DES SORTIE
752c    ========================================================
753
754! save temperature profile
755      if(saveprofile)then
756         OPEN(12,file='profile.out',form='formatted')
757         DO ilayer=1,nlayermx
758            write(12,*) temp(ilayer), play(ilayer)
759         ENDDO
760         CLOSE(12)
761      endif
762
763
764c    ========================================================
765      end                       !rcm1d
766 
767c***********************************************************************
768c***********************************************************************
769c     Subroutines Bidons utilise seulement en 3D, mais
770c     necessaire a la compilation de rcm1d en 1D
771
772      subroutine gr_fi_dyn
773      RETURN
774      END
775 
776c***********************************************************************
777c***********************************************************************
778
779#include "../dyn3d/disvert.F"
Note: See TracBrowser for help on using the repository browser.