phylmd/Interface_surf/cdrag.f

module cdrag_m

  IMPLICIT NONE

contains

  SUBROUTINE cdrag(nsrf, speed, t, q, zgeop, psol, ts, qsurf, rugos, cdragm, &
       cdragh, pref)

    ! From LMDZ4/libf/phylmd/clcdrag.F90 and
    ! LMDZ4/libf/phylmd/coefcdrag.F90, version 1.1.1.1, 2004/05/19
    ! 12:53:07

    ! Objet : calcul des drag coefficients au sol pour le moment et
    ! les flux de chaleurs sensible et latente et calcul de la
    ! pression au niveau de reference.

    ! Ionela MUSAT, July, 1st, 2002

    ! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short
    ! history of the operational PBL parametrization at
    ! ECMWF'. Workshop on boundary layer parametrization, November
    ! 1981, ECMWF, Reading, England. Page: 19. Equations in Table 1.

    use nr_util, only: assert_eq

    use clesphys, only: f_cdrag_oce, f_cdrag_ter
    use indicesol, only: is_oce
    use SUPHEC_M, only: rcpd, rd, retv, rg
    USE yoethf_m, ONLY: rvtmp2

    INTEGER, intent(in):: nsrf ! indice pour le type de surface

    REAL, intent(in):: speed(:) ! (knon)
    ! norm of the wind at the first model level

    REAL, intent(in):: t(:) ! (knon)
    ! temperature de l'air au 1er niveau du modele

    REAL, intent(in):: q(:) ! (knon) ! humidite de l'air au 1er niveau du modele

    REAL, intent(in):: zgeop(:) ! (knon)
    ! g\'eopotentiel au 1er niveau du mod\`ele

    REAL, intent(in) :: psol(:) ! (knon) pression au sol
    REAL, intent(in):: ts(:) ! (knon) temperature de l'air a la surface
    REAL, intent(in):: qsurf(:) ! (knon) humidite de l'air a la surface
    REAL, intent(in):: rugos(:) ! (knon) rugosit\'e
    REAL, intent(out):: cdragm(:) ! (knon) drag coefficient pour le moment 

    REAL, intent(out):: cdragh(:) ! (knon)
    ! drag coefficient pour les flux de chaleur latente et sensible

    REAL, intent(out), optional:: pref(:) ! (knon) pression au niveau zgeop / RG

    ! Local:

    REAL, PARAMETER:: ckap = 0.40, cb = 5., cc = 5., cd = 5., cepdu2 = 0.1**2
    real, parameter:: f_ri_cd_min = 0.1
    INTEGER i, knon
    REAL zdu2, ztsolv, ztvd, zscf, zucf
    real zcdn ! drag coefficient neutre

    REAL zri
    ! nombre de Richardson entre la surface et le niveau de reference
    ! zgeop / RG

    !-------------------------------------------------------------------------

    knon = assert_eq([size(speed), size(t), size(q), size(zgeop), size(ts), &
         size(qsurf), size(rugos), size(cdragm), size(cdragh)], &
         "cdrag knon")

    DO i = 1, knon
       zdu2 = max(cepdu2, speed(i)**2)
       ztsolv = ts(i) * (1. + RETV * max(qsurf(i), 0.))
       ztvd = (t(i) + zgeop(i) / RCPD / (1. + RVTMP2 * q(i))) &
            * (1. + RETV * q(i))
       zri = zgeop(i) * (ztvd - ztsolv) / (zdu2 * ztvd)
       zcdn = (ckap / log(1. + zgeop(i) / (RG * rugos(i))))**2

       IF (zri > 0.) THEN
          ! Situation stable. Pour \'eviter les incoh\'erences dans
          ! les cas tr\`es stables, on limite zri \`a 20. Cf Hess et
          ! al. (1995).
          zri = min(20., zri)
          zscf = SQRT(1. + cd * ABS(zri))
          cdragm(i) = zcdn * max(1. / (1. + 2. * CB * zri / zscf), f_ri_cd_min)
          cdragh(i) = merge(f_cdrag_oce, f_cdrag_ter, nsrf == is_oce) * zcdn &
               * max(1. / (1. + 3. * CB * zri * zscf), f_ri_cd_min)
       ELSE
          ! situation instable
          zucf = 1. / (1. + 3. * cb * cc * zcdn &
               * SQRT(ABS(zri) * (1. + zgeop(i) / (RG * rugos(i)))))
          cdragm(i) = zcdn * max((1. - 2. * cb * zri * zucf), f_ri_cd_min)

          IF (nsrf == is_oce) then
             ! Cf. Miller et al. (1992).
             cdragh(i) = f_cdrag_oce * zcdn * (1. + ((0.0016 &
                  / (zcdn * SQRT(zdu2))) * ABS(ztvd - ztsolv)**(1. &
                  / 3.))**1.25)**(1. / 1.25)
          else
             cdragh(i) = f_cdrag_ter * zcdn &
                  * max((1. - 3. * cb * zri * zucf), f_ri_cd_min)
          end IF
       ENDIF
    END DO

    if (present(pref)) &
         pref = exp(log(psol) - zgeop / (RD * t * (1. + RETV * max(q, 0.))))

  END SUBROUTINE cdrag

end module cdrag_m
1	guez	275	module cdrag_m
2	guez	40
3	guez	62	IMPLICIT NONE
4	guez	40
5	guez	62	contains
6	guez	40
7	guez	291	SUBROUTINE cdrag(nsrf, speed, t, q, zgeop, psol, ts, qsurf, rugos, cdragm, &
8			cdragh, pref)
9	guez	40
10	guez	275	! From LMDZ4/libf/phylmd/clcdrag.F90 and
11			! LMDZ4/libf/phylmd/coefcdrag.F90, version 1.1.1.1, 2004/05/19
12			! 12:53:07
13	guez	40
14	guez	274	! Objet : calcul des drag coefficients au sol pour le moment et
15	guez	291	! les flux de chaleurs sensible et latente et calcul de la
16			! pression au niveau de reference.
17	guez	40
18	guez	291	! Ionela MUSAT, July, 1st, 2002
19	guez	274
20	guez	291	! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short
21			! history of the operational PBL parametrization at
22			! ECMWF'. Workshop on boundary layer parametrization, November
23			! 1981, ECMWF, Reading, England. Page: 19. Equations in Table 1.
24
25			use nr_util, only: assert_eq
26
27	guez	274	use clesphys, only: f_cdrag_oce, f_cdrag_ter
28			use indicesol, only: is_oce
29			use SUPHEC_M, only: rcpd, rd, retv, rg
30	guez	221	USE yoethf_m, ONLY: rvtmp2
31
32			INTEGER, intent(in):: nsrf ! indice pour le type de surface
33
34	guez	271	REAL, intent(in):: speed(:) ! (knon)
35	guez	274	! norm of the wind at the first model level
36	guez	40
37	guez	248	REAL, intent(in):: t(:) ! (knon)
38			! temperature de l'air au 1er niveau du modele
39
40			REAL, intent(in):: q(:) ! (knon) ! humidite de l'air au 1er niveau du modele
41	guez	274
42			REAL, intent(in):: zgeop(:) ! (knon)
43			! g\'eopotentiel au 1er niveau du mod\`ele
44	guez	291
45	guez	272	REAL, intent(in) :: psol(:) ! (knon) pression au sol
46	guez	221	REAL, intent(in):: ts(:) ! (knon) temperature de l'air a la surface
47			REAL, intent(in):: qsurf(:) ! (knon) humidite de l'air a la surface
48	guez	248	REAL, intent(in):: rugos(:) ! (knon) rugosit\'e
49	guez	291	REAL, intent(out):: cdragm(:) ! (knon) drag coefficient pour le moment
50	guez	221
51	guez	291	REAL, intent(out):: cdragh(:) ! (knon)
52	guez	271	! drag coefficient pour les flux de chaleur latente et sensible
53	guez	47
54	guez	291	REAL, intent(out), optional:: pref(:) ! (knon) pression au niveau zgeop / RG
55	guez	272
56	guez	221	! Local:
57	guez	291
58			REAL, PARAMETER:: ckap = 0.40, cb = 5., cc = 5., cd = 5., cepdu2 = 0.1**2
59			real, parameter:: f_ri_cd_min = 0.1
60	guez	274	INTEGER i, knon
61	guez	291	REAL zdu2, ztsolv, ztvd, zscf, zucf
62	guez	274	real zcdn ! drag coefficient neutre
63	guez	40
64	guez	274	REAL zri
65	guez	291	! nombre de Richardson entre la surface et le niveau de reference
66			! zgeop / RG
67	guez	40
68	guez	274	!-------------------------------------------------------------------------
69	guez	40
70	guez	271	knon = assert_eq([size(speed), size(t), size(q), size(zgeop), size(ts), &
71	guez	291	size(qsurf), size(rugos), size(cdragm), size(cdragh)], &
72	guez	275	"cdrag knon")
73	guez	291
74	guez	248	DO i = 1, knon
75	guez	274	zdu2 = max(cepdu2, speed(i)**2)
76			ztsolv = ts(i) * (1. + RETV * max(qsurf(i), 0.))
77	guez	291	ztvd = (t(i) + zgeop(i) / RCPD / (1. + RVTMP2 * q(i))) &
78			* (1. + RETV * q(i))
79			zri = zgeop(i) * (ztvd - ztsolv) / (zdu2 * ztvd)
80			zcdn = (ckap / log(1. + zgeop(i) / (RG * rugos(i))))**2
81	guez	62
82	guez	274	IF (zri > 0.) THEN
83	guez	291	! Situation stable. Pour \'eviter les incoh\'erences dans
84			! les cas tr\`es stables, on limite zri \`a 20. Cf Hess et
85			! al. (1995).
86	guez	274	zri = min(20., zri)
87	guez	291	zscf = SQRT(1. + cd * ABS(zri))
88			cdragm(i) = zcdn * max(1. / (1. + 2. * CB * zri / zscf), f_ri_cd_min)
89			cdragh(i) = merge(f_cdrag_oce, f_cdrag_ter, nsrf == is_oce) * zcdn &
90			* max(1. / (1. + 3. * CB * zri * zscf), f_ri_cd_min)
91	guez	274	ELSE
92	guez	62	! situation instable
93	guez	291	zucf = 1. / (1. + 3. * cb * cc * zcdn &
94			* SQRT(ABS(zri) * (1. + zgeop(i) / (RG * rugos(i)))))
95			cdragm(i) = zcdn * max((1. - 2. * cb * zri * zucf), f_ri_cd_min)
96	guez	274
97	guez	291	IF (nsrf == is_oce) then
98			! Cf. Miller et al. (1992).
99			cdragh(i) = f_cdrag_oce * zcdn * (1. + ((0.0016 &
100			/ (zcdn * SQRT(zdu2))) * ABS(ztvd - ztsolv)**(1. &
101			/ 3.))1.25)(1. / 1.25)
102			else
103			cdragh(i) = f_cdrag_ter * zcdn &
104			* max((1. - 3. * cb * zri * zucf), f_ri_cd_min)
105			end IF
106	guez	62	ENDIF
107			END DO
108
109	guez	274	if (present(pref)) &
110			pref = exp(log(psol) - zgeop / (RD * t * (1. + RETV * max(q, 0.))))
111
112	guez	275	END SUBROUTINE cdrag
113	guez	62
114	guez	275	end module cdrag_m