source: trunk/LMDZ.GENERIC/libf/phystd/profile.F @ 1477

Last change on this file since 1477 was 1150, checked in by aslmd, 11 years ago

LMDZ.GENERIC. put surface temperature update under .not.nosurf flag. commented and added prints for ichoice=8. added a modification of iradia in rcm1d.def to be able to still use same .def as subsequent 3D runs for consistency
File size: 6.4 KB
Line 
1      SUBROUTINE profile(nlev,zkm,temp)
2! to use  'getin'
3      USE ioipsl_getincom
4      IMPLICIT NONE
5c=======================================================================
6c     Subroutine utilisee dans "rcm1d"
7c     pour l'initialisation du profil atmospherique
8c=======================================================================
9c
10c   differents profils d'atmospheres. T=f(z)
11c   entree:
12c     nlev    nombre de niveaux
13c     zkm     alititudes en km
14c     ichoice choix de l'atmosphere:
15c             1 Temperature constante
16c             2 profil Savidjari
17c             3 Lindner (profil polaire)
18c             4 Inversion pour Francois
19c             5 Seiff (moyen)
20c             6 T constante + perturbation gauss (level) (christophe 10/98)
21c             7 T constante + perturbation gauss   (km)  (christophe 10/98)
22c             8 Lecture du profile dans un fichier ASCII (profile)
23c     tref    temperature de reference
24c     isin    ajout d'une perturbation (isin=1)
25c     pic     pic perturbation gauss pour ichoice = 6 ou 7
26c     largeur largeur de la perturbation gauss pour ichoice = 6 ou 7
27c     hauteur hauteur de la perturbation gauss pour ichoice = 6 ou 7
28c
29c   sortie:
30c     temp    temperatures en K
31c     
32c=======================================================================
33c-----------------------------------------------------------------------
34c   declarations:
35c   -------------
36
37c   arguments:
38c   ----------
39
40       INTEGER nlev, unit
41       REAL zkm(nlev),temp(nlev)
42
43c   local:
44c   ------
45
46      INTEGER il,ichoice,nseiff,iseiff,isin,iter
47      REAL pi
48      PARAMETER(nseiff=37)
49      REAL tref,t1,t2,t3,ww
50      REAL tseiff(nseiff)
51      DATA tseiff/214.,213.8,213.4,212.4,209.3,205.,201.4,197.8,
52     $           194.6,191.4,188.2,185.2,182.5,180.,177.5,175,
53     $           172.5,170.,167.5,164.8,162.4,160.,158.,156.,
54     $           154.1,152.2,150.3,148.7,147.2,145.7,144.2,143.,
55     $           142.,141.,140,139.5,139./
56      REAL pic,largeur
57      REAL hauteur,tmp
58
59c-----------------------------------------------------------------------
60c   read input profile type:
61c--------------------------
62
63      ichoice=1 ! default value for ichoice
64      call getin("ichoice",ichoice)
65      tref=200 ! default value for tref
66      call getin("tref",tref)
67      isin=0 ! default value for isin (=0 means no perturbation)
68      call getin("isin",isin)
69      pic=26.522 ! default value for pic
70      call getin("pic",pic)
71      largeur=10 ! default value for largeur
72      call getin("largeur",largeur)
73      hauteur=30 ! default value for hauteur
74      call getin("hauteur",hauteur)
75
76c-----------------------------------------------------------------------
77c   ichoice=1 temperature constante:
78c   --------------------------------
79
80      IF(ichoice.EQ.1) THEN
81         DO il=1,nlev
82            temp(il)=tref
83         ENDDO
84
85c-----------------------------------------------------------------------
86c   ichoice=2 savidjari:
87c   --------------------
88
89      ELSE IF(ichoice.EQ.2) THEN
90         DO il=1,nlev
91            temp(il)=AMAX1(219.-2.5*zkm(il),140.)
92         ENDDO
93
94c-----------------------------------------------------------------------
95c   ichoice=3 Lindner:
96c   ------------------
97
98      ELSE IF(ichoice.EQ.3) THEN
99         DO il=1,nlev
100            IF(zkm(il).LT.2.5) THEN
101               temp(il)=150.+30.*zkm(il)/2.5
102            ELSE IF(zkm(il).LT.5.) THEN
103               temp(il)=180.
104            ELSE
105               temp(il)=AMAX1(180.-2.*(zkm(il)-5.),130.)
106            ENDIF
107         ENDDO
108
109c-----------------------------------------------------------------------
110c   ichoice=4 Inversion pour Francois:
111c   ------------------
112
113      ELSE IF(ichoice.EQ.4) THEN
114         DO il=1,nlev
115            IF(zkm(il).LT.20.) THEN
116               temp(il)=135.
117            ELSE
118               temp(il)=AMIN1(135.+5.*(zkm(il)-20.),200.)
119            ENDIF
120         ENDDO
121
122
123c-----------------------------------------------------------------------
124c   ichoice=5 Seiff:
125c   ----------------
126
127      ELSE IF(ichoice.EQ.5) THEN
128         DO il=1,nlev
129            iseiff=INT(zkm(il)/2.)+1
130            IF(iseiff.LT.nseiff-1) THEN
131               temp(il)=tseiff(iseiff)+(zkm(il)-2.*(iseiff-1))*
132     $         (tseiff(iseiff+1)-tseiff(iseiff))/2.
133            ELSE
134               temp(il)=tseiff(nseiff)
135            ENDIF
136         ENDDO
137c IF(ichoice.EQ.6) THEN
138c           DO iter=1,6
139c           t2=temp(1)
140c           t3=temp(2)
141c           DO il=2,nlev-1
142c              t1=t2
143c              t2=t3
144c              t3=temp(il+1)
145c              ww=tanh(zkm(il)/20.)
146c              ww=ww*ww*ww
147c              temp(il)=t2+ww*.5*(t1-2.*t2+t3)
148c           ENDDO
149c           ENDDO
150c        ENDIF
151
152c-----------------------------------------------------------------------
153c   ichoice=6
154c   ---------
155
156      ELSE IF(ichoice.EQ.6) THEN
157      DO il=1,nlev
158        tmp=il-pic
159        temp(il)=tref + hauteur*exp(-tmp*tmp/largeur/largeur)
160      ENDDO
161
162
163c-----------------------------------------------------------------------
164c   ichoice=7
165c   ---------
166
167      ELSE IF(ichoice.EQ.7) THEN
168      DO il=1,nlev
169        tmp=zkm(il)-pic
170        temp(il)=tref + hauteur*exp(-tmp*tmp*4/largeur/largeur)
171      ENDDO
172
173c-----------------------------------------------------------------------
174c   ichoice=8
175c   ---------
176
177      ! first value is surface temperature
178      ! then profile of atmospheric temperature
179      ELSE IF(ichoice.GE.8) THEN
180      OPEN(11,file='profile',status='old',form='formatted',err=101)
181      DO il=1,nlev
182        READ (11,*) temp(il)
183      ENDDO
184
185      GOTO 201
186101   STOP'fichier profile inexistant'
187201   CONTINUE
188      CLOSE(10)
189
190c-----------------------------------------------------------------------
191
192      ENDIF
193
194c-----------------------------------------------------------------------
195c   rajout eventuel d'une perturbation:
196c   -----------------------------------
197
198      IF(isin.EQ.1) THEN
199         pi=2.*ASIN(1.)
200         DO il=1,nlev
201c       if (nlev.EQ.501) then
202c         if (zkm(il).LE.70.5) then
203c       temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*(
204c    s      6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
205c         endif
206c       else
207        temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*(
208     s      6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
209c       endif
210         ENDDO
211      ENDIF
212
213
214c-----------------------------------------------------------------------
215c   Ecriture du profil de temperature dans un fichier profile.out
216c   -------------------------------------------------------------
217
218
219      OPEN(12,file='profile.out',form='formatted')
220          DO il=1,nlev
221            write(12,*) temp(il)
222          ENDDO
223      CLOSE(12)
224
225      RETURN
226      END
Note: See TracBrowser for help on using the repository browser.