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	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	guez	342	! LMDZ4/libf/phylmd/coefcdrag.F90, version 1.1.1.1, 2004 May 19th
12	guez	40
13	guez	274	! Objet : calcul des drag coefficients au sol pour le moment et
14	guez	291	! les flux de chaleurs sensible et latente et calcul de la
15	guez	342	! pression au niveau de r\'ef\'erence.
16	guez	40
17	guez	291	! Ionela MUSAT, July, 1st, 2002
18	guez	274
19	guez	292	! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982. A short
20	guez	291	! history of the operational PBL parametrization at
21	guez	292	! ECMWF. Workshop on boundary layer parametrization, November
22	guez	291	! 1981, ECMWF, Reading, England. Page: 19. Equations in Table 1.
23
24	guez	292	! 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	guez	342	! Library:
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	343	! norm of the wind in the first model layer, in m s-1
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	343	REAL, intent(out):: cdragm(:) ! (knon) drag coefficient for momentum
54	guez	221
55	guez	291	REAL, intent(out):: cdragh(:) ! (knon)
56	guez	343	! drag coefficient for latent and sensible heat fluxes
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	343	REAL, PARAMETER:: ckap = 0.4, cb = 5., cc = 5., cd = 5.
63			REAL, PARAMETER:: cepdu2 = 0.1**2 ! in m2 s-2
64	guez	291	real, parameter:: f_ri_cd_min = 0.1
65	guez	274	INTEGER i, knon
66	guez	343	REAL du2 ! in m2 s-2
67			real tsolv, tvd, zscf, zucf
68			real cdn ! drag coefficient for neutral conditions
69	guez	40
70	guez	343	REAL ri
71	guez	291	! nombre de Richardson entre la surface et le niveau de reference
72			! zgeop / RG
73	guez	40
74	guez	274	!-------------------------------------------------------------------------
75	guez	40
76	guez	271	knon = assert_eq([size(speed), size(t), size(q), size(zgeop), size(ts), &
77	guez	292	size(qsurf), size(rugos), size(cdragm), size(cdragh)], "cdrag knon")
78	guez	291
79	guez	248	DO i = 1, knon
80	guez	343	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	guez	291	* (1. + RETV * q(i))
84	guez	343	ri = zgeop(i) * (tvd - tsolv) / (du2 * tvd)
85	guez	341	cdn = (ckap / log(1. + zgeop(i) / (RG * rugos(i))))**2
86	guez	62
87	guez	343	IF (ri < 0.) THEN
88			! Situation instable
89	guez	341	zucf = 1. / (1. + 3. * cb * cc * cdn &
90	guez	343	* SQRT(ABS(ri) * (1. + zgeop(i) / (RG * rugos(i)))))
91			cdragm(i) = cdn * max((1. - 2. * cb * ri * zucf), f_ri_cd_min)
92	guez	274
93	guez	291	IF (nsrf == is_oce) then
94	guez	343	! Cf. Miller et al. (1992)
95	guez	341	cdragh(i) = f_cdrag_oce * cdn * (1. + ((0.0016 / (cdn &
96	guez	343	* SQRT(du2))) &
97			* ABS(tvd - tsolv)(1. / 3.))1.25)**(1. / 1.25)
98	guez	291	else
99	guez	341	cdragh(i) = f_cdrag_ter * cdn &
100	guez	343	* max((1. - 3. * cb * ri * zucf), f_ri_cd_min)
101	guez	291	end IF
102	guez	292	ELSE
103			! Situation stable. Pour \'eviter les incoh\'erences dans
104	guez	343	! les cas tr\`es stables, on limite ri \`a 20. Cf Hess et
105	guez	292	! al. (1995).
106	guez	343	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	guez	341	cdragh(i) = merge(f_cdrag_oce, f_cdrag_ter, nsrf == is_oce) * cdn &
110	guez	343	* max(1. / (1. + 3. * CB * ri * zscf), f_ri_cd_min)
111	guez	62	ENDIF
112			END DO
113
114	guez	274	if (present(pref)) &
115	guez	343	pref = psol * exp(- zgeop / (RD * t * (1. + RETV * max(q, 0.))))
116	guez	274
117	guez	275	END SUBROUTINE cdrag
118	guez	62
119	guez	275	end module cdrag_m