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

Last change on this file since 533 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

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