source: LMDZ5/trunk/libf/phylmd/aaam_bud.F90 @ 4799

Last change on this file since 4799 was 2350, checked in by Ehouarn Millour, 9 years ago

Corrections for 1D case where nbp_lon=nbp_lat=1.
EM + MPL

  • 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
RevLine 
[1992]1
[1403]2! $Id: aaam_bud.F90 2350 2015-08-25 11:40:19Z nfevrier $
[644]3
[1992]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)
[644]6
[1992]7  USE dimphy
[2350]8  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo
[1992]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.
[1403]15
[1992]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.
[1403]20
[1992]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)
[1403]29
[1992]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))
[644]57
[1992]58  ! Implicit Arguments:
59  ! ===================
[644]60
[2346]61  ! nbp_lon--common-I: Number of longitude intervals
62  ! (nbp_lat-1)--common-I: Number of latitude intervals
[1992]63  ! klon-common-I: Number of points seen by the physics
[2346]64  ! nbp_lon*(nbp_lat-2)+2 for instance
[1992]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)
[644]76
[1992]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  ! ======================================================================
[644]86
[1992]87  ! ARGUMENTS
[644]88
[1992]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)
[644]95
[1992]96  ! Variables locales:
[644]97
[1992]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
[644]107
[1992]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)
[644]113
[1992]114  CHARACTER (LEN=20) :: modname = 'aaam_bud'
115  CHARACTER (LEN=80) :: abort_message
[644]116
117
118
[1992]119  ! PUT AAM QUANTITIES AT ZERO:
[644]120
[2346]121  IF (nbp_lon+1>801 .OR. nbp_lat>401) THEN
[1992]122    abort_message = 'Pb de dimension dans aaam_bud'
[2311]123    CALL abort_physic(modname, abort_message, 1)
[1992]124  END IF
[644]125
[1992]126  xpi = acos(-1.)
127  hadley = 1.E18
128  hadday = 1.E18*24.*3600.
[2350]129  IF(klon_glo.EQ.1) THEN
130    dlat = xpi
131  ELSE
132    dlat = xpi/real(nbp_lat-1)
133  ENDIF
[2346]134  dlon = 2.*xpi/real(nbp_lon)
[644]135
[1992]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
[644]143
[1992]144  ! MOUNTAIN HEIGHT, PRESSURE AND BAROTROPIC WIND:
[644]145
[1992]146  ! North pole values (j=1):
[644]147
[1992]148  l = 1
[644]149
[1992]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
[644]156
[1992]157  zlat(1) = plat(l)*xpi/180.
[644]158
[2346]159  DO i = 1, nbp_lon + 1
[644]160
[1992]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)
[644]169
[1992]170  END DO
[644]171
172
[2346]173  DO j = 2, nbp_lat-1
[644]174
[1992]175    ! Values at Greenwich (Periodicity)
[644]176
[2346]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.
[1992]184    zlat(j) = plat(l+1)*xpi/180.
[644]185
[2346]186    ub(nbp_lon+1, j) = 0.
187    vb(nbp_lon+1, j) = 0.
[1992]188    DO k = 1, nlev
[2346]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
[1992]191    END DO
[644]192
193
[2346]194    DO i = 1, nbp_lon
[644]195
[1992]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.
[644]204
[1992]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
[644]211
[1992]212    END DO
[644]213
[1992]214  END DO
[644]215
216
[1992]217  ! South Pole
[644]218
[2346]219  IF (nbp_lat-1>1) THEN
[1992]220    l = l + 1
[2346]221    ub(1, nbp_lat) = 0.
222    vb(1, nbp_lat) = 0.
[1992]223    DO k = 1, nlev
[2346]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
[1992]226    END DO
[2346]227    zlat(nbp_lat) = plat(l)*xpi/180.
[644]228
[2346]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)
[1992]238    END DO
239  END IF
[644]240
241
[1992]242  ! MOMENT ANGULAIRE
[644]243
[2346]244  DO j = 1, nbp_lat-1
245    DO i = 1, nbp_lon
[1992]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
[2346]278  DO j = 1, nbp_lat-1
279    DO i = 1, nbp_lon
[1992]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
[2346]287  DO j = 2, nbp_lat-1
288    DO i = 1, nbp_lon
[1992]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
[2346]302  DO j = 2, nbp_lat-1
303    DO i = 1, nbp_lon
[1992]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
[2346]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)
[1992]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.