/[lmdze]/trunk/Sources/phylmd/Orography/orodrag.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/Orography/orodrag.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 178 - (hide annotations)
Fri Mar 11 18:47:26 2016 UTC (8 years, 3 months ago) by guez
File size: 6887 byte(s) 
Moved variables date0, deltat, datasz_max, ncvar_ids, point, buff_pos,
buffer, regular from module histcom_var to modules where they are
defined.

Removed procedure ioipslmpp, useless for a sequential program.

Added argument datasz_max to histwrite_real (to avoid circular
dependency with histwrite).

Removed useless variables and computations everywhere.

Changed real litteral constants from default kind to double precision
in lwb, lwu, lwvn, sw1s, swtt, swtt1, swu.

Removed unused arguments: paer of sw, sw1s, sw2s, swclr; pcldsw of
sw1s, sw2s; pdsig, prayl of swr; co2_ppm of clmain, clqh; tsol of
transp_lay; nsrf of screenp; kcrit and kknu of gwstress; pstd of
orosetup.

Added output of relative humidity.
1 guez 150 SUBROUTINE orodrag(nlon,nlev,ktest,ptsphy,paphm1,papm1,pgeom1, &
2 guez 23 ptm1,pum1,pvm1,pmea,pstd,psig,pgamma,ptheta,ppic,pval &
3     ,pulow,pvlow,pvom,pvol,pte)
4    
5     USE dimens_m
6     USE dimphy
7 guez 178 use gwstress_m, only: gwstress
8 guez 38 USE suphec_m
9 guez 23 USE yoegwd
10 guez 54 use gwprofil_m, only: gwprofil
11 guez 178 use orosetup_m, only: orosetup
12 guez 23 IMPLICIT NONE
13    
14    
15    
16     !**** *gwdrag* - does the gravity wave parametrization.
17    
18     ! purpose.
19     ! --------
20    
21     ! this routine computes the physical tendencies of the
22     ! prognostic variables u,v and t due to vertical transports by
23     ! subgridscale orographically excited gravity waves
24    
25     !** interface.
26     ! ----------
27     ! called from *callpar*.
28    
29     ! the routine takes its input from the long-term storage:
30     ! u,v,t and p at t-1.
31    
32     ! explicit arguments :
33     ! --------------------
34     ! ==== inputs ===
35     ! ==== outputs ===
36    
37     ! implicit arguments : none
38     ! --------------------
39    
40     ! implicit logical (l)
41    
42     ! method.
43     ! -------
44    
45     ! reference.
46     ! ----------
47    
48     ! author.
49     ! -------
50     ! m.miller + b.ritter e.c.m.w.f. 15/06/86.
51    
52     ! f.lott + m. miller e.c.m.w.f. 22/11/94
53     !-----------------------------------------------------------------------
54    
55    
56     !-----------------------------------------------------------------------
57    
58     !* 0.1 arguments
59     ! ---------
60    
61    
62 guez 178 INTEGER nlon, nlev
63 guez 150 INTEGER jl, ilevp1, jk, ji
64 guez 23 REAL zdelp, ztemp, zforc, ztend
65     REAL rover, zb, zc, zconb, zabsv
66     REAL zzd1, ratio, zbet, zust, zvst, zdis
67     REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(klon), &
68     pvlow(klon)
69     REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev), pmea(nlon)
70     REAL, INTENT (IN) :: pstd(nlon)
71     REAL, INTENT (IN) :: psig(nlon)
72     REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &
73     pgeom1(nlon,nlev), papm1(nlon,nlev), paphm1(nlon,nlev+1)
74    
75 guez 150 INTEGER ktest(nlon)
76 guez 23 !-----------------------------------------------------------------------
77    
78     !* 0.2 local arrays
79     ! ------------
80 guez 150 INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &
81     iknu(klon), iknu2(klon), ikcrit(klon)
82 guez 23
83 guez 150 REAL ztau(klon,klev+1), zstab(klon,klev+1), &
84 guez 23 zvph(klon,klev+1), zrho(klon,klev+1), zri(klon,klev+1), &
85     zpsi(klon,klev+1), zzdep(klon,klev)
86 guez 178 REAL zdudt(klon), zdvdt(klon), zvidis(klon), &
87 guez 23 znu(klon), zd1(klon), zd2(klon), zdmod(klon)
88 guez 178 REAL ztmst
89 guez 23 REAL, INTENT (IN) :: ptsphy
90    
91     !------------------------------------------------------------------
92    
93     !* 1. initialization
94     ! --------------
95    
96     !* 1.1 computational constants
97     ! -----------------------
98    
99     ztmst = ptsphy
100     ! ------------------------------------------------------------------
101    
102     !* 1.3 check whether row contains point for printing
103     ! ---------------------------------------------
104    
105     !* 2. precompute basic state variables.
106     !* ---------- ----- ----- ----------
107     !* define low level wind, project winds in plane of
108     !* low level wind, determine sector in which to take
109     !* the variance and set indicator for critical levels.
110    
111    
112     CALL orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &
113 guez 178 papm1,pum1,pvm1,ptm1,pgeom1,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &
114 guez 23 pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)
115    
116    
117    
118     !***********************************************************
119    
120    
121     !* 3. compute low level stresses using subcritical and
122     !* supercritical forms.computes anisotropy coefficient
123     !* as measure of orographic twodimensionality.
124    
125 guez 178 CALL gwstress(nlon,nlev,ktest,ikenvh,zrho,zstab,zvph,pstd, &
126 guez 23 psig,pmea,ppic,ztau,pgeom1,zdmod)
127    
128    
129     !* 4. compute stress profile.
130     !* ------- ------ --------
131    
132 guez 150 CALL gwprofil(nlon,nlev,ktest,ikcrith,icrit,paphm1,zrho,zstab, &
133 guez 23 zvph,zri,ztau,zdmod,psig,pstd)
134    
135    
136     !* 5. compute tendencies.
137     !* -------------------
138    
139     ! explicit solution at all levels for the gravity wave
140     ! implicit solution for the blocked levels
141    
142     DO 510 jl = 1, klon
143     zvidis(jl) = 0.0
144     zdudt(jl) = 0.0
145     zdvdt(jl) = 0.0
146     510 CONTINUE
147    
148     ilevp1 = klev + 1
149    
150    
151     DO 524 jk = 1, klev
152    
153    
154     ! Modif vectorisation 02/04/2004
155     DO 523 ji = 1, klon
156     IF (ktest(ji)==1) THEN
157    
158     zdelp = paphm1(ji,jk+1) - paphm1(ji,jk)
159     ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
160     zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
161     zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
162    
163     ! controle des overshoots:
164    
165     zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) + 1.E-12
166     ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst + 1.E-12
167     rover = 0.25
168     IF (zforc>=rover*ztend) THEN
169     zdudt(ji) = rover*ztend/zforc*zdudt(ji)
170     zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
171     END IF
172    
173     ! fin du controle des overshoots
174    
175     IF (jk>=ikenvh(ji)) THEN
176     zb = 1.0 - 0.18*pgamma(ji) - 0.04*pgamma(ji)**2
177     zc = 0.48*pgamma(ji) + 0.3*pgamma(ji)**2
178     zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
179     zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
180     zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2
181     ratio = (cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji, &
182     jk))**2)/(pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
183     zbet = max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
184    
185     ! simplement oppose au vent
186    
187     zdudt(ji) = -pum1(ji,jk)/ztmst
188     zdvdt(ji) = -pvm1(ji,jk)/ztmst
189    
190     ! projection dans la direction de l'axe principal de l'orographie
191     !mod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
192     !mod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
193     !mod * *cos(ptheta(ji)*rpi/180.)/ztmst
194     !mod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
195     !mod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
196     !mod * *sin(ptheta(ji)*rpi/180.)/ztmst
197     zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
198     zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
199     END IF
200     pvom(ji,jk) = zdudt(ji)
201     pvol(ji,jk) = zdvdt(ji)
202     zust = pum1(ji,jk) + ztmst*zdudt(ji)
203     zvst = pvm1(ji,jk) + ztmst*zdvdt(ji)
204     zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
205     zvidis(ji) = zvidis(ji) + zdis*zdelp
206    
207     ! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
208    
209     pte(ji,jk) = 0.0
210    
211     END IF
212     523 CONTINUE
213    
214     524 CONTINUE
215    
216    
217     RETURN
218     END
  ViewVC Help
Powered by ViewVC 1.1.21