/[lmdze]/trunk/libf/phylmd/Orography/orodrag.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/Orography/orodrag.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (hide annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
File size: 7668 byte(s) 
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

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