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.

    ! Miller, M. J., A. C. M. Beljaars, T. N. Palmer, 1992. The
    ! sensitivity of the ECMWF model to the parameterization of
    ! evaporation from the tropical oceans. J. Climate, 5:418-434.

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