phylmd/Interface_surf/cdrag.f90

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 May 19th

    ! 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 r\'ef\'erence.

    ! 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.

    ! Library:
    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 in the first model layer, in m s-1

    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 for momentum

    REAL, intent(out):: cdragh(:) ! (knon)
    ! drag coefficient for latent and sensible heat fluxes

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

    ! Local:

    REAL, PARAMETER:: ckap = 0.4, cb = 5., cc = 5., cd = 5.
    REAL, PARAMETER:: cepdu2 = 0.1**2 ! in m2 s-2
    real, parameter:: f_ri_cd_min = 0.1
    INTEGER i, knon
    REAL du2 ! in m2 s-2
    real tsolv, tvd, zscf, zucf
    real cdn ! drag coefficient for neutral conditions

    REAL ri
    ! 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
       du2 = max(cepdu2, speed(i)**2)
       tsolv = ts(i) * (1. + RETV * max(qsurf(i), 0.))
       tvd = (t(i) + zgeop(i) / RCPD / (1. + RVTMP2 * q(i))) &
            * (1. + RETV * q(i))
       ri = zgeop(i) * (tvd - tsolv) / (du2 * tvd)
       cdn = (ckap / log(1. + zgeop(i) / (RG * rugos(i))))**2

       IF (ri < 0.) THEN
          ! Situation instable
          zucf = 1. / (1. + 3. * cb * cc * cdn &
               * SQRT(ABS(ri) * (1. + zgeop(i) / (RG * rugos(i)))))
          cdragm(i) = cdn * max((1. - 2. * cb * ri * zucf), f_ri_cd_min)

          IF (nsrf == is_oce) then
             ! Cf. Miller et al. (1992)
             cdragh(i) = f_cdrag_oce * cdn * (1. + ((0.0016 / (cdn &
                  * SQRT(du2))) &
                  * ABS(tvd - tsolv)**(1. / 3.))**1.25)**(1. / 1.25)
          else
             cdragh(i) = f_cdrag_ter * cdn &
                  * max((1. - 3. * cb * ri * zucf), f_ri_cd_min)
          end IF
       ELSE
          ! Situation stable. Pour \'eviter les incoh\'erences dans
          ! les cas tr\`es stables, on limite ri \`a 20. Cf Hess et
          ! al. (1995).
          ri = min(20., ri)
          zscf = SQRT(1. + cd * ABS(ri))
          cdragm(i) = cdn * max(1. / (1. + 2. * CB * ri / zscf), f_ri_cd_min)
          cdragh(i) = merge(f_cdrag_oce, f_cdrag_ter, nsrf == is_oce) * cdn &
               * max(1. / (1. + 3. * CB * ri * zscf), f_ri_cd_min)
       ENDIF
    END DO

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

  END SUBROUTINE cdrag

end module cdrag_m
1	module cdrag_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE cdrag(nsrf, speed, t, q, zgeop, psol, ts, qsurf, rugos, cdragm, &
8	cdragh, pref)
9
10	! From LMDZ4/libf/phylmd/clcdrag.F90 and
11	! LMDZ4/libf/phylmd/coefcdrag.F90, version 1.1.1.1, 2004 May 19th
12
13	! Objet : calcul des drag coefficients au sol pour le moment et
14	! les flux de chaleurs sensible et latente et calcul de la
15	! pression au niveau de r\'ef\'erence.
16
17	! Ionela MUSAT, July, 1st, 2002
18
19	! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982. A short
20	! history of the operational PBL parametrization at
21	! ECMWF. Workshop on boundary layer parametrization, November
22	! 1981, ECMWF, Reading, England. Page: 19. Equations in Table 1.
23
24	! Miller, M. J., A. C. M. Beljaars, T. N. Palmer, 1992. The
25	! sensitivity of the ECMWF model to the parameterization of
26	! evaporation from the tropical oceans. J. Climate, 5:418-434.
27
28	! Library:
29	use nr_util, only: assert_eq
30
31	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	USE yoethf_m, ONLY: rvtmp2
35
36	INTEGER, intent(in):: nsrf ! indice pour le type de surface
37
38	REAL, intent(in):: speed(:) ! (knon)
39	! norm of the wind in the first model layer, in m s-1
40
41	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
46	REAL, intent(in):: zgeop(:) ! (knon)
47	! g\'eopotentiel au 1er niveau du mod\`ele
48
49	REAL, intent(in) :: psol(:) ! (knon) pression au sol
50	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	REAL, intent(in):: rugos(:) ! (knon) rugosit\'e
53	REAL, intent(out):: cdragm(:) ! (knon) drag coefficient for momentum
54
55	REAL, intent(out):: cdragh(:) ! (knon)
56	! drag coefficient for latent and sensible heat fluxes
57
58	REAL, intent(out), optional:: pref(:) ! (knon) pression au niveau zgeop / RG
59
60	! Local:
61
62	REAL, PARAMETER:: ckap = 0.4, cb = 5., cc = 5., cd = 5.
63	REAL, PARAMETER:: cepdu2 = 0.1**2 ! in m2 s-2
64	real, parameter:: f_ri_cd_min = 0.1
65	INTEGER i, knon
66	REAL du2 ! in m2 s-2
67	real tsolv, tvd, zscf, zucf
68	real cdn ! drag coefficient for neutral conditions
69
70	REAL ri
71	! nombre de Richardson entre la surface et le niveau de reference
72	! zgeop / RG
73
74	!-------------------------------------------------------------------------
75
76	knon = assert_eq([size(speed), size(t), size(q), size(zgeop), size(ts), &
77	size(qsurf), size(rugos), size(cdragm), size(cdragh)], "cdrag knon")
78
79	DO i = 1, knon
80	du2 = max(cepdu2, speed(i)**2)
81	tsolv = ts(i) * (1. + RETV * max(qsurf(i), 0.))
82	tvd = (t(i) + zgeop(i) / RCPD / (1. + RVTMP2 * q(i))) &
83	* (1. + RETV * q(i))
84	ri = zgeop(i) * (tvd - tsolv) / (du2 * tvd)
85	cdn = (ckap / log(1. + zgeop(i) / (RG * rugos(i))))**2
86
87	IF (ri < 0.) THEN
88	! Situation instable
89	zucf = 1. / (1. + 3. * cb * cc * cdn &
90	* SQRT(ABS(ri) * (1. + zgeop(i) / (RG * rugos(i)))))
91	cdragm(i) = cdn * max((1. - 2. * cb * ri * zucf), f_ri_cd_min)
92
93	IF (nsrf == is_oce) then
94	! Cf. Miller et al. (1992)
95	cdragh(i) = f_cdrag_oce * cdn * (1. + ((0.0016 / (cdn &
96	* SQRT(du2))) &
97	* ABS(tvd - tsolv)(1. / 3.))1.25)**(1. / 1.25)
98	else
99	cdragh(i) = f_cdrag_ter * cdn &
100	* max((1. - 3. * cb * ri * zucf), f_ri_cd_min)
101	end IF
102	ELSE
103	! Situation stable. Pour \'eviter les incoh\'erences dans
104	! les cas tr\`es stables, on limite ri \`a 20. Cf Hess et
105	! al. (1995).
106	ri = min(20., ri)
107	zscf = SQRT(1. + cd * ABS(ri))
108	cdragm(i) = cdn * max(1. / (1. + 2. * CB * ri / zscf), f_ri_cd_min)
109	cdragh(i) = merge(f_cdrag_oce, f_cdrag_ter, nsrf == is_oce) * cdn &
110	* max(1. / (1. + 3. * CB * ri * zscf), f_ri_cd_min)
111	ENDIF
112	END DO
113
114	if (present(pref)) &
115	pref = psol * exp(- zgeop / (RD * t * (1. + RETV * max(q, 0.))))
116
117	END SUBROUTINE cdrag
118
119	end module cdrag_m