source: LMDZ5/branches/testing/libf/phylmd/aaam_bud.F90 @ 3078

Last change on this file since 3078 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch

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