trunk/phylmd/yamada.f

module yamada_m

  IMPLICIT NONE

contains

  SUBROUTINE yamada(ngrid, g, rconst, plev, temp, zlev, zlay, u, v, teta, q2, &
       km, kn, ustar, l_mix)

    ! From LMDZ4/libf/phylmd/yamada.F,v 1.1 2004/06/22 11:45:36

    USE dimens_m
    USE dimphy
    ! .......................................................................
    ! .......................................................................

    ! g  : g
    ! zlev : altitude a chaque niveau (interface inferieure de la couche
    ! de meme indice)
    ! zlay : altitude au centre de chaque couche
    ! u,v : vitesse au centre de chaque couche
    ! (en entree : la valeur au debut du pas de temps)
    ! teta : temperature potentielle au centre de chaque couche
    ! (en entree : la valeur au debut du pas de temps)
    ! q2 : $q^2$ au bas de chaque couche
    ! (en entree : la valeur au debut du pas de temps)
    ! (en sortie : la valeur a la fin du pas de temps)
    ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
    ! couche)
    ! (en sortie : la valeur a la fin du pas de temps)
    ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
    ! (en sortie : la valeur a la fin du pas de temps)

    ! .......................................................................
    REAL, INTENT (IN) :: g
    REAL rconst
    REAL plev(klon, klev+1), temp(klon, klev)
    REAL ustar(klon), snstable
    REAL zlev(klon, klev+1)
    REAL zlay(klon, klev)
    REAL u(klon, klev)
    REAL v(klon, klev)
    REAL teta(klon, klev)
    REAL q2(klon, klev+1)
    REAL km(klon, klev+1)
    REAL kn(klon, klev+1)
    INTEGER l_mix, ngrid


    INTEGER nlay, nlev
    PARAMETER (nlay=klev)
    PARAMETER (nlev=klev+1)

    LOGICAL first
    SAVE first
    DATA first/.TRUE./


    INTEGER ig, k

    REAL ri, zrif, zalpha, zsm
    REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)

    REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
    REAL l(klon, klev+1), l0(klon)

    REAL sq(klon), sqz(klon), zz(klon, klev+1)
    INTEGER iter

    REAL ric, rifc, b1, kap
    SAVE ric, rifc, b1, kap
    DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.3/

    IF (0==1 .AND. first) THEN
       DO ig = 1, 1000
          ri = (ig-800.)/500.
          IF (ri<ric) THEN
             zrif = frif(ri)
          ELSE
             zrif = rifc
          END IF
          IF (zrif<0.16) THEN
             zalpha = falpha(zrif)
             zsm = fsm(zrif)
          ELSE
             zalpha = 1.12
             zsm = 0.085
          END IF
          PRINT *, ri, rif, zalpha, zsm
       END DO
       first = .FALSE.
    END IF

    ! Correction d'un bug sauvage a verifier.
    ! do k=2,nlev
    DO k = 2, nlay
       DO ig = 1, ngrid
          dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
          m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
               k-1))**2)/(dz(ig,k)*dz(ig,k))
          n2(ig, k) = g*2.*(teta(ig,k)-teta(ig,k-1))/(teta(ig,k-1)+teta(ig,k))/ &
               dz(ig, k)
          ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
          IF (ri<ric) THEN
             rif(ig, k) = frif(ri)
          ELSE
             rif(ig, k) = rifc
          END IF
          IF (rif(ig,k)<0.16) THEN
             alpha(ig, k) = falpha(rif(ig,k))
             sm(ig, k) = fsm(rif(ig,k))
          ELSE
             alpha(ig, k) = 1.12
             sm(ig, k) = 0.085
          END IF
          zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
       END DO
    END DO

    ! iterration pour determiner la longueur de melange

    DO ig = 1, ngrid
       l0(ig) = 100.
    END DO
    DO k = 2, klev - 1
       DO ig = 1, ngrid
          l(ig, k) = l0(ig)*kap*zlev(ig, k)/(kap*zlev(ig,k)+l0(ig))
       END DO
    END DO

    DO iter = 1, 10
       DO ig = 1, ngrid
          sq(ig) = 1.E-10
          sqz(ig) = 1.E-10
       END DO
       DO k = 2, klev - 1
          DO ig = 1, ngrid
             q2(ig, k) = l(ig, k)**2*zz(ig, k)
             l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
                  k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
             zq = sqrt(q2(ig,k))
             sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
             sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
          END DO
       END DO
       DO ig = 1, ngrid
          l0(ig) = 0.2*sqz(ig)/sq(ig)
       END DO
       ! (abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0

    END DO

    DO k = 2, klev
       DO ig = 1, ngrid
          l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
               k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
          q2(ig, k) = l(ig, k)**2*zz(ig, k)
          km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
          kn(ig, k) = km(ig, k)*alpha(ig, k)
       END DO
    END DO

  contains

    REAL function frif(ri)
      real ri
      frif = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    end function frif

    REAL function falpha(ri)
      real ri
      falpha = 1.318*(0.2231-ri)/(0.2341-ri)
    end function falpha

    REAL function fsm(ri)
      real ri
      fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
    end function fsm

  END SUBROUTINE yamada

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