Sources/phylmd/yamada.f

module yamada_m

  IMPLICIT NONE

contains

  SUBROUTINE yamada(ngrid, g, zlev, zlay, u, v, teta, q2, km, kn)

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