GCC Code Coverage Report


Directory: ./
File: phys/yamada.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 44 0.0%
Branches: 0 28 0.0%

Line Branch Exec Source
1
2 ! $Header$
3
4 SUBROUTINE yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
5 cd, q2, km, kn, ustar, l_mix)
6 USE dimphy
7 IMPLICIT NONE
8
9 ! dt : pas de temps
10 ! g : g
11 ! zlev : altitude a chaque niveau (interface inferieure de la couche
12 ! de meme indice)
13 ! zlay : altitude au centre de chaque couche
14 ! u,v : vitesse au centre de chaque couche
15 ! (en entree : la valeur au debut du pas de temps)
16 ! teta : temperature potentielle au centre de chaque couche
17 ! (en entree : la valeur au debut du pas de temps)
18 ! cd : cdrag
19 ! (en entree : la valeur au debut du pas de temps)
20 ! q2 : $q^2$ au bas de chaque couche
21 ! (en entree : la valeur au debut du pas de temps)
22 ! (en sortie : la valeur a la fin du pas de temps)
23 ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
24 ! couche)
25 ! (en sortie : la valeur a la fin du pas de temps)
26 ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
27 ! (en sortie : la valeur a la fin du pas de temps)
28
29 ! .......................................................................
30 REAL dt, g, rconst
31 REAL plev(klon, klev+1), temp(klon, klev)
32 REAL ustar(klon), snstable
33 REAL zlev(klon, klev+1)
34 REAL zlay(klon, klev)
35 REAL u(klon, klev)
36 REAL v(klon, klev)
37 REAL teta(klon, klev)
38 REAL cd(klon)
39 REAL q2(klon, klev+1)
40 REAL km(klon, klev+1)
41 REAL kn(klon, klev+1)
42 INTEGER l_mix, ngrid
43
44
45 INTEGER nlay, nlev
46 ! ym PARAMETER (nlay=klev)
47 ! ym PARAMETER (nlev=klev+1)
48
49 LOGICAL first
50 SAVE first
51 DATA first/.TRUE./
52 !$OMP THREADPRIVATE(first)
53
54 INTEGER ig, k
55
56 REAL ri, zrif, zalpha, zsm
57 REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
58
59 REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
60 REAL l(klon, klev+1), l0(klon)
61
62 REAL sq(klon), sqz(klon), zz(klon, klev+1)
63 INTEGER iter
64
65 REAL ric, rifc, b1, kap
66 SAVE ric, rifc, b1, kap
67 DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.3/
68 !$OMP THREADPRIVATE(ric,rifc,b1,kap)
69
70 REAL frif, falpha, fsm
71
72 frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
73 falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
74 fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
75
76 nlay = klev
77 nlev = klev + 1
78
79 IF (0==1 .AND. first) THEN
80 DO ig = 1, 1000
81 ri = (ig-800.)/500.
82 IF (ri<ric) THEN
83 zrif = frif(ri)
84 ELSE
85 zrif = rifc
86 END IF
87 IF (zrif<0.16) THEN
88 zalpha = falpha(zrif)
89 zsm = fsm(zrif)
90 ELSE
91 zalpha = 1.12
92 zsm = 0.085
93 END IF
94 PRINT *, ri, rif, zalpha, zsm
95 END DO
96 first = .FALSE.
97 END IF
98
99 ! Correction d'un bug sauvage a verifier.
100 ! do k=2,nlev
101 DO k = 2, nlay
102 DO ig = 1, ngrid
103 dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
104 m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
105 k-1))**2)/(dz(ig,k)*dz(ig,k))
106 n2(ig, k) = g*2.*(teta(ig,k)-teta(ig,k-1))/(teta(ig,k-1)+teta(ig,k))/ &
107 dz(ig, k)
108 ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
109 IF (ri<ric) THEN
110 rif(ig, k) = frif(ri)
111 ELSE
112 rif(ig, k) = rifc
113 END IF
114 IF (rif(ig,k)<0.16) THEN
115 alpha(ig, k) = falpha(rif(ig,k))
116 sm(ig, k) = fsm(rif(ig,k))
117 ELSE
118 alpha(ig, k) = 1.12
119 sm(ig, k) = 0.085
120 END IF
121 zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
122 END DO
123 END DO
124
125 ! iterration pour determiner la longueur de melange
126
127 DO ig = 1, ngrid
128 l0(ig) = 100.
129 END DO
130 DO k = 2, klev - 1
131 DO ig = 1, ngrid
132 l(ig, k) = l0(ig)*kap*zlev(ig, k)/(kap*zlev(ig,k)+l0(ig))
133 END DO
134 END DO
135
136 DO iter = 1, 10
137 DO ig = 1, ngrid
138 sq(ig) = 1.E-10
139 sqz(ig) = 1.E-10
140 END DO
141 DO k = 2, klev - 1
142 DO ig = 1, ngrid
143 q2(ig, k) = l(ig, k)**2*zz(ig, k)
144 l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
145 k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
146 zq = sqrt(q2(ig,k))
147 sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
148 sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
149 END DO
150 END DO
151 DO ig = 1, ngrid
152 l0(ig) = 0.2*sqz(ig)/sq(ig)
153 END DO
154 ! (abd 3 5 2) print*,'ITER=',iter,' L0=',l0
155
156 END DO
157
158 DO k = 2, klev
159 DO ig = 1, ngrid
160 l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
161 k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
162 q2(ig, k) = l(ig, k)**2*zz(ig, k)
163 km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
164 kn(ig, k) = km(ig, k)*alpha(ig, k)
165 END DO
166 END DO
167
168 RETURN
169 END SUBROUTINE yamada
170