/[lmdze]/trunk/Sources/phylmd/yamada.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/yamada.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 81 - (hide annotations)
Wed Mar 5 14:38:41 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/phylmd/yamada.f90
File size: 4600 byte(s) 
 Converted to free source form files which were still in fixed source
form. The conversion was done using the polish mode of the NAG Fortran
Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

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