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

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

All GCMs:
Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation (up to rev r2420 of LMDZ5)

  • all physics packages:
  • added module callphysiq_mod.F90 in dynphy_lonlat/phy* which contains the routine "call_physiq" which is called by calfis* and calls the physics. This way different "physiq" routine from different physics packages may be called: The calfis* routines now exposes all available fields that might be transmitted to physiq but which is actually send (ie: expected/needed by physiq) is decided in call_physiq.
  • turned "physiq.F[90]" into module "physiq_mod.F[90]" for better control of "physiq" arguments. for phyvenus/phytitan, extracted gr_fi_ecrit from physiq.F as gr_fi_ecrit.F90 (note that it can only work in serial).
  • misc:
  • updated wxios.F90 to keep up with LMDZ5 modifications.
  • dyn3d_common:
  • infotrac.F90 keep up with LMDZ5 modifications (cosmetics)
  • dyn3d:
  • gcm.F90: cosmetic cleanup.
  • leapfrog.F90: fix computation of date as function of itau.
  • dyn3dpar:
  • gcm.F: cosmetic cleanup.
  • leapfrog_p.F90: fix computation of date as function of itau.

NB: physics are given the date corresponding to the end of the
physics step.

  • dynphy_lonlat:
  • calfis.F : added computation of relative wind vorticity.
  • calfis_p.F: added computation of relative wind vorticity (input required by Earth physics)

EM
File size: 14.3 KB
Line 
1      PROGRAM rcm1d
2     
3      USE infotrac
4      use control_mod, only: planet_type, day_step
5      USE phys_state_var_mod
6      use chemparam_mod
7      USE comconst_mod, ONLY: cpp,t0_venus,nu_venus
8      use cpdet_mod, only: ini_cpdet
9      use moyzon_mod, only: tmoy
10      USE comvert_mod, ONLY: ap,bp,presnivs,pa,preff,nivsigs,nivsig,
11     .                       aps,bps,scaleheight,pseudoalt,
12     .                       disvert_type,pressure_exner
13      use conc, only: rho
14      USE iniphysiq_mod, ONLY: iniphysiq
15      USE mod_const_mpi, ONLY: comm_lmdz
16      USE physiq_mod, ONLY: physiq
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(1)   
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,
213     &            1,comm_lmdz,
214     &            daysec,day0,dtphys,
215     &            (/lati(1),0./),(/0./),
216     &            (/0.,0./),(/long(1),0./),
217     &            (/ (/area,0./),(/0.,0./) /),
218     &            (/cufi,0.,0.,0./),
219     &            (/cvfi,0./),
220     &            rad,g,r,cpp,1)
221
222      call ini_cpdet
223
224c   le geopotentiel au sol est inutile en 1D car tout est controle
225c   par la pression de surface --->
226      phisfi(1)=0.E+0
227
228c   Initialisation pour prendre en compte les vents en 1-D
229c   ------------------------------------------------------
230 
231c    vent geostrophique
232      PRINT *,'composante vers l est du vent geostrophique (U) ?'
233      READ(unit,*) gru
234      PRINT *,'composante vers le nord du vent geostrophique (V) ?'
235      READ(unit,*) grv
236
237c     Initialisation des vents  au premier pas de temps
238      DO ilayer=1,nlayer
239         u(ilayer)=gru
240         v(ilayer)=grv
241      ENDDO
242
243c  calcul des pressions et altitudes en utilisant les niveaux sigma
244c  ----------------------------------------------------------------
245
246c    Vertical Coordinates  (hybrids)
247c    """"""""""""""""""""
248      CALL  disvert_noterre
249     
250c     Calcul au milieu des couches : Vient de la version Mars
251c     WARNING : le choix de placer le milieu des couches au niveau de
252c     pression intermédiaire est arbitraire et pourrait etre modifié.
253c     C'est fait de la meme facon dans disvert
254
255      DO l = 1, llm
256       aps(l) =  0.5 *( ap(l) +ap(l+1))
257       bps(l) =  0.5 *( bp(l) +bp(l+1))
258      ENDDO
259
260      DO ilevel=1,nlevel
261        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
262      ENDDO
263
264      DO ilayer=1,nlayer
265        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
266        pk(ilayer)  =cpp*(play(ilayer)/preff)**rcp
267c       write(120,*) ilayer,plev(ilayer),play(ilayer)
268      ENDDO
269c     write(120,*) nlevel,plev(nlevel)
270c     stop
271     
272      pks=cpp*(psurf/preff)**rcp
273
274c  init des variables pour phyredem
275c  --------------------------------
276      call phys_state_var_init
277
278c  profil de temperature et altitude au premier appel
279c  --------------------------------------------------
280
281c modif par rapport a Mars:
282c   on envoie dz/T=-log(play/psurf)*r/g dans profile
283      tmp1(0)=0.0
284      tmp1(1)= -log(play(1)/psurf)*r/g
285      DO ilayer=2,nlayer
286        tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
287      ENDDO
288      call profile(unit,nlayer+1,tmp1,tmp2)
289      CLOSE(unit)
290
291      print*,"               Pression        Altitude     Temperature"
292      ilayer=1
293      ftsol(1)=tmp2(0)
294       temp(1)=tmp2(1)
295       zlay(1)=tmp2(1)*tmp1(1)
296      print*,"           0",ftsol(1)
297      print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
298      DO ilayer=2,nlayer
299        temp(ilayer)=tmp2(ilayer)
300        zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer)
301        print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
302      ENDDO
303
304      allocate(tmoy(llm))
305      tmoy(:)=temp(:)
306     
307c     temperature du sous-sol
308c     ~~~~~~~~~~~~~~~~~~~~~~~
309      DO isoil=1,nsoil
310         ftsoil(1,isoil)=ftsol(1)
311      ENDDO
312
313c    Initialisation des traceurs
314c    ---------------------------
315
316      DO iq=1,nqtot
317        DO ilayer=1,nlayer
318           q(ilayer,iq) = 0.
319        ENDDO
320      ENDDO
321
322c FULL CHEMISTRY !! AJOUTER INIT AURELIEN...
323C Faudrait lire les cles avant pour mettre ca en option....
324c ou alors mettre ca dans physiq
325
326c    Initialisation des parametres d'oro
327c    -----------------------------------
328
329      zmea(1) = 0.
330      zstd(1) = 0.
331      zsig(1) = 0.
332      zgam(1) = 0.
333      zthe(1) = 0.
334      zpic(1) = 0.
335      zval(1) = 0.
336
337c  Initialisation albedo
338c  ----------------------
339
340      falbe(1)=0.1
341
342c  Ecriture de "startphy.nc"
343c  -------------------------
344c  (Ce fichier sera aussitot relu au premier
345c   appel de "physiq", mais il est necessaire pour passer
346c   les variables purement physiques a "physiq"...
347
348      solsw(1)    = 0.
349      sollw(1)    = 0.
350      fder(1)      = 0.
351      radsol(1)   = 0.
352     
353      radpas      = NINT(1.*day_step/nbapp_rad)
354      soil_model  = .true.
355
356      call phyredem("startphy.nc")
357
358c  deallocation des variables phyredem
359c  -----------------------------------
360      call phys_state_var_end
361
362c=======================================================================
363c  BOUCLE TEMPORELLE DU MODELE 1D
364c=======================================================================
365c
366      firstcall=.true.
367      lastcall=.false.
368
369      DO idt=1,ndt
370        IF (idt.eq.ndt) then
371         lastcall=.true.
372c toujours nulle dans le cas de Venus, pour l'instant...
373         zls = 0.0
374c        write(103,*) 'Ls=',zls*180./pi
375c        write(103,*) 'Lat=', lati(1)
376c        write(103,*) 'RunEnd - Atmos. Temp. File'
377c        write(103,*) 'RunEnd - Atmos. Temp. File'
378c        write(104,*) 'Ls=',zls*180./pi
379c        write(104,*) 'Lat=', lati(1)
380c        write(104,*) 'RunEnd - Atmos. Temp. File'
381        ENDIF
382
383c    calcul du geopotentiel
384c     ~~~~~~~~~~~~~~~~~~~~~
385! ADAPTATION GCM POUR CP(T)
386      DO ilayer=1,nlayer
387        s(ilayer)=(play(ilayer)/psurf)**rcp
388      ENDDO
389      phi(1)=cpp*temp(1)*(1.E+0-s(1))
390      DO ilayer=2,nlayer
391         phi(ilayer)=phi(ilayer-1)+
392     &     cpp*(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))*0.5
393     &        *(s(ilayer-1)-s(ilayer))
394
395      ENDDO
396
397c       appel de la physique
398c       --------------------
399
400      CALL physiq (1,llm,nqtot,
401     ,     firstcall,lastcall,
402     ,     day,time,dtphys,
403     ,     plev,play,pk,phi,phisfi,
404     ,     presnivs,
405     ,     u,v,temp,q, 
406     ,     w,
407C - sorties
408     s     du,dv,dtemp,dq,dpsurf)
409
410c     calcul de rho
411       rho = 0.
412c     print*,rho
413
414
415c     print*,"DT APRES PHYSIQ=",day,time,dtime
416c     print*,dtemp
417c     print*,temp
418c     print*," "
419c     stop
420
421c       evolution du vent : modele 1D
422c       -----------------------------
423 
424c       la physique calcule les derivees temporelles de u et v.
425c       Pas de coriolis
426          DO ilayer=1,nlayer
427             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
428             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
429          ENDDO
430c     
431c       Calcul du temps au pas de temps suivant
432c       ---------------------------------------
433        firstcall=.false.
434        time=time+dtphys/daysec
435        IF (time.gt.1.E+0) then
436            time=time-1.E+0
437            day=day+1
438        ENDIF
439
440c       calcul des vitesses et temperature au pas de temps suivant
441c       ----------------------------------------------------------
442
443        DO ilayer=1,nlayer
444           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
445           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
446           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
447        ENDDO
448
449c       calcul des pressions au pas de temps suivant
450c       ----------------------------------------------------------
451
452           psurf=psurf+dtphys*dpsurf(1)   ! evolution de la pression de surface
453           DO ilevel=1,nlevel
454             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
455           ENDDO
456           DO ilayer=1,nlayer
457             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
458           ENDDO
459
460      ENDDO   ! fin de la boucle temporelle
461
462c    ========================================================
463c    GESTION DES SORTIE
464c    ========================================================
465
466        print*,"Temperature finale:"
467        print*,temp
468       
469c stabilite
470      DO ilayer=1,nlayer
471        zlay(ilayer) = phi(ilayer)/g/1000.  !en km
472      ENDDO
473      DO ilayer=2,nlayer
474        tmp1(ilayer) =
475     .    (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1))
476     .   + 1000.*g/cpp
477      ENDDO
478
479      OPEN(11,file='profile.new')
480      DO ilayer=1,nlayer
481        write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer)
482      ENDDO
483
484c    ========================================================
485      END
486 
487c***********************************************************************
488c***********************************************************************
489
490!#include "../dyn3d_common/disvert_noterre.F"
491!#include "../dyn3d/abort_gcm.F"
492
Note: See TracBrowser for help on using the repository browser.