source: LMDZ6/branches/IPSLCM6.0.14/libf/dyn3d/sw_case_williamson91_6.F @ 5444

Last change on this file since 5444 was 2600, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90
EM

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