libf/phylmd/coefcdrag.f90

!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/coefcdrag.F90,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
!
!
!
!
      SUBROUTINE coefcdrag (klon, knon, nsrf, zxli, &
                            speed, t, q, zgeop, psol, &
                            ts, qsurf, rugos, okri, ri1, &
                            cdram, cdrah, cdran, zri1, pref)
      use indicesol
      use YOMCST
      use yoethf
      IMPLICIT none
!-------------------------------------------------------------------------
! Objet : calcul des cdrags pour le moment (cdram) et les flux de chaleur 
!         sensible et latente (cdrah), du cdrag neutre (cdran), 
!         du nombre de Richardson entre la surface et le niveau de reference 
!         (zri1) et de la pression au niveau de reference (pref).    
!
! I. Musat, 01.07.2002
!-------------------------------------------------------------------------
!
! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
! knon----input-I- nombre de points pour un type de surface
! nsrf----input-I- indice pour le type de surface; voir indicesol.inc
! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
! speed---input-R- module du vent au 1er niveau du modele
! t-------input-R- temperature de l'air au 1er niveau du modele
! q-------input-R- humidite de l'air au 1er niveau du modele
! zgeop---input-R- geopotentiel au 1er niveau du modele
! psol----input-R- pression au sol 
! ts------input-R- temperature de l'air a la surface
! qsurf---input-R- humidite de l'air a la surface
! rugos---input-R- rugosite
! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce 
!                  et zref par rapport au Ri entre la sfce et la 1ere couche
! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche
!
! cdram--output-R- cdrag pour le moment
! cdrah--output-R- cdrag pour les flux de chaleur latente et sensible
! cdran--output-R- cdrag neutre
! zri1---output-R- nb. Richardson entre la surface et la couche zgeop/RG
! pref---output-R- pression au niveau zgeop/RG
!
      INTEGER, intent(in) :: klon, knon, nsrf
      LOGICAL, intent(in) :: zxli 
      REAL, dimension(klon), intent(in) :: speed, t, q, zgeop, psol
      REAL, dimension(klon), intent(in) :: ts, qsurf, rugos, ri1 
      LOGICAL, intent(in) :: okri    
!
      REAL, dimension(klon), intent(out) :: cdram, cdrah, cdran, zri1, pref
!-------------------------------------------------------------------------
!
! Quelques constantes :
      REAL, parameter :: RKAR=0.40, CB=5.0, CC=5.0, CD=5.0
!
! Variables locales :
      INTEGER :: i
      REAL, dimension(klon) :: zdu2, zdphi, ztsolv, ztvd
      REAL, dimension(klon) :: zscf, friv, frih, zucf, zcr
      REAL, dimension(klon) :: zcfm1, zcfh1
      REAL, dimension(klon) :: zcfm2, zcfh2
      REAL, dimension(klon) :: trm0, trm1
!-------------------------------------------------------------------------
      REAL :: fsta, fins, x
      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
      fins(x) = SQRT(1.0-18.0*x)
!-------------------------------------------------------------------------
!
      DO i = 1, knon
!
       zdphi(i) = zgeop(i)
       zdu2(i) = speed(i)**2
       pref(i) = exp(log(psol(i)) - zdphi(i)/(RD*t(i)* &
                 (1.+ RETV * max(q(i),0.0))))
       ztsolv(i) = ts(i)
       ztvd(i) = t(i) * (psol(i)/pref(i))**RKAPPA
       trm0(i) = 1. + RETV * max(qsurf(i),0.0)
       trm1(i) = 1. + RETV * max(q(i),0.0)
       ztsolv(i) = ztsolv(i) * trm0(i)
       ztvd(i) = ztvd(i) * trm1(i)
       zri1(i) = zdphi(i)*(ztvd(i)-ztsolv(i))/(zdu2(i)*ztvd(i))
!
! on teste zri1 par rapport au Richardson de la 1ere couche ri1 
!
!IM +++
       IF(1.EQ.0) THEN
       IF (okri) THEN
         IF (ri1(i).GE.0.0.AND.zri1(i).LT.0.0) THEN
           zri1(i) = ri1(i)
         ELSE IF(ri1(i).LT.0.0.AND.zri1(i).GE.0.0) THEN
           zri1(i) = ri1(i)
         ENDIF 
       ENDIF
       ENDIF
!IM ---
! 
       cdran(i) = (RKAR/log(1.+zdphi(i)/(RG*rugos(i))))**2

       IF (zri1(i) .ge. 0.) THEN 
!
! situation stable : pour eviter les inconsistances dans les cas 
! tres stables on limite zri1 a 20. cf Hess et al. (1995)
!
         zri1(i) = min(20.,zri1(i))
!
         IF (.NOT.zxli) THEN
           zscf(i) = SQRT(1.+CD*ABS(zri1(i)))
           friv(i) = max(1. / (1.+2.*CB*zri1(i)/ zscf(i)), 0.1)
           zcfm1(i) = cdran(i) * friv(i)
           frih(i) = max(1./ (1.+3.*CB*zri1(i)*zscf(i)), 0.1 )
           zcfh1(i) = cdran(i) * frih(i)
           cdram(i) = zcfm1(i)
           cdrah(i) = zcfh1(i)
         ELSE
           cdram(i) = cdran(i)* fsta(zri1(i))
           cdrah(i) = cdran(i)* fsta(zri1(i))
         ENDIF
!
       ELSE
! 
! situation instable
!
         IF (.NOT.zxli) THEN
           zucf(i) = 1./(1.+3.0*CB*CC*cdran(i)*SQRT(ABS(zri1(i)) &
                 *(1.0+zdphi(i)/(RG*rugos(i)))))
           zcfm2(i) = cdran(i)*max((1.-2.0*CB*zri1(i)*zucf(i)),0.1)
           zcfh2(i) = cdran(i)*max((1.-3.0*CB*zri1(i)*zucf(i)),0.1)
           cdram(i) = zcfm2(i)
           cdrah(i) = zcfh2(i)
         ELSE
           cdram(i) = cdran(i)* fins(zri1(i))
           cdrah(i) = cdran(i)* fins(zri1(i))
         ENDIF
!
! cdrah sur l'ocean cf. Miller et al. (1992)
!
         zcr(i) = (0.0016/(cdran(i)*SQRT(zdu2(i))))*ABS(ztvd(i)-ztsolv(i)) &
               **(1./3.)
         IF (nsrf.EQ.is_oce) cdrah(i) = cdran(i)*(1.0+zcr(i)**1.25) &
                  **(1./1.25)
       ENDIF
!
      END DO
      RETURN 
      END SUBROUTINE coefcdrag
1	guez	3	!
2			! $Header: /home/cvsroot/LMDZ4/libf/phylmd/coefcdrag.F90,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
3			!
4			!
5			!
6			!
7			SUBROUTINE coefcdrag (klon, knon, nsrf, zxli, &
8			speed, t, q, zgeop, psol, &
9			ts, qsurf, rugos, okri, ri1, &
10			cdram, cdrah, cdran, zri1, pref)
11			use indicesol
12			use YOMCST
13			use yoethf
14			IMPLICIT none
15			!-------------------------------------------------------------------------
16			! Objet : calcul des cdrags pour le moment (cdram) et les flux de chaleur
17			! sensible et latente (cdrah), du cdrag neutre (cdran),
18			! du nombre de Richardson entre la surface et le niveau de reference
19			! (zri1) et de la pression au niveau de reference (pref).
20			!
21			! I. Musat, 01.07.2002
22			!-------------------------------------------------------------------------
23			!
24			! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
25			! knon----input-I- nombre de points pour un type de surface
26			! nsrf----input-I- indice pour le type de surface; voir indicesol.inc
27			! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
28			! speed---input-R- module du vent au 1er niveau du modele
29			! t-------input-R- temperature de l'air au 1er niveau du modele
30			! q-------input-R- humidite de l'air au 1er niveau du modele
31			! zgeop---input-R- geopotentiel au 1er niveau du modele
32			! psol----input-R- pression au sol
33			! ts------input-R- temperature de l'air a la surface
34			! qsurf---input-R- humidite de l'air a la surface
35			! rugos---input-R- rugosite
36			! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce
37			! et zref par rapport au Ri entre la sfce et la 1ere couche
38			! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche
39			!
40			! cdram--output-R- cdrag pour le moment
41			! cdrah--output-R- cdrag pour les flux de chaleur latente et sensible
42			! cdran--output-R- cdrag neutre
43			! zri1---output-R- nb. Richardson entre la surface et la couche zgeop/RG
44			! pref---output-R- pression au niveau zgeop/RG
45			!
46			INTEGER, intent(in) :: klon, knon, nsrf
47			LOGICAL, intent(in) :: zxli
48			REAL, dimension(klon), intent(in) :: speed, t, q, zgeop, psol
49			REAL, dimension(klon), intent(in) :: ts, qsurf, rugos, ri1
50			LOGICAL, intent(in) :: okri
51			!
52			REAL, dimension(klon), intent(out) :: cdram, cdrah, cdran, zri1, pref
53			!-------------------------------------------------------------------------
54			!
55			! Quelques constantes :
56			REAL, parameter :: RKAR=0.40, CB=5.0, CC=5.0, CD=5.0
57			!
58			! Variables locales :
59			INTEGER :: i
60			REAL, dimension(klon) :: zdu2, zdphi, ztsolv, ztvd
61			REAL, dimension(klon) :: zscf, friv, frih, zucf, zcr
62			REAL, dimension(klon) :: zcfm1, zcfh1
63			REAL, dimension(klon) :: zcfm2, zcfh2
64			REAL, dimension(klon) :: trm0, trm1
65			!-------------------------------------------------------------------------
66			REAL :: fsta, fins, x
67			fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
68			fins(x) = SQRT(1.0-18.0*x)
69			!-------------------------------------------------------------------------
70			!
71			DO i = 1, knon
72			!
73			zdphi(i) = zgeop(i)
74			zdu2(i) = speed(i)**2
75			pref(i) = exp(log(psol(i)) - zdphi(i)/(RDt(i) &
76			(1.+ RETV * max(q(i),0.0))))
77			ztsolv(i) = ts(i)
78			ztvd(i) = t(i) * (psol(i)/pref(i))**RKAPPA
79			trm0(i) = 1. + RETV * max(qsurf(i),0.0)
80			trm1(i) = 1. + RETV * max(q(i),0.0)
81			ztsolv(i) = ztsolv(i) * trm0(i)
82			ztvd(i) = ztvd(i) * trm1(i)
83			zri1(i) = zdphi(i)(ztvd(i)-ztsolv(i))/(zdu2(i)ztvd(i))
84			!
85			! on teste zri1 par rapport au Richardson de la 1ere couche ri1
86			!
87			!IM +++
88			IF(1.EQ.0) THEN
89			IF (okri) THEN
90			IF (ri1(i).GE.0.0.AND.zri1(i).LT.0.0) THEN
91			zri1(i) = ri1(i)
92			ELSE IF(ri1(i).LT.0.0.AND.zri1(i).GE.0.0) THEN
93			zri1(i) = ri1(i)
94			ENDIF
95			ENDIF
96			ENDIF
97			!IM ---
98			!
99			cdran(i) = (RKAR/log(1.+zdphi(i)/(RGrugos(i))))*2
100
101			IF (zri1(i) .ge. 0.) THEN
102			!
103			! situation stable : pour eviter les inconsistances dans les cas
104			! tres stables on limite zri1 a 20. cf Hess et al. (1995)
105			!
106			zri1(i) = min(20.,zri1(i))
107			!
108			IF (.NOT.zxli) THEN
109			zscf(i) = SQRT(1.+CD*ABS(zri1(i)))
110			friv(i) = max(1. / (1.+2.CBzri1(i)/ zscf(i)), 0.1)
111			zcfm1(i) = cdran(i) * friv(i)
112			frih(i) = max(1./ (1.+3.CBzri1(i)*zscf(i)), 0.1 )
113			zcfh1(i) = cdran(i) * frih(i)
114			cdram(i) = zcfm1(i)
115			cdrah(i) = zcfh1(i)
116			ELSE
117			cdram(i) = cdran(i)* fsta(zri1(i))
118			cdrah(i) = cdran(i)* fsta(zri1(i))
119			ENDIF
120			!
121			ELSE
122			!
123			! situation instable
124			!
125			IF (.NOT.zxli) THEN
126			zucf(i) = 1./(1.+3.0CBCCcdran(i)SQRT(ABS(zri1(i)) &
127			(1.0+zdphi(i)/(RGrugos(i)))))
128			zcfm2(i) = cdran(i)max((1.-2.0CBzri1(i)zucf(i)),0.1)
129			zcfh2(i) = cdran(i)max((1.-3.0CBzri1(i)zucf(i)),0.1)
130			cdram(i) = zcfm2(i)
131			cdrah(i) = zcfh2(i)
132			ELSE
133			cdram(i) = cdran(i)* fins(zri1(i))
134			cdrah(i) = cdran(i)* fins(zri1(i))
135			ENDIF
136			!
137			! cdrah sur l'ocean cf. Miller et al. (1992)
138			!
139			zcr(i) = (0.0016/(cdran(i)SQRT(zdu2(i))))ABS(ztvd(i)-ztsolv(i)) &
140			**(1./3.)
141			IF (nsrf.EQ.is_oce) cdrah(i) = cdran(i)(1.0+zcr(i)*1.25) &
142			**(1./1.25)
143			ENDIF
144			!
145			END DO
146			RETURN
147			END SUBROUTINE coefcdrag