source: trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F @ 899

Last change on this file since 899 was 894, checked in by slebonnois, 12 years ago

SL: small modifications on rcm1d (Venus, Titan)

File size: 13.8 KB
Line 
1
2      PROGRAM rcm1d
3     
4      USE infotrac
5      use control_mod
6      use comgeomphy
7      IMPLICIT NONE
8
9c=======================================================================
10c   subject:
11c   --------
12c   PROGRAM useful to run physical part of the venusian GCM in a 1D column
13c       
14c Can be compiled with a command like (e.g. for 50 layers)
15c  "makelmdz -p venus -d 50 rcm1d"
16
17c It requires the files "rcm1d.def" "physiq.def"
18c      and a file describing the sigma layers (e.g. "z2sig.def")
19c
20c   author: Frederic Hourdin, R.Fournier,F.Forget (original Mars version)
21c   ------- Sebastien Lebonnois (Venus version)
22c   
23c=======================================================================
24
25#include "dimensions.h"
26#include "dimsoil.h"
27#include "comcstfi.h"
28#include "comvert.h"
29#include "netcdf.inc"
30#include "logic.h"
31#include "clesphys.h"
32#include "iniprint.h"
33#include "tabcontrol.h"
34
35c --------------------------------------------------------------
36c  Declarations
37c --------------------------------------------------------------
38c
39      INTEGER unit           ! unite de lecture de "rcm1d.def"
40      INTEGER unitstart      ! unite d'ecriture de "startphy.nc"
41      INTEGER nlayer,nlevel,nsoil,ndt
42      INTEGER ilayer,ilevel,isoil,idt,iq
43      LOGICAl firstcall,lastcall
44c
45      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
46      REAL day              ! date durant le run
47      REAL time             ! time (0<time<1 ; time=0.5 a midi)
48      REAL play(llm)   ! Pressure at the middle of the layers (Pa)
49      REAL plev(llm+1) ! intermediate pressure levels (pa)
50      REAL psurf,tsurf(1)     
51      REAL u(llm),v(llm)  ! zonal, meridional wind
52      REAL gru,grv   ! prescribed "geostrophic" background wind
53      REAL temp(llm)   ! temperature at the middle of the layers
54      REAL,allocatable :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
55      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
56      REAL zlay(llm)   ! altitude estimee dans les couches (km)
57      REAL long(1),lati(1),area(1)
58      REAL cufi(1),cvfi(1)
59      REAL phisfi(1),albedo(1)
60      REAL solsw(1),sollwdown(1),dlw(1),radsol(1)
61      REAL zmea(1), zstd(1)
62      REAL zsig(1), zgam(1), zthe(1)
63      REAL zpic(1), zval(1)
64
65c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
66      REAL du(llm),dv(llm),dtemp(llm)
67      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
68      REAL dpsurf   
69      REAL,allocatable :: dq(:,:)
70
71c   Various intermediate variables
72      REAL zls
73      REAL phi(llm),s(llm)
74      REAL pk(llm),pks, w(llm)
75      INTEGER l, ierr, aslun
76      REAL tmp1(0:llm),tmp2(0:llm)                       
77
78      character*2 str2
79
80c normalement dans dyn3d/comconst.h
81      COMMON/cpdetvenus/cppdyn,nu_venus,t0_venus
82      REAL cppdyn,nu_venus,t0_venus
83      real pi
84
85c=======================================================================
86
87c=======================================================================
88c INITIALISATION
89c=======================================================================
90
91      lunout = 6
92
93c ------------------------------------------------------
94c  Constantes prescrites ICI
95c ------------------------------------------------------
96
97      pi=2.E+0*asin(1.E+0)
98
99c     Constante de la Planete Venus
100c     -----------------------------
101      planet_type = "venus"
102      rad=6051300.               ! rayon de Venus (m)  ~6051300 m
103      daysec= 1.0087e7           ! duree du sol (s)  ~1.e7 s
104      omeg=4.*asin(1.)/19.4141e6 ! vitesse de rotation (rad.s-1)
105      g= 8.87                    ! gravite (m.s-2) ~8.87
106      mugaz=43.44                ! Masse molaire de l'atm (g.mol-1) ~43.44
107! ADAPTATION GCM POUR CP(T)
108! VENUS: Cp(T) = cpp*(T/T0)^nu
109! avec T0=460. et nu=0.35
110      cpp=1.0e3
111!     cpp=9.0e2      ! version constante
112      cppdyn=cpp
113      r= 8.314511E+0 *1000.E+0/mugaz
114      rcp= r/cpp
115
116c-----------------------------------------------------------------------
117c   Initialisation des traceurs
118c   ---------------------------
119c  Choix du nombre de traceurs et du schema pour l'advection
120c  dans fichier traceur.def
121      call infotrac_init
122
123c Allocation de la tableau q : champs advectes   
124      allocate(q(llm,nqtot))
125      allocate(dq(llm,nqtot))
126
127c ------------------------------------------------------
128c  Lecture des parametres dans "rcm1d.def"
129c ------------------------------------------------------
130
131c   Opening parameters file "rcm1d.def"
132c   ---------------------------------------
133      unit =97
134      OPEN(unit,file='rcm1d.def',status='old',form='formatted'
135     .     ,iostat=ierr)
136
137      IF(ierr.ne.0) THEN
138        write(*,*) 'Problem to open "rcm1d.def'
139        write(*,*) 'Is it there ?'
140        stop
141      END IF
142
143c  Date et heure locale du debut du run
144c  ------------------------------------
145c    Date (en sols depuis le solstice de printemps) du debut du run
146      day0 = 0
147      PRINT *,'date de depart ?'
148      READ(unit,*) day0
149      day=REAL(day0)
150      PRINT *,day0
151c  Heure de demarrage
152      PRINT *,'heure de debut de simulation (entre 0 et 24) ?'
153      READ(unit,*) time
154      time=time/24.E+0
155
156c  Discretisation (Definition de la grille et des pas de temps)
157c  --------------
158c
159      nlayer=llm
160      nlevel=nlayer+1
161      nsoil=nsoilmx
162      PRINT *,'nombre de pas de temps par jour ?'
163      READ(unit,*) day_step
164      print*,day_step
165
166c     PRINT *,'nombre d appel au rayonnement par jour ?'
167c     READ(unit,*) nbapp_rad
168c     print*,nbapp_rad
169c LU DANS PHYSIQ.DEF...
170      nbapp_rad = 1000.
171
172      PRINT *,'nombre de jours simules ?'
173      READ(unit,*) ndt
174      print*,ndt
175
176      ndt=ndt*day_step     
177      dtphys=daysec/day_step 
178
179c Pression de surface sur la planete
180c ------------------------------------
181c
182      PRINT *,'pression au sol'
183      READ(unit,*) psurf
184      PRINT *,psurf
185c Pression de reference  ! voir dyn3d/etat0_venus
186c     pa     =  5.e4
187      pa     =  1.e6
188      preff  = 9.2e6 ! 92 bars
189c     preff  = psurf
190 
191c  latitude/longitude
192c  -------------------
193      PRINT *,'latitude en degres ?'
194      READ(unit,*) lati(1)
195      PRINT *,lati(1)
196      long(1)=0.E+0
197
198c  Initialisation albedo
199c  ----------------------
200c ne sert pas ici...
201      albedo(1)=0.1
202c      PRINT *,'Albedo du sol nu ?'
203c      READ(unit,*) albedo(1)
204c      PRINT *,albedo(1)
205
206c   Initialisation speciales "physiq"
207c   ---------------------------------
208
209      CALL init_phys_lmdz(iim,jjm,llm,1,(/1/))
210      call initcomgeomphy
211      call infotrac_init
212      call ini_cpdet
213
214c   la surface de chaque maille est inutile en 1D --->
215      area(1)=1.E+0
216c de meme ?
217      cufi(1)=1.E+0
218      cvfi(1)=1.E+0
219
220c   le geopotentiel au sol est inutile en 1D car tout est controle
221c   par la pression de surface --->
222      phisfi(1)=0.E+0
223
224      CALL iniphysiq(1,llm,daysec,day0,dtphys,
225     .            lati,long,area,cufi,cvfi,rad,g,r,cpp)
226
227c   Initialisation pour prendre en compte les vents en 1-D
228c   ------------------------------------------------------
229 
230c    vent geostrophique
231      PRINT *,'composante vers l est du vent geostrophique (U) ?'
232      READ(unit,*) gru
233      PRINT *,'composante vers le nord du vent geostrophique (V) ?'
234      READ(unit,*) grv
235
236c     Initialisation des vents  au premier pas de temps
237      DO ilayer=1,nlayer
238         u(ilayer)=gru
239         v(ilayer)=grv
240      ENDDO
241
242c  calcul des pressions et altitudes en utilisant les niveaux sigma
243c  ----------------------------------------------------------------
244
245c    Vertical Coordinates  (hybrids)
246c    """"""""""""""""""""
247      CALL  disvert_noterre
248     
249c     Calcul au milieu des couches : Vient de la version Mars
250c     WARNING : le choix de placer le milieu des couches au niveau de
251c     pression intermédiaire est arbitraire et pourrait etre modifié.
252c     C'est fait de la meme facon dans disvert
253
254      DO l = 1, llm
255       aps(l) =  0.5 *( ap(l) +ap(l+1))
256       bps(l) =  0.5 *( bp(l) +bp(l+1))
257      ENDDO
258
259      DO ilevel=1,nlevel
260        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
261      ENDDO
262
263      DO ilayer=1,nlayer
264        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
265        pk(ilayer)  =cpp*(play(ilayer)/preff)**rcp
266c       write(120,*) ilayer,plev(ilayer),play(ilayer)
267      ENDDO
268c     write(120,*) nlevel,plev(nlevel)
269c     stop
270     
271      pks=cpp*(psurf/preff)**rcp
272
273c  profil de temperature et altitude au premier appel
274c  --------------------------------------------------
275
276c modif par rapport a Mars:
277c   on envoie dz/T=-log(play/psurf)*r/g dans profile
278      tmp1(0)=0.0
279      tmp1(1)= -log(play(1)/psurf)*r/g
280      DO ilayer=2,nlayer
281        tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
282      ENDDO
283      call profile(unit,nlayer+1,tmp1,tmp2)
284      CLOSE(unit)
285
286      print*,"               Pression        Altitude     Temperature"
287      ilayer=1
288      tsurf(1)=tmp2(0)
289       temp(1)=tmp2(1)
290       zlay(1)=tmp2(1)*tmp1(1)
291      print*,"           0",tsurf(1)
292      print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
293      DO ilayer=2,nlayer
294        temp(ilayer)=tmp2(ilayer)
295        zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer)
296        print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
297      ENDDO
298     
299c     temperature du sous-sol
300c     ~~~~~~~~~~~~~~~~~~~~~~~
301      DO isoil=1,nsoil
302         tsoil(isoil)=tsurf(1)
303      ENDDO
304
305c    Initialisation des traceurs
306c    ---------------------------
307
308      DO iq=1,nqtot
309        DO ilayer=1,nlayer
310           q(ilayer,iq) = 0.
311        ENDDO
312      ENDDO
313
314c    Initialisation des parametres d'oro
315c    -----------------------------------
316
317      zmea(1) = 0.
318      zstd(1) = 0.
319      zsig(1) = 0.
320      zgam(1) = 0.
321      zthe(1) = 0.
322      zpic(1) = 0.
323      zval(1) = 0.
324
325c  Ecriture de "startphy.nc"
326c  -------------------------
327c  (Ce fichier sera aussitot relu au premier
328c   appel de "physiq", mais il est necessaire pour passer
329c   les variables purement physiques a "physiq"...
330
331      solsw(1)    = 0.
332      sollwdown(1)= 0.
333      dlw(1)      = 0.
334      radsol(1)   = 0.
335     
336      radpas      = NINT(1.*day_step/nbapp_rad)
337      soil_model  = .true.
338
339      call phyredem("startphy.nc  ",
340     .              lati,long,
341     .              tsurf,tsoil,albedo,
342     .              solsw,sollwdown,dlw,radsol,
343     .    zmea, zstd, zsig, zgam, zthe, zpic, zval,
344     .              temp)
345
346c=======================================================================
347c  BOUCLE TEMPORELLE DU MODELE 1D
348c=======================================================================
349c
350      firstcall=.true.
351      lastcall=.false.
352
353      DO idt=1,ndt
354        IF (idt.eq.ndt) then
355         lastcall=.true.
356c toujours nulle dans le cas de Venus, pour l'instant...
357         zls = 0.0
358c        write(103,*) 'Ls=',zls*180./pi
359c        write(103,*) 'Lat=', lati(1)
360c        write(103,*) 'RunEnd - Atmos. Temp. File'
361c        write(103,*) 'RunEnd - Atmos. Temp. File'
362c        write(104,*) 'Ls=',zls*180./pi
363c        write(104,*) 'Lat=', lati(1)
364c        write(104,*) 'RunEnd - Atmos. Temp. File'
365        ENDIF
366
367c    calcul du geopotentiel
368c     ~~~~~~~~~~~~~~~~~~~~~
369! ADAPTATION GCM POUR CP(T)
370      DO ilayer=1,nlayer
371        s(ilayer)=(play(ilayer)/psurf)**rcp
372      ENDDO
373      phi(1)=cpp*temp(1)*(1.E+0-s(1))
374      DO ilayer=2,nlayer
375         phi(ilayer)=phi(ilayer-1)+
376     &     cpp*(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))*0.5
377     &        *(s(ilayer-1)-s(ilayer))
378
379      ENDDO
380
381c       appel de la physique
382c       --------------------
383
384      CALL physiq (1,llm,nqtot,
385     ,     firstcall,lastcall,
386     ,     day,time,dtphys,
387     ,     plev,play,pk,phi,phisfi,
388     ,     presnivs,
389     ,     u,v,temp,q, 
390     ,     w,
391C - sorties
392     s     du,dv,dtemp,dq,dpsurf)
393
394c     print*,"DT APRES PHYSIQ=",day,time
395c     print*,dtemp
396c     print*,temp
397c     print*," "
398c     stop
399
400c       evolution du vent : modele 1D
401c       -----------------------------
402 
403c       la physique calcule les derivees temporelles de u et v.
404c       Pas de coriolis
405          DO ilayer=1,nlayer
406             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
407             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
408          ENDDO
409c     
410c       Calcul du temps au pas de temps suivant
411c       ---------------------------------------
412        firstcall=.false.
413        time=time+dtphys/daysec
414        IF (time.gt.1.E+0) then
415            time=time-1.E+0
416            day=day+1
417        ENDIF
418
419c       calcul des vitesses et temperature au pas de temps suivant
420c       ----------------------------------------------------------
421
422        DO ilayer=1,nlayer
423           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
424           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
425           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
426        ENDDO
427
428c       calcul des pressions au pas de temps suivant
429c       ----------------------------------------------------------
430
431           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
432           DO ilevel=1,nlevel
433             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
434           ENDDO
435           DO ilayer=1,nlayer
436             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
437           ENDDO
438
439      ENDDO   ! fin de la boucle temporelle
440
441c    ========================================================
442c    GESTION DES SORTIE
443c    ========================================================
444
445        print*,"Temperature finale:"
446        print*,temp
447       
448c stabilite
449      DO ilayer=1,nlayer
450        zlay(ilayer) = phi(ilayer)/g/1000.  !en km
451      ENDDO
452      DO ilayer=2,nlayer
453        tmp1(ilayer) =
454     .    (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1))
455     .   + 1000.*g/cpp
456      ENDDO
457
458      OPEN(11,file='profile.new')
459      write (11,*) tsurf
460      DO ilayer=1,nlayer
461        write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer)
462      ENDDO
463
464c    ========================================================
465      END
466 
467c***********************************************************************
468c***********************************************************************
469
470#include "../dyn3d/disvert_noterre.F"
471#include "../dyn3d/abort_gcm.F"
472#include "../dyn3d/dump2d.F"
473#include "../dyn3d/cpdet.F"
474
Note: See TracBrowser for help on using the repository browser.