source: trunk/libf/phyvenus/testphys1d.F @ 97

Last change on this file since 97 was 97, checked in by slebonnois, 14 years ago

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