source: trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F @ 1048

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

SL: update pour divers details titan + quelques modifs arch et makelmdz

File size: 14.1 KB
Line 
1      PROGRAM rcm1d
2     
3      USE infotrac
4      use control_mod
5      use comgeomphy
6      use cpdet_mod, only: ini_cpdet
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 55 layers)
15c  "makelmdz -p titan -d 55 rcm1d"
16
17c It requires the files "rcm1d.def" "physiq.def" "traceur.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
25c Version TITAN a tester et verifier 
26c  - verifier pour Ls...
27c  - faire un profile.F ...
28
29#include "dimensions.h"
30#include "dimsoil.h"
31#include "comcstfi.h"
32#include "comvert.h"
33#include "netcdf.inc"
34#include "logic.h"
35#include "clesphys.h"
36#include "iniprint.h"
37#include "tabcontrol.h"
38
39c --------------------------------------------------------------
40c  Declarations
41c --------------------------------------------------------------
42c
43      INTEGER unit           ! unite de lecture de "rcm1d.def"
44      INTEGER unitstart      ! unite d'ecriture de "startphy.nc"
45      INTEGER nlayer,nlevel,nsoil,ndt
46      INTEGER ilayer,ilevel,isoil,idt,iq
47      LOGICAl firstcall,lastcall
48c
49      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
50      REAL day              ! date durant le run
51      REAL time             ! time (0<time<1 ; time=0.5 a midi)
52      REAL play(llm)   ! Pressure at the middle of the layers (Pa)
53      REAL plev(llm+1) ! intermediate pressure levels (pa)
54      REAL psurf,tsurf(1)     
55      REAL u(llm),v(llm)  ! zonal, meridional wind
56      REAL gru,grv   ! prescribed "geostrophic" background wind
57      REAL temp(llm)   ! temperature at the middle of the layers
58      REAL,allocatable :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
59      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
60      REAL zlay(llm)   ! altitude estimee dans les couches (km)
61      REAL long(1),lati(1),area(1)
62      REAL cufi(1),cvfi(1)
63      REAL phisfi(1),albedo(1)
64      REAL solsw(1),sollwdown(1),dlw(1),radsol(1)
65      REAL zmea(1), zstd(1)
66      REAL zsig(1), zgam(1), zthe(1)
67      REAL zpic(1), zval(1)
68
69c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
70      REAL du(llm),dv(llm),dtemp(llm)
71      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
72      REAL dpsurf   
73      REAL,allocatable :: dq(:,:)
74
75c   Various intermediate variables
76      REAL zls
77      REAL phi(llm),s(llm)
78      REAL pk(llm),pks, w(llm)
79      INTEGER l, ierr, aslun
80      REAL tmp1(0:llm),tmp2(0:llm),tmp3(0:llm)                       
81
82      character*2 str2
83
84c normalement dans dyn3d/comconst.h
85      COMMON/cpdetvenus/cppdyn,nu_venus,t0_venus
86      REAL cppdyn,nu_venus,t0_venus
87
88c=======================================================================
89
90c=======================================================================
91c INITIALISATION
92c=======================================================================
93
94      lunout = 6
95
96c ------------------------------------------------------
97c  Constantes prescrites ICI
98c ------------------------------------------------------
99
100      pi=2.E+0*asin(1.E+0)
101
102c     Constante de Titan
103c     -------------------
104      planet_type = "titan"
105      rad=2575000.               ! rayon de Venus (m)
106      daysec=1.37889e6           ! duree du sol (s) 
107      omeg=2*pi/daysec           ! vitesse de rotation (rad.s-1)
108      g= 1.35                    ! gravite (m.s-2)
109      mugaz=28.                  ! Masse molaire de l'atm (g.mol-1)
110      cpp=1.039e3
111      cppdyn=cpp
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 = 10.
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
178c Pression de surface sur la planete
179c ------------------------------------
180c
181      PRINT *,'pression au sol'
182      READ(unit,*) psurf
183      PRINT *,psurf
184c Pression de reference 
185      pa     =  1.e4
186      preff  = 1.4e5
187c     preff  = psurf
188 
189c  latitude/longitude
190c  -------------------
191      PRINT *,'latitude en degres ?'
192      READ(unit,*) lati(1)
193      PRINT *,lati(1)
194      long(1)=0.E+0
195
196c  Initialisation albedo
197c  ----------------------
198c ne sert pas ici...
199      albedo(1)=0.3
200c      PRINT *,'Albedo du sol nu ?'
201c      READ(unit,*) albedo(1)
202c      PRINT *,albedo(1)
203
204c   Initialisation speciales "physiq"
205c   ---------------------------------
206
207      CALL init_phys_lmdz(iim,jjm,llm,1,(/1/))
208      call initcomgeomphy
209      call infotrac_init
210      call ini_cpdet
211
212c   la surface de chaque maille est inutile en 1D --->
213      area(1)=1.E+0
214c de meme ?
215      cufi(1)=1.E+0
216      cvfi(1)=1.E+0
217
218c   le geopotentiel au sol est inutile en 1D car tout est controle
219c   par la pression de surface --->
220      phisfi(1)=0.E+0
221
222      CALL iniphysiq(1,llm,daysec,day0,dtphys,
223     .            lati,long,area,cufi,cvfi,rad,g,r,cpp)
224
225c   Initialisation pour prendre en compte les vents en 1-D
226c   ------------------------------------------------------
227 
228c    vent geostrophique
229      PRINT *,'composante vers l est du vent geostrophique (U) ?'
230      READ(unit,*) gru
231      PRINT *,'composante vers le nord du vent geostrophique (V) ?'
232      READ(unit,*) grv
233
234c     Initialisation des vents  au premier pas de temps
235      DO ilayer=1,nlayer
236         u(ilayer)=gru
237         v(ilayer)=grv
238      ENDDO
239
240c  calcul des pressions et altitudes en utilisant les niveaux sigma
241c  ----------------------------------------------------------------
242
243c    Vertical Coordinates  (hybrids)
244c    """"""""""""""""""""
245      CALL  disvert_noterre
246     
247c     Calcul au milieu des couches : Vient de la version Mars
248c     WARNING : le choix de placer le milieu des couches au niveau de
249c     pression intermédiaire est arbitraire et pourrait etre modifié.
250c     C'est fait de la meme facon dans disvert
251
252      DO l = 1, llm
253       aps(l) =  0.5 *( ap(l) +ap(l+1))
254       bps(l) =  0.5 *( bp(l) +bp(l+1))
255      ENDDO
256
257      DO ilevel=1,nlevel
258        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
259      ENDDO
260
261      DO ilayer=1,nlayer
262        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
263        pk(ilayer)  =cpp*(play(ilayer)/preff)**rcp
264c       write(120,*) ilayer,plev(ilayer),play(ilayer)
265      ENDDO
266c     write(120,*) nlevel,plev(nlevel)
267c     stop
268     
269      pks=cpp*(psurf/preff)**rcp
270
271c  profil de temperature et altitude au premier appel
272c  --------------------------------------------------
273
274c modif par rapport a Mars:
275c   on envoie dz/T=-log(play/psurf)*r/g dans profile
276      tmp1(0)=0.0
277      tmp1(1)= -log(play(1)/psurf)*r/g
278      DO ilayer=2,nlayer
279        tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
280      ENDDO
281      DO ilayer=0,nlayer
282        tmp2(ilayer)=plev(ilayer+1)
283      ENDDO
284      call profile(unit,nlayer+1,tmp1,tmp2,tmp3)
285      CLOSE(unit)
286
287      print*,"               Pression        Altitude     Temperature"
288      ilayer=1
289      tsurf(1)=tmp3(0)
290       temp(1)=tmp3(1)
291       zlay(1)=tmp3(1)*tmp1(1)
292      print*,"           0",tsurf(1)
293      print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
294      DO ilayer=2,nlayer
295        temp(ilayer)=tmp3(ilayer)
296        zlay(ilayer)=zlay(ilayer-1)+tmp3(ilayer)*tmp1(ilayer)
297        print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
298      ENDDO
299     
300c     temperature du sous-sol
301c     ~~~~~~~~~~~~~~~~~~~~~~~
302      DO isoil=1,nsoil
303         tsoil(isoil)=93.
304      ENDDO
305
306c    Initialisation des traceurs
307c    ---------------------------
308
309      DO iq=1,nqtot
310        DO ilayer=1,nlayer
311           q(ilayer,iq) = 0.
312        ENDDO
313      ENDDO
314
315c    Initialisation des parametres d'oro
316c    -----------------------------------
317
318      zmea(1) = 0.
319      zstd(1) = 0.
320      zsig(1) = 0.
321      zgam(1) = 0.
322      zthe(1) = 0.
323      zpic(1) = 0.
324      zval(1) = 0.
325
326c  Initialisation Ls
327c ------------------
328         zls=0.
329         print*,'Ls=',zls*180./pi
330
331c  Ecriture de "startphy.nc"
332c  -------------------------
333c  (Ce fichier sera aussitot relu au premier
334c   appel de "physiq", mais il est necessaire pour passer
335c   les variables purement physiques a "physiq"...
336
337      solsw(1)    = 0.
338      sollwdown(1)= 0.
339      dlw(1)      = 0.
340      radsol(1)   = 0.
341     
342      radpas      = NINT(1.*day_step/nbapp_rad)
343      soil_model  = .true.
344
345      call phyredem("startphy.nc  ",
346     .              lati,long,
347     .              tsurf,tsoil,albedo,
348     .              solsw,sollwdown,dlw,radsol,
349     .    zmea, zstd, zsig, zgam, zthe, zpic, zval,
350     .              temp)
351
352c=======================================================================
353c  BOUCLE TEMPORELLE DU MODELE 1D
354c=======================================================================
355c
356      firstcall=.true.
357      lastcall=.false.
358
359      DO idt=1,ndt
360        IF (idt.eq.ndt) then
361         lastcall=.true.
362c        write(103,*) 'Ls=',zls*180./pi
363c        write(103,*) 'Lat=', lati(1)
364c        write(103,*) 'RunEnd - Atmos. Temp. File'
365c        write(103,*) 'RunEnd - Atmos. Temp. File'
366c        write(104,*) 'Ls=',zls*180./pi
367c        write(104,*) 'Lat=', lati(1)
368c        write(104,*) 'RunEnd - Atmos. Temp. File'
369        ENDIF
370
371c    calcul du geopotentiel
372c     ~~~~~~~~~~~~~~~~~~~~~
373! ADAPTATION GCM POUR CP(T)
374      DO ilayer=1,nlayer
375        s(ilayer)=(play(ilayer)/psurf)**rcp
376      ENDDO
377      phi(1)=cpp*temp(1)*(1.E+0-s(1))
378      DO ilayer=2,nlayer
379         phi(ilayer)=phi(ilayer-1)+
380     &     cpp*(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))*0.5
381     &        *(s(ilayer-1)-s(ilayer))
382
383      ENDDO
384
385c       appel de la physique
386c       --------------------
387
388      CALL physiq (1,llm,nqtot,
389     ,     firstcall,lastcall,
390     ,     day,time,dtphys,
391     ,     plev,play,pk,phi,phisfi,
392     ,     presnivs,
393     ,     u,v,temp,q, 
394     ,     w,
395C - sorties
396     s     du,dv,dtemp,dq,dpsurf)
397
398c     print*,"DT APRES PHYSIQ=",day,time
399c     print*,dtemp
400c     print*,temp
401c     print*," "
402c     stop
403
404c       evolution du vent : modele 1D
405c       -----------------------------
406 
407c       la physique calcule les derivees temporelles de u et v.
408c       Pas de coriolis
409          DO ilayer=1,nlayer
410             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
411             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
412          ENDDO
413c     
414c       Calcul du temps au pas de temps suivant
415c       ---------------------------------------
416        firstcall=.false.
417        time=time+dtphys/daysec
418        IF (time.gt.1.E+0) then
419            time=time-1.E+0
420            day=day+1
421        ENDIF
422
423c       calcul des vitesses et temperature au pas de temps suivant
424c       ----------------------------------------------------------
425
426        DO ilayer=1,nlayer
427           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
428           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
429           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
430        ENDDO
431
432c       calcul des pressions au pas de temps suivant
433c       ----------------------------------------------------------
434
435           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
436           DO ilevel=1,nlevel
437             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
438           ENDDO
439           DO ilayer=1,nlayer
440             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
441           ENDDO
442
443      ENDDO   ! fin de la boucle temporelle
444
445c    ========================================================
446c    GESTION DES SORTIE
447c    ========================================================
448
449        print*,"Temperature finale:"
450        print*,temp
451       
452c stabilite
453      DO ilayer=1,nlayer
454        zlay(ilayer) = phi(ilayer)/g/1000.  !en km
455      ENDDO
456      DO ilayer=2,nlayer
457        tmp1(ilayer) =
458     .    (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1))
459     .   + 1000.*g/cpp
460      ENDDO
461
462      OPEN(11,file='profile.new')
463      write (11,*) tsurf
464      DO ilayer=1,nlayer
465        write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer)
466      ENDDO
467
468c    ========================================================
469      END
470 
471c***********************************************************************
472c***********************************************************************
473
474#include "../dyn3d/disvert_noterre.F"
475#include "../dyn3d/abort_gcm.F"
476!#include "../dyn3d/dump2d.F"
477
478c***********************************************************************
479      function ssum(n,sx,incx)
480c
481      IMPLICIT NONE
482c
483      integer n,incx,i,ix
484      real ssum,sx((n-1)*incx+1)
485c
486      ssum=0.
487      ix=1
488      do 10 i=1,n
489         ssum=ssum+sx(ix)
490         ix=ix+incx
49110    continue
492c
493      return
494      end
Note: See TracBrowser for help on using the repository browser.