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

Last change on this file since 1523 was 1523, checked in by emillour, 9 years ago

All models: More updates to make planetary codes (+Earth) setups converge.

  • in dyn3d_common:
  • convmas.F => convmas.F90
  • enercin.F => enercin.F90
  • flumass.F => flumass.F90
  • massbar.F => massbar.F90
  • tourpot.F => tourpot.F90
  • vitvert.F => vitvert.F90
  • in misc:
  • move "q_sat" from "dyn3d_common" to "misc" (in Earth model, it is also called by the physics)
  • move "write_field" from "dyn3d_common" to "misc"(may be called from physics or dynamics and depends on neither).
  • in phy_common:
  • move "write_field_phy" here since it may be called from any physics package)
  • add module "regular_lonlat_mod" to store global information on lon-lat grid
  • in dynlonlat_phylonlat/phy*:
  • turn "iniphysiq.F90" into module "iniphysiq_mod.F90" (and of course adapt gcm.F[90] and 1D models accordingly)

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