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

Last change on this file since 5157 was 5136, checked in by abarral, 8 weeks ago

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