source: trunk/LMDZ.MARS/libf/phymars/dyn1d/profile.F @ 3027

Last change on this file since 3027 was 2274, checked in by emillour, 5 years ago

Mars GCM:
Some code cleanup:

  • move profile.F to dyn1d since it is only used by the 1D model
  • remove scatter.F, which is not used (and moreover a synonym of the scatter routine from the "mod_phys_lmdz_para" module)
  • remove obsolete "multipl.F" routine, replace the calls to it in vdifc_mod by modern Fortran syntax.

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