source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/aaam_bud.F90 @ 2404

Last change on this file since 2404 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.8 KB
Line 
1
2! $Id: aaam_bud.F90 1992 2014-03-05 13:19:12Z fhourdin $
3
4SUBROUTINE aaam_bud(iam, nlon, nlev, rjour, rsec, rea, rg, ome, plat, plon, &
5    phis, dragu, liftu, phyu, dragv, liftv, phyv, p, u, v, aam, torsfc)
6
7  USE dimphy
8  IMPLICIT NONE
9  ! ======================================================================
10  ! Auteur(s): F.Lott (LMD/CNRS) date: 20031020
11  ! Object: Compute different terms of the axial AAAM Budget.
12  ! No outputs, every AAM quantities are written on the IAM
13  ! File.
14
15  ! Modif : I.Musat (LMD/CNRS) date : 20041020
16  ! Outputs : axial components of wind AAM "aam" and total surface torque
17  ! "torsfc",
18  ! but no write in the iam file.
19
20  ! WARNING: Only valid for regular rectangular grids.
21  ! REMARK: CALL DANS PHYSIQ AFTER lift_noro:
22  ! CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
23  ! C               ra,rg,romega,
24  ! C               rlat,rlon,pphis,
25  ! C               zustrdr,zustrli,zustrph,
26  ! C               zvstrdr,zvstrli,zvstrph,
27  ! C               paprs,u,v)
28
29  ! ======================================================================
30  ! Explicit Arguments:
31  ! ==================
32  ! iam-----input-I-File number where AAMs and torques are written
33  ! It is a formatted file that has been opened
34  ! in physiq.F
35  ! nlon----input-I-Total number of horizontal points that get into physics
36  ! nlev----input-I-Number of vertical levels
37  ! rjour        -R-Jour compte depuis le debut de la simu (run.def)
38  ! rsec         -R-Seconde de la journee
39  ! rea          -R-Earth radius
40  ! rg           -R-gravity constant
41  ! ome          -R-Earth rotation rate
42  ! plat ---input-R-Latitude en degres
43  ! plon ---input-R-Longitude en degres
44  ! phis ---input-R-Geopotential at the ground
45  ! dragu---input-R-orodrag stress (zonal)
46  ! liftu---input-R-orolift stress (zonal)
47  ! phyu----input-R-Stress total de la physique (zonal)
48  ! dragv---input-R-orodrag stress (Meridional)
49  ! liftv---input-R-orolift stress (Meridional)
50  ! phyv----input-R-Stress total de la physique (Meridional)
51  ! p-------input-R-Pressure (Pa) at model half levels
52  ! u-------input-R-Horizontal wind (m/s)
53  ! v-------input-R-Meridional wind (m/s)
54  ! aam-----output-R-Axial Wind AAM (=raam(3))
55  ! torsfc--output-R-Total surface torque (=tmou(3)+tsso(3)+tbls(3))
56
57  ! Implicit Arguments:
58  ! ===================
59
60  ! iim--common-I: Number of longitude intervals
61  ! jjm--common-I: Number of latitude intervals
62  ! klon-common-I: Number of points seen by the physics
63  ! iim*(jjm-1)+2 for instance
64  ! klev-common-I: Number of vertical layers
65  ! ======================================================================
66  ! Local Variables:
67  ! ================
68  ! dlat-----R: Latitude increment (Radians)
69  ! dlon-----R: Longitude increment (Radians)
70  ! raam  ---R: Wind AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
71  ! oaam  ---R: Mass AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
72  ! tmou-----R: Resolved Mountain torque (3 components)
73  ! tsso-----R: Parameterised Moutain drag torque (3 components)
74  ! tbls-----R: Parameterised Boundary layer torque (3 components)
75
76  ! LOCAL ARRAY:
77  ! ===========
78  ! zs    ---R: Topographic height
79  ! ps    ---R: Surface Pressure
80  ! ub    ---R: Barotropic wind zonal
81  ! vb    ---R: Barotropic wind meridional
82  ! zlat  ---R: Latitude in radians
83  ! zlon  ---R: Longitude in radians
84  ! ======================================================================
85
86  include "dimensions.h"
87  ! cc#include "dimphy.h"
88
89  ! ARGUMENTS
90
91  INTEGER iam, nlon, nlev
92  REAL, INTENT (IN) :: rjour, rsec, rea, rg, ome
93  REAL plat(nlon), plon(nlon), phis(nlon)
94  REAL dragu(nlon), liftu(nlon), phyu(nlon)
95  REAL dragv(nlon), liftv(nlon), phyv(nlon)
96  REAL p(nlon, nlev+1), u(nlon, nlev), v(nlon, nlev)
97
98  ! Variables locales:
99
100  INTEGER i, j, k, l
101  REAL xpi, hadley, hadday
102  REAL dlat, dlon
103  REAL raam(3), oaam(3), tmou(3), tsso(3), tbls(3)
104  INTEGER iax
105  ! IM ajout aam, torsfc
106  ! aam = composante axiale du Wind AAM raam
107  ! torsfc = composante axiale de (tmou+tsso+tbls)
108  REAL aam, torsfc
109
110  REAL zs(801, 401), ps(801, 401)
111  REAL ub(801, 401), vb(801, 401)
112  REAL ssou(801, 401), ssov(801, 401)
113  REAL blsu(801, 401), blsv(801, 401)
114  REAL zlon(801), zlat(401)
115
116  CHARACTER (LEN=20) :: modname = 'aaam_bud'
117  CHARACTER (LEN=80) :: abort_message
118
119
120
121  ! PUT AAM QUANTITIES AT ZERO:
122
123  IF (iim+1>801 .OR. jjm+1>401) THEN
124    abort_message = 'Pb de dimension dans aaam_bud'
125    CALL abort_gcm(modname, abort_message, 1)
126  END IF
127
128  xpi = acos(-1.)
129  hadley = 1.E18
130  hadday = 1.E18*24.*3600.
131  dlat = xpi/real(jjm)
132  dlon = 2.*xpi/real(iim)
133
134  DO iax = 1, 3
135    oaam(iax) = 0.
136    raam(iax) = 0.
137    tmou(iax) = 0.
138    tsso(iax) = 0.
139    tbls(iax) = 0.
140  END DO
141
142  ! MOUNTAIN HEIGHT, PRESSURE AND BAROTROPIC WIND:
143
144  ! North pole values (j=1):
145
146  l = 1
147
148  ub(1, 1) = 0.
149  vb(1, 1) = 0.
150  DO k = 1, nlev
151    ub(1, 1) = ub(1, 1) + u(l, k)*(p(l,k)-p(l,k+1))/rg
152    vb(1, 1) = vb(1, 1) + v(l, k)*(p(l,k)-p(l,k+1))/rg
153  END DO
154
155  zlat(1) = plat(l)*xpi/180.
156
157  DO i = 1, iim + 1
158
159    zs(i, 1) = phis(l)/rg
160    ps(i, 1) = p(l, 1)
161    ub(i, 1) = ub(1, 1)
162    vb(i, 1) = vb(1, 1)
163    ssou(i, 1) = dragu(l) + liftu(l)
164    ssov(i, 1) = dragv(l) + liftv(l)
165    blsu(i, 1) = phyu(l) - dragu(l) - liftu(l)
166    blsv(i, 1) = phyv(l) - dragv(l) - liftv(l)
167
168  END DO
169
170
171  DO j = 2, jjm
172
173    ! Values at Greenwich (Periodicity)
174
175    zs(iim+1, j) = phis(l+1)/rg
176    ps(iim+1, j) = p(l+1, 1)
177    ssou(iim+1, j) = dragu(l+1) + liftu(l+1)
178    ssov(iim+1, j) = dragv(l+1) + liftv(l+1)
179    blsu(iim+1, j) = phyu(l+1) - dragu(l+1) - liftu(l+1)
180    blsv(iim+1, j) = phyv(l+1) - dragv(l+1) - liftv(l+1)
181    zlon(iim+1) = -plon(l+1)*xpi/180.
182    zlat(j) = plat(l+1)*xpi/180.
183
184    ub(iim+1, j) = 0.
185    vb(iim+1, j) = 0.
186    DO k = 1, nlev
187      ub(iim+1, j) = ub(iim+1, j) + u(l+1, k)*(p(l+1,k)-p(l+1,k+1))/rg
188      vb(iim+1, j) = vb(iim+1, j) + v(l+1, k)*(p(l+1,k)-p(l+1,k+1))/rg
189    END DO
190
191
192    DO i = 1, iim
193
194      l = l + 1
195      zs(i, j) = phis(l)/rg
196      ps(i, j) = p(l, 1)
197      ssou(i, j) = dragu(l) + liftu(l)
198      ssov(i, j) = dragv(l) + liftv(l)
199      blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
200      blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
201      zlon(i) = plon(l)*xpi/180.
202
203      ub(i, j) = 0.
204      vb(i, j) = 0.
205      DO k = 1, nlev
206        ub(i, j) = ub(i, j) + u(l, k)*(p(l,k)-p(l,k+1))/rg
207        vb(i, j) = vb(i, j) + v(l, k)*(p(l,k)-p(l,k+1))/rg
208      END DO
209
210    END DO
211
212  END DO
213
214
215  ! South Pole
216
217  IF (jjm>1) THEN
218    l = l + 1
219    ub(1, jjm+1) = 0.
220    vb(1, jjm+1) = 0.
221    DO k = 1, nlev
222      ub(1, jjm+1) = ub(1, jjm+1) + u(l, k)*(p(l,k)-p(l,k+1))/rg
223      vb(1, jjm+1) = vb(1, jjm+1) + v(l, k)*(p(l,k)-p(l,k+1))/rg
224    END DO
225    zlat(jjm+1) = plat(l)*xpi/180.
226
227    DO i = 1, iim + 1
228      zs(i, jjm+1) = phis(l)/rg
229      ps(i, jjm+1) = p(l, 1)
230      ssou(i, jjm+1) = dragu(l) + liftu(l)
231      ssov(i, jjm+1) = dragv(l) + liftv(l)
232      blsu(i, jjm+1) = phyu(l) - dragu(l) - liftu(l)
233      blsv(i, jjm+1) = phyv(l) - dragv(l) - liftv(l)
234      ub(i, jjm+1) = ub(1, jjm+1)
235      vb(i, jjm+1) = vb(1, jjm+1)
236    END DO
237  END IF
238
239
240  ! MOMENT ANGULAIRE
241
242  DO j = 1, jjm
243    DO i = 1, iim
244
245      raam(1) = raam(1) - rea**3*dlon*dlat*0.5*(cos(zlon(i))*sin(zlat(j))*cos &
246        (zlat(j))*ub(i,j)+cos(zlon(i))*sin(zlat(j+1))*cos(zlat(j+ &
247        1))*ub(i,j+1)) + rea**3*dlon*dlat*0.5*(sin(zlon(i))*cos(zlat(j))*vb(i &
248        ,j)+sin(zlon(i))*cos(zlat(j+1))*vb(i,j+1))
249
250      oaam(1) = oaam(1) - ome*rea**4*dlon*dlat/rg*0.5*(cos(zlon(i))*cos(zlat( &
251        j))**2*sin(zlat(j))*ps(i,j)+cos(zlon(i))*cos(zlat(j+ &
252        1))**2*sin(zlat(j+1))*ps(i,j+1))
253
254      raam(2) = raam(2) - rea**3*dlon*dlat*0.5*(sin(zlon(i))*sin(zlat(j))*cos &
255        (zlat(j))*ub(i,j)+sin(zlon(i))*sin(zlat(j+1))*cos(zlat(j+ &
256        1))*ub(i,j+1)) - rea**3*dlon*dlat*0.5*(cos(zlon(i))*cos(zlat(j))*vb(i &
257        ,j)+cos(zlon(i))*cos(zlat(j+1))*vb(i,j+1))
258
259      oaam(2) = oaam(2) - ome*rea**4*dlon*dlat/rg*0.5*(sin(zlon(i))*cos(zlat( &
260        j))**2*sin(zlat(j))*ps(i,j)+sin(zlon(i))*cos(zlat(j+ &
261        1))**2*sin(zlat(j+1))*ps(i,j+1))
262
263      raam(3) = raam(3) + rea**3*dlon*dlat*0.5*(cos(zlat(j))**2*ub(i,j)+cos( &
264        zlat(j+1))**2*ub(i,j+1))
265
266      oaam(3) = oaam(3) + ome*rea**4*dlon*dlat/rg*0.5*(cos(zlat(j))**3*ps(i,j &
267        )+cos(zlat(j+1))**3*ps(i,j+1))
268
269    END DO
270  END DO
271
272
273  ! COUPLE DES MONTAGNES:
274
275
276  DO j = 1, jjm
277    DO i = 1, iim
278      tmou(1) = tmou(1) - rea**2*dlon*0.5*sin(zlon(i))*(zs(i,j)-zs(i,j+1))*( &
279        cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))
280      tmou(2) = tmou(2) + rea**2*dlon*0.5*cos(zlon(i))*(zs(i,j)-zs(i,j+1))*( &
281        cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))
282    END DO
283  END DO
284
285  DO j = 2, jjm
286    DO i = 1, iim
287      tmou(1) = tmou(1) + rea**2*dlat*0.5*sin(zlat(j))*(zs(i+1,j)-zs(i,j))*( &
288        cos(zlon(i+1))*ps(i+1,j)+cos(zlon(i))*ps(i,j))
289      tmou(2) = tmou(2) + rea**2*dlat*0.5*sin(zlat(j))*(zs(i+1,j)-zs(i,j))*( &
290        sin(zlon(i+1))*ps(i+1,j)+sin(zlon(i))*ps(i,j))
291      tmou(3) = tmou(3) - rea**2*dlat*0.5*cos(zlat(j))*(zs(i+1,j)-zs(i,j))*( &
292        ps(i+1,j)+ps(i,j))
293    END DO
294  END DO
295
296
297  ! COUPLES DES DIFFERENTES FRICTION AU SOL:
298
299  l = 1
300  DO j = 2, jjm
301    DO i = 1, iim
302      l = l + 1
303      tsso(1) = tsso(1) - rea**3*cos(zlat(j))*dlon*dlat*ssou(i, j)*sin(zlat(j &
304        ))*cos(zlon(i)) + rea**3*cos(zlat(j))*dlon*dlat*ssov(i, j)*sin(zlon(i &
305        ))
306
307      tsso(2) = tsso(2) - rea**3*cos(zlat(j))*dlon*dlat*ssou(i, j)*sin(zlat(j &
308        ))*sin(zlon(i)) - rea**3*cos(zlat(j))*dlon*dlat*ssov(i, j)*cos(zlon(i &
309        ))
310
311      tsso(3) = tsso(3) + rea**3*cos(zlat(j))*dlon*dlat*ssou(i, j)*cos(zlat(j &
312        ))
313
314      tbls(1) = tbls(1) - rea**3*cos(zlat(j))*dlon*dlat*blsu(i, j)*sin(zlat(j &
315        ))*cos(zlon(i)) + rea**3*cos(zlat(j))*dlon*dlat*blsv(i, j)*sin(zlon(i &
316        ))
317
318      tbls(2) = tbls(2) - rea**3*cos(zlat(j))*dlon*dlat*blsu(i, j)*sin(zlat(j &
319        ))*sin(zlon(i)) - rea**3*cos(zlat(j))*dlon*dlat*blsv(i, j)*cos(zlon(i &
320        ))
321
322      tbls(3) = tbls(3) + rea**3*cos(zlat(j))*dlon*dlat*blsu(i, j)*cos(zlat(j &
323        ))
324
325    END DO
326  END DO
327
328
329  ! write(*,*) 'AAM',rsec,
330  ! write(*,*) 'AAM',rjour+rsec/86400.,
331  ! c      raam(3)/hadday,oaam(3)/hadday,
332  ! c      tmou(3)/hadley,tsso(3)/hadley,tbls(3)/hadley
333
334  ! write(iam,100)rjour+rsec/86400.,
335  ! c      raam(1)/hadday,oaam(1)/hadday,
336  ! c      tmou(1)/hadley,tsso(1)/hadley,tbls(1)/hadley,
337  ! c      raam(2)/hadday,oaam(2)/hadday,
338  ! c      tmou(2)/hadley,tsso(2)/hadley,tbls(2)/hadley,
339  ! c      raam(3)/hadday,oaam(3)/hadday,
340  ! c      tmou(3)/hadley,tsso(3)/hadley,tbls(3)/hadley
341100 FORMAT (F12.5, 15(1X,F12.5))
342
343  ! write(iam+1,*)((zs(i,j),i=1,iim),j=1,jjm+1)
344  ! write(iam+1,*)((ps(i,j),i=1,iim),j=1,jjm+1)
345  ! write(iam+1,*)((ub(i,j),i=1,iim),j=1,jjm+1)
346  ! write(iam+1,*)((vb(i,j),i=1,iim),j=1,jjm+1)
347  ! write(iam+1,*)((ssou(i,j),i=1,iim),j=1,jjm+1)
348  ! write(iam+1,*)((ssov(i,j),i=1,iim),j=1,jjm+1)
349  ! write(iam+1,*)((blsu(i,j),i=1,iim),j=1,jjm+1)
350  ! write(iam+1,*)((blsv(i,j),i=1,iim),j=1,jjm+1)
351
352  aam = raam(3)
353  torsfc = tmou(3) + tsso(3) + tbls(3)
354
355  RETURN
356END SUBROUTINE aaam_bud
Note: See TracBrowser for help on using the repository browser.