source: LMDZ6/trunk/libf/phylmd/surf_param_mod.F90 @ 5926

Last change on this file since 5926 was 5627, checked in by amaison, 10 months ago

Representation of heterogeneous continental subsurfaces with parameter or flux aggregation in the simplified surface model (bucket) for 1D case studies.
File size: 6.9 KB
Line 
1MODULE surf_param_mod
2!
3  IMPLICIT NONE
4!
5CONTAINS
6!
7!-------------------------------------------------------------------------------
8!
9FUNCTION eff_surf_param(klon, nbtersurf, x, frac, hatype, Zref) RESULT(eff_param)
10!
11!-------------------------------------------------------------------------------
12!
13! Arguments:
14  INTEGER, INTENT(IN)                          :: klon       ! grid point
15  INTEGER, INTENT(IN)                          :: nbtersurf  ! number of land hetero. subsurfaces
16  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: x          ! variable or parameter to integrate
17  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: frac       ! fraction of each land hetero. subsurface
18  CHARACTER(LEN=3), INTENT(IN)                 :: hatype     ! method to integrate the parameter
19  REAL, OPTIONAL, DIMENSION(klon), INTENT(IN)  :: Zref       ! reference height for CDN averaging (m)
20  REAL, DIMENSION(klon)                        :: eff_param  ! effective parameter
21!-------------------------------------------------------------------------------
22! Local variables:
23  INTEGER ik, is
24  REAL :: zrefd = 10.  ! default reference height for CDN averaging (m)
25  REAL :: Cdref        ! reference height for CDN averaging (m)
26!-------------------------------------------------------------------------------
27!
28  eff_param(:) = 0.
29  DO ik = 1, klon
30    DO is = 1, nbtersurf
31      !
32      ! arithmetic averaging
33      IF (hatype == 'ARI') THEN
34        eff_param(ik) = eff_param(ik) + frac(ik,is) * x(ik,is)
35      !
36      ! inverse averaging
37      ELSEIF (hatype == 'INV') THEN
38        IF (x(ik,is) .NE. 0.) THEN
39          eff_param(ik) = eff_param(ik) + frac(ik,is) * 1./x(ik,is)
40        ENDIF
41      !
42      ! inverse of square logarithm averaging
43      ELSEIF (hatype == 'CDN') THEN
44        IF (PRESENT(Zref)) THEN
45          Cdref = Zref(ik)
46        ELSE
47          Cdref = zrefd
48        ENDIF
49      !
50        IF (x(ik,is) .NE. 0.) THEN
51          eff_param(ik) = eff_param(ik) + frac(ik,is) * 1./(LOG(Cdref/x(ik,is)))**2
52        ENDIF
53      !
54      ELSE
55        PRINT*, 'eff_surf_param: invalid averaging type: ', hatype
56      ENDIF
57    ENDDO
58    !
59    IF (hatype == 'CDN') THEN
60      eff_param(ik) = Cdref * exp(-sqrt(1./eff_param(ik)))
61    ENDIF
62  !
63  ENDDO
64!
65END FUNCTION eff_surf_param
66!
67!
68!-------------------------------------------------------------------------------
69!
70FUNCTION average_surf_var(klon, nbtersurf, x, frac, hatype) RESULT(x_avg)
71!
72!-------------------------------------------------------------------------------
73!
74! Arguments:
75  INTEGER, INTENT(IN)                          :: klon       ! grid point
76  INTEGER, INTENT(IN)                          :: nbtersurf  ! number of land hetero. subsurfaces
77  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: x          ! variable or parameter to integrate
78  REAL, DIMENSION(klon, nbtersurf), INTENT(IN) :: frac       ! fraction of each land hetero. subsurface
79  CHARACTER(LEN=3), INTENT(IN)                 :: hatype     ! method to integrate the parameter
80  REAL, DIMENSION(klon)                        :: x_avg      ! average variable
81!
82! Local variables:
83  INTEGER ik, is
84!
85!-------------------------------------------------------------------------------
86!
87  x_avg(:) = 0.
88  DO ik = 1, klon
89    DO is = 1, nbtersurf
90      !
91      ! arithmetic averaging
92      IF (hatype == 'ARI') THEN
93        x_avg(ik) = x_avg(ik) + frac(ik,is) * x(ik,is)
94      ELSE
95        PRINT*, 'average_surf_var: invalid averaging type: ', hatype
96      ENDIF
97    ENDDO
98  ENDDO
99!
100END FUNCTION average_surf_var
101!
102!-------------------------------------------------------------------------------
103!
104FUNCTION interpol_tsoil(klon, nbtersurf, nsoilmx, nbtsoildepths, alpha, period, inertie, hcond, tsoil_depth, tsurf, tsoil_) RESULT(tsoil)
105!
106!-------------------------------------------------------------------------------
107!
108! Arguments:
109  INTEGER, INTENT(IN)                                         :: klon          ! grid point
110  INTEGER, INTENT(IN)                                         :: nbtersurf     ! number of land hetero. subsurfaces
111  INTEGER, INTENT(IN)                                         :: nsoilmx       ! number of soil layers in the model grid
112  INTEGER, INTENT(IN)                                         :: nbtsoildepths ! number of soil depths for soil temperature initialization
113  REAL, INTENT(IN)                                            :: alpha         ! parameter for soil discretization
114  REAL, INTENT(IN)                                            :: period        ! parameter for soil discretization
115  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: inertie       ! soil thermal inertia
116  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: hcond         ! soil heat conductivity
117  REAL, DIMENSION(klon, nbtsoildepths, nbtersurf), INTENT(IN) :: tsoil_depth   ! soil depth at which temperature is given (m)
118  REAL, DIMENSION(klon, nbtersurf), INTENT(IN)                :: tsurf         ! surface temperature
119  REAL, DIMENSION(klon, nbtsoildepths, nbtersurf), INTENT(IN) :: tsoil_        ! soil temperature given at tsoil_depths
120  REAL, DIMENSION(klon, nsoilmx, nbtersurf)                   :: tsoil         ! soil temperature interpolated in the model grid
121!
122! Local variables:
123  INTEGER ik, is, iq, it
124  REAL pi, slope, inter
125  REAL, DIMENSION(klon, nbtersurf)          :: z1, hcap  ! first layer depth and soil heat capacity
126  REAL, DIMENSION(klon, nsoilmx, nbtersurf) :: z         ! depth of the middle of the soil layer in the model grid (m)
127!
128!-------------------------------------------------------------------------------
129!
130  pi = ACOS(-1.)
131  !
132  DO ik = 1, klon
133    DO is = 1, nbtersurf
134      !
135      hcap(ik,is) = inertie(ik,is)*inertie(ik,is)/hcond(ik,is)
136      z1(ik,is) = SQRT(period*hcond(ik,is)/(pi*hcap(ik,is)))
137      !
138      DO iq = 1, nsoilmx
139        ! compute depth of middle soil layer (in m)
140        z(ik,iq,is) = z1(ik,is) * (alpha**(iq-0.5) - 1.) / (alpha - 1.)
141        ! if z is between the surface and first tsoil_depth
142        IF ((z(ik,iq,is) .GT. 0.) .AND. (z(ik,iq,is) .LE. tsoil_depth(ik,1,is))) THEN
143          slope = (tsoil_(ik,1,is) - tsurf(ik,is)) / (tsoil_depth(ik,1,is) - 0.)
144          inter = tsurf(ik,is)
145          tsoil(ik,iq,is) = slope * z(ik,iq,is) + inter
146        ENDIF
147        ! other levels
148        DO it = 1, nbtsoildepths-1
149          IF ((z(ik,iq,is) .GT. tsoil_depth(ik,it,is)) .AND. (z(ik,iq,is) .LE. tsoil_depth(ik,it+1,is))) THEN
150            slope = (tsoil_(ik,it+1,is) - tsoil_(ik,it,is)) / (tsoil_depth(ik,it+1,is) - tsoil_depth(ik,it,is))
151            inter = tsoil_(ik,it,is) - slope * tsoil_depth(ik,it,is)
152            tsoil(ik,iq,is) = slope * z(ik,iq,is) + inter
153          ENDIF
154        ENDDO
155        ! for layers below
156        IF (z(ik,iq,is) .GT. tsoil_depth(ik,nbtsoildepths,is)) THEN
157          tsoil(ik,iq,is) = tsoil_(ik,nbtsoildepths,is)
158        ENDIF
159      ENDDO
160      !
161    ENDDO
162  ENDDO
163!
164END FUNCTION interpol_tsoil
165!
166!-------------------------------------------------------------------------------
167!
168END MODULE surf_param_mod
Note: See TracBrowser for help on using the repository browser.