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

Last change on this file since 1198 was 1121, checked in by slebonnois, 11 years ago

SL: bugs in rcm1d (venus, titan)

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