source: LMDZ6/branches/Amaury_dev/libf/dyn3d/sw_case_williamson91_6.F90 @ 5218

Last change on this file since 5218 was 5159, checked in by abarral, 5 months ago

Put dimensions.h and paramet.h into modules

  • 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
File size: 4.2 KB
Line 
1! $Id $
2
3SUBROUTINE sw_case_williamson91_6(vcov, ucov, teta, masse, ps)
4
5  !=======================================================================
6
7  !   Author:    Thomas Dubos      original: 26/01/2010
8  !   -------
9
10  !   Subject:
11  !   ------
12  !   Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz
13
14  !   Method:
15  !   --------
16
17  !   Interface:
18  !   ----------
19
20  !  Input:
21  !  ------
22
23  !  Output:
24  !  -------
25
26  !=======================================================================
27  USE comconst_mod, ONLY: cpp, omeg, rad
28  USE comvert_mod, ONLY: ap, bp, preff
29  USE lmdz_iniprint, ONLY: lunout, prt_level
30  USE lmdz_comgeom
31
32  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
33  USE lmdz_paramet
34  IMPLICIT NONE
35  !-----------------------------------------------------------------------
36  !   Declararations:
37  !   ---------------
38
39
40
41
42  !   Arguments:
43  !   ----------
44
45  !   variables dynamiques
46  REAL :: vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
47  REAL :: teta(ip1jmp1, llm)                 ! temperature potentielle
48  REAL :: ps(ip1jmp1)                       ! pression  au sol
49  REAL :: masse(ip1jmp1, llm)                ! masse d'air
50  REAL :: phis(ip1jmp1)                     ! geopotentiel au sol
51
52  !   Local:
53  !   ------
54
55  REAL :: p (ip1jmp1, llmp1)               ! pression aux interfac.des couches
56  REAL :: pks(ip1jmp1)                      ! exner au  sol
57  REAL :: pk(ip1jmp1, llm)                   ! exner au milieu des couches
58  REAL :: pkf(ip1jmp1, llm)                  ! exner filt.au milieu des couches
59  REAL :: alpha(ip1jmp1, llm), beta(ip1jmp1, llm)
60
61  REAL :: sinth, costh, costh2, Ath, Bth, Cth, lon, dps
62  INTEGER :: i, j, ij
63
64  REAL, PARAMETER :: rho = 1 ! masse volumique de l'air (arbitraire)
65  REAL, PARAMETER :: K = 7.848e-6  ! K = \omega
66  REAL, PARAMETER :: gh0 = 9.80616 * 8e3
67  INTEGER, PARAMETER :: R0 = 4, R1 = R0 + 1, R2 = R0 + 2         ! mode 4
68  ! NB : rad = 6371220 dans W91 (6371229 dans LMDZ)
69  ! omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ)
70
71  IF(0==0) THEN
72    ! Williamson et al. (1991) : onde de Rossby-Haurwitz
73    teta = preff / rho / cpp
74    ! geopotentiel (pression de surface)
75    DO j = 1, jjp1
76      costh2 = cos(rlatu(j))**2
77      Ath = (R0 + 1) * (costh2**2) + (2 * R0 * R0 - R0 - 2) * costh2 - 2 * R0 * R0
78      Ath = .25 * (K**2) * (costh2**(R0 - 1)) * Ath
79      Ath = .5 * K * (2 * omeg + K) * costh2 + Ath
80      Bth = (R1 * R1 + 1) - R1 * R1 * costh2
81      Bth = 2 * (omeg + K) * K / (R1 * R2) * (costh2**(R0 / 2)) * Bth
82      Cth = R1 * costh2 - R2
83      Cth = .25 * K * K * (costh2**R0) * Cth
84      DO i = 1, iip1
85        ij = (j - 1) * iip1 + i
86        lon = rlonv(i)
87        dps = Ath + Bth * cos(R0 * lon) + Cth * cos(2 * R0 * lon)
88        ps(ij) = rho * (gh0 + (rad**2) * dps)
89      enddo
90    enddo
91    WRITE(lunout, *) 'W91 ps', MAXVAL(ps), MINVAL(ps)
92    ! vitesse zonale ucov
93    DO j = 1, jjp1
94      costh = cos(rlatu(j))
95      costh2 = costh**2
96      Ath = rad * K * costh
97      Bth = R0 * (1 - costh2) - costh2
98      Bth = rad * K * Bth * (costh**(R0 - 1))
99      DO i = 1, iip1
100        ij = (j - 1) * iip1 + i
101        lon = rlonu(i)
102        ucov(ij, 1) = (Ath + Bth * cos(R0 * lon))
103      enddo
104    enddo
105    WRITE(lunout, *) 'W91 u', MAXVAL(ucov(:, 1)), MINVAL(ucov(:, 1))
106    ucov(:, 1) = ucov(:, 1) * cu
107    ! vitesse meridienne vcov
108    DO j = 1, jjm
109      sinth = sin(rlatv(j))
110      costh = cos(rlatv(j))
111      Ath = -rad * K * R0 * sinth * (costh**(R0 - 1))
112      DO i = 1, iip1
113        ij = (j - 1) * iip1 + i
114        lon = rlonv(i)
115        vcov(ij, 1) = Ath * sin(R0 * lon)
116      enddo
117    enddo
118    WRITE(lunout, *) 'W91 v', MAXVAL(vcov(:, 1)), MINVAL(vcov(:, 1))
119    vcov(:, 1) = vcov(:, 1) * cv
120
121    ! ucov=0
122    ! vcov=0
123  ELSE
124    ! test non-tournant, onde se propageant en latitude
125    DO j = 1, jjp1
126      DO i = 1, iip1
127        ij = (j - 1) * iip1 + i
128        ps(ij) = 1e5 * (1 + .1 * exp(-100 * (1 + sin(rlatu(j)))**2))
129      enddo
130    enddo
131
132    ! rho = preff/(cpp*teta)
133    teta = .01 * preff / cpp   ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j
134    ucov = 0.
135    vcov = 0.
136  END IF
137
138  CALL pression (ip1jmp1, ap, bp, ps, p)
139  CALL massdair(p, masse)
140
141END SUBROUTINE sw_case_williamson91_6
142!-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.