GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: dyn3d/sw_case_williamson91_6.F Lines: 0 42 0.0 %
Date: 2023-06-30 12:51:15 Branches: 0 68 0.0 %

Line Branch Exec Source
1
!
2
! $Id $
3
!
4
      SUBROUTINE sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
5
6
c=======================================================================
7
c
8
c   Author:    Thomas Dubos      original: 26/01/2010
9
c   -------
10
c
11
c   Subject:
12
c   ------
13
c   Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz
14
c
15
c   Method:
16
c   --------
17
c
18
c   Interface:
19
c   ----------
20
c
21
c      Input:
22
c      ------
23
c
24
c      Output:
25
c      -------
26
c
27
c=======================================================================
28
      USE comconst_mod, ONLY: cpp, omeg, rad
29
      USE comvert_mod, ONLY: ap, bp, preff
30
31
      IMPLICIT NONE
32
c-----------------------------------------------------------------------
33
c   Declararations:
34
c   ---------------
35
36
      include "dimensions.h"
37
      include "paramet.h"
38
      include "comgeom.h"
39
      include "iniprint.h"
40
41
c   Arguments:
42
c   ----------
43
44
c   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
51
c   Local:
52
c   ------
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
67
c NB : rad = 6371220 dans W91 (6371229 dans LMDZ)
68
c      omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ)
69
70
      IF(0==0) THEN
71
c Williamson et al. (1991) : onde de Rossby-Haurwitz
72
         teta = preff/rho/cpp
73
c 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)
91
c 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
106
c 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
120
c         ucov=0
121
c         vcov=0
122
      ELSE
123
c 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
131
c     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
141
c-----------------------------------------------------------------------