source: trunk/libf/dyn3dpar/sw_case_williamson91_6.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5

File size: 4.3 KB
RevLine 
[1]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      IMPLICIT NONE
29c-----------------------------------------------------------------------
30c   Declararations:
31c   ---------------
32
33#include "dimensions.h"
34#include "paramet.h"
35#include "comvert.h"
36#include "comconst.h"
37#include "comgeom.h"
38#include "iniprint.h"
39
40c   Arguments:
41c   ----------
42
43c   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
50c   Local:
51c   ------
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
66c NB : rad = 6371220 dans W91 (6371229 dans LMDZ)
67c      omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ)
68 
69      IF(0==0) THEN
70c Williamson et al. (1991) : onde de Rossby-Haurwitz
71         teta = preff/rho/cpp
72c 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)
90c 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
105c 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         
119c         ucov=0
120c         vcov=0
121      ELSE
122c 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         
130c     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
139      END
140c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.