source: trunk/LMDZ.VENUS/libf/phyvenus/profile.F @ 3556

Last change on this file since 3556 was 894, checked in by slebonnois, 12 years ago

SL: small modifications on rcm1d (Venus, Titan)
File size: 5.9 KB
Line 
1      SUBROUTINE profile(unit,nlev,dzst,temp)
2      IMPLICIT NONE
3c=======================================================================
4c     Subroutine utilisee dans le modele 1-D  "rcm1d"
5c     pour l'initialisation du profil atmospherique
6c=======================================================================
7c
8c   VERSION VENUS
9c
10c   differents profils d'atmospheres. T=f(z)
11c   entree:
12c     unit    unite de lecture de "rcm1d.def"
13c     nlev    nombre de niveaux (nlev=llm+1, surf + couches 1 a llm)
14c     dzst    dz/T (avec dz = epaisseur de la couche en m)
15c     ichoice choix de l'atmosphere:
16c             1 Temperature constante
17c             2 profil VIRA lisse
18c             3
19c             4
20c             5
21c             6 T constante + perturbation gauss (level) (christophe 10/98)
22c             7 T constante + perturbation gauss   (km)  (christophe 10/98)
23c             8 Lecture du profile dans un fichier ASCII (profile)
24c     tref    temperature de reference
25c     isin    ajout d'une perturbation (isin=1)
26c     pic     pic perturbation gauss pour ichoice = 6 ou 7
27c     largeur largeur de la perturbation gauss pour ichoice = 6 ou 7
28c     hauteur hauteur de la perturbation gauss pour ichoice = 6 ou 7
29c
30c   sortie:
31c     temp    temperatures en K
32c     
33c=======================================================================
34c-----------------------------------------------------------------------
35c   declarations:
36c   -------------
37
38c   arguments:
39c   ----------
40
41       INTEGER nlev, unit
42       REAL dzst(nlev),temp(nlev)
43
44c   local:
45c   ------
46
47      INTEGER il,ichoice,isin,iter
48      REAL pi
49      REAL tref,t1,t2,t3,ww
50      REAL pic,largeur
51      REAL hauteur,tmp
52      REAL zkm(nlev)    ! altitude en km
53
54      isin = 0
55
56c-----------------------------------------------------------------------
57c   choix du profil:
58c   ----------------
59
60c la lecture se fait dans le rcm1d.def, ouvert par rcm1d.F
61      READ(unit,*)
62      READ(unit,*)
63      READ(unit,*)
64      READ(unit,*) ichoice
65      READ(unit,*) tref
66      READ(unit,*) isin
67      READ(unit,*) pic
68      READ(unit,*) largeur
69      READ(unit,*) hauteur
70
71c-----------------------------------------------------------------------
72c   ichoice=1 temperature constante:
73c   --------------------------------
74
75      IF(ichoice.EQ.1) THEN
76         temp(1) = tref
77         zkm(1)  = 0.0
78         DO il=2,nlev
79            temp(il)= tref
80            zkm(il) = zkm(il-1)+temp(il)*dzst(il)/1000.
81         ENDDO
82
83c-----------------------------------------------------------------------
84c   ichoice=2 VIRA lisse:
85c   ---------------------
86
87      ELSE IF(ichoice.EQ.2) THEN
88         temp(1) = 735.
89         zkm(1)  = 0.0
90         DO il=2,nlev
91            zkm(il) = zkm(il-1)+temp(il-1)*dzst(il)/1000. ! approx avec T(l-1)
92            if(zkm(il).lt.60.) then
93              temp(il)=735.-7.95*zkm(il)
94            else
95              temp(il)=AMAX1(258.-3.*(zkm(il)-60.),168.)
96            endif
97            zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
98         ENDDO
99
100c-----------------------------------------------------------------------
101c   ichoice=3 VIRA lisse - 135K:
102c   ----------------------------
103
104      ELSE IF(ichoice.EQ.3) THEN
105         temp(1) = 600.
106         zkm(1)  = 0.0
107         DO il=2,nlev
108            zkm(il) = zkm(il-1)+temp(il-1)*dzst(il)/1000. ! approx avec T(l-1)
109            temp(il)=AMAX1(600.-7.95*zkm(il),168.)
110            zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
111         ENDDO
112
113c-----------------------------------------------------------------------
114c   ichoice=4 :
115c   ------------------
116
117      ELSE IF(ichoice.EQ.4) THEN
118         print*,"Cas non defini..."
119         print*,"Stop dans profile.F"
120         STOP
121
122c-----------------------------------------------------------------------
123c   ichoice=5 :
124c   ----------------
125
126      ELSE IF(ichoice.EQ.5) THEN
127         print*,"Cas non defini..."
128         print*,"Stop dans profile.F"
129         STOP
130
131c-----------------------------------------------------------------------
132c   ichoice=6
133c   ---------
134
135      ELSE IF(ichoice.EQ.6) THEN
136      temp(1) = tref
137      zkm(1)  = 0.0
138      DO il=2,nlev
139        tmp=il-pic
140        temp(il)= tref + hauteur*exp(-tmp*tmp/largeur/largeur)
141        zkm(il) = zkm(il-1)+temp(il)*dzst(il)/1000.
142      ENDDO
143
144
145c-----------------------------------------------------------------------
146c   ichoice=7
147c   ---------
148
149      ELSE IF(ichoice.EQ.7) THEN
150      temp(1) = tref
151      zkm(1)  = 0.0
152      DO il=2,nlev
153        zkm(il) = zkm(il-1)+tref*dzst(il)/1000. ! approx
154        tmp=zkm(il)-pic
155        temp(il)= tref + hauteur*exp(-tmp*tmp*4/largeur/largeur)
156        zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
157      ENDDO
158
159c-----------------------------------------------------------------------
160c   ichoice=8
161c   ---------
162
163      ELSE IF(ichoice.GE.8) THEN
164      OPEN(11,file='profile',status='old',form='formatted',err=101)
165      DO il=1,nlev
166        READ (11,*) temp(il)
167      ENDDO
168      zkm(1) = 0.0
169      DO il=2,nlev
170        zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
171      ENDDO
172
173      GOTO 201
174101   STOP'fichier profile inexistant'
175201   CONTINUE
176      CLOSE(10)
177
178c-----------------------------------------------------------------------
179
180      ENDIF
181
182c-----------------------------------------------------------------------
183c   rajout eventuel d'une perturbation:
184c   -----------------------------------
185
186      IF(isin.EQ.1) THEN
187         pi=2.*ASIN(1.)
188         DO il=1,nlev
189        temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*(
190     s      6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
191         ENDDO
192      ENDIF
193
194
195c-----------------------------------------------------------------------
196c   Ecriture du profil de temperature dans un fichier profile.out
197c   -------------------------------------------------------------
198
199
200c     OPEN(12,file='profile.out')
201c         DO il=1,nlev
202c            write(12,*) temp(il)
203c           write(12,*) temp(il),zkm(il)
204c         ENDDO
205c     CLOSE(12)
206
207      RETURN
208      END
Note: See TracBrowser for help on using the repository browser.