source: LMDZ6/trunk/libf/phylmd/surf_land_bucket_hetero_mod.F90 @ 5627

Last change on this file since 5627 was 5627, checked in by amaison, 6 weeks ago

Representation of heterogeneous continental subsurfaces with parameter or flux aggregation in the simplified surface model (bucket) for 1D case studies.
File size: 15.4 KB
Line 
1!
2MODULE surf_land_bucket_hetero_mod
3!
4! 2025/04 A. Maison (adapted from surf_land_bucket_mod)
5! Surface land bucket module with heterogeneous continental sub-surfaces
6! This module is used when no external land model is choosen and iflag_hetero_surf = 1 or 2.
7!
8  IMPLICIT NONE
9
10CONTAINS
11
12  SUBROUTINE surf_land_bucket_hetero(itime, jour, knon, knindex, debut, dtime,&
13       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
14       spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, pref, plev, &
15       u1, v1, gustiness, rugoro, swnet, lwnet, &
16       snow, qsol, agesno, tsoil, &
17       qsurf, z0m, z0h, alb1_new, alb2_new, evap, &
18       fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l, &
19       tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
20       cdragm_tersrf, cdragh_tersrf, &
21       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)
22
23    USE clesphys_mod_h
24    USE dimsoil_mod_h, ONLY: nsoilmx
25    USE YOMCST, ONLY: RD, RG, RCPD, RSIGMA
26    USE compbl_mod_h
27    USE dimpft_mod_h
28    USE limit_read_mod
29    USE surface_data
30    USE fonte_neige_mod
31    USE calcul_fluxs_mod
32    USE cpl_mod
33    USE dimphy
34    USE geometry_mod, ONLY: longitude,latitude
35    USE mod_grid_phy_lmdz
36    USE mod_phys_lmdz_para
37    USE indice_sol_mod
38    USE phys_state_var_mod, ONLY: frac_tersrf, ratio_z0m_z0h_tersrf, z0m_tersrf, &
39                                  albedo_tersrf, beta_tersrf, inertie_tersrf, &
40                                  hcond_tersrf
41    USE surf_param_mod, ONLY: eff_surf_param, average_surf_var
42    USE cdrag_mod
43
44!****************************************************************************************
45! Bucket calculations for surface.
46!****************************************************************************************
47!
48! Input variables 
49!****************************************************************************************
50    INTEGER, INTENT(IN)                     :: itime, jour, knon
51    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
52    LOGICAL, INTENT(IN)                     :: debut
53    REAL, INTENT(IN)                        :: dtime
54    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
55    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
56    REAL, DIMENSION(klon), INTENT(IN)       :: tq_cdrag
57    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
58    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
59    REAL, DIMENSION(klon), INTENT(IN)       :: petAcoef, peqAcoef
60    REAL, DIMENSION(klon), INTENT(IN)       :: petBcoef, peqBcoef
61    REAL, DIMENSION(klon), INTENT(IN)       :: pref
62    REAL, DIMENSION(klon), INTENT(IN)       :: u1, v1, gustiness
63    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
64    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
65    REAL, DIMENSION(klon, nbtersrf), INTENT(IN) :: tsurf_tersrf
66
67
68! In/Output variables
69!****************************************************************************************
70    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
71    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
72    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
73    REAL, DIMENSION(klon), INTENT(INOUT)          :: plev
74    REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf
75
76
77! Output variables
78!****************************************************************************************
79    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
80    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
81    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new, alb2_new
82    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
83    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
84    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
85    !
86    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: tsurf_new_tersrf
87    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: qsurf_tersrf
88    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragm_tersrf, cdragh_tersrf
89    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: swnet_tersrf, lwnet_tersrf
90    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxsens_tersrf, fluxlat_tersrf
91
92! Local variables
93!****************************************************************************************
94    REAL, DIMENSION(klon) :: soilcap, soilflux
95    REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
96    REAL, DIMENSION(klon) :: alb_neig, alb_lim, icesub
97    REAL, DIMENSION(klon) :: zfra
98    REAL, DIMENSION(klon) :: radsol
99    REAL, DIMENSION(klon) :: u0, v0, u1_lay, v1_lay
100    REAL, DIMENSION(klon) :: dummy_riverflow,dummy_coastalflow
101    INTEGER               :: i, j, k
102    !
103    REAL, DIMENSION(klon)           :: zlev, geop1, speed, pblh, ri_in, sst
104    REAL, DIMENSION(klon)           :: beta_eff, inertie_eff, conv_ratio_eff
105    REAL, DIMENSION(klon)           :: meansqT
106    REAL, DIMENSION(klon, nbtersrf) :: z0h_tersrf, emis_tersrf, conv_ratio_tersrf
107    REAL, DIMENSION(klon, nbtersrf) :: evap_tersrf
108    REAL, DIMENSION(klon, nbtersrf) :: dflux_s_tersrf, dflux_l_tersrf
109    REAL, DIMENSION(klon, nbtersrf) :: radsol_tersrf
110    REAL, DIMENSION(klon, nbtersrf) :: zri_tersrf
111    REAL, PARAMETER                 :: klon_1D = 1
112
113!
114!****************************************************************************************
115
116    ! *** Calculations common to the two flag values ***
117
118    ! average albedo
119    alb_lim = eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI')
120
121    ! suppose zero surface speed
122    u0(:)=0.0
123    v0(:)=0.0
124    u1_lay(:) = u1(:) - u0(:)
125    v1_lay(:) = v1(:) - v0(:)
126    speed(:) = (u1_lay(:)**2 + v1_lay(:)**2)**0.5
127    !
128    geop1(1:knon) = RD * temp_air(1:knon) / (0.5*(pref(1:knon)+p1lay(1:knon))) &
129            * (pref(1:knon)-p1lay(1:knon))
130    zlev(1:knon) = (plev(1:knon)*RD*temp_air(1:knon)/(pref(1:knon)*RG))/2.
131    !
132    ! compute roughness lengths
133    DO i=1, knon
134      DO j=1, nbtersrf
135        IF (ratio_z0m_z0h_tersrf(i,j) .NE. 0.) THEN
136          z0h_tersrf(i,j) = z0m_tersrf(i,j) / ratio_z0m_z0h_tersrf(i,j)
137        ELSE
138          z0h_tersrf(i,j) = 1.E-12
139        ENDIF
140      ENDDO
141    ENDDO
142
143    z0m = eff_surf_param(klon, nbtersrf, z0m_tersrf, frac_tersrf, 'CDN', zlev)
144    z0h = eff_surf_param(klon, nbtersrf, z0h_tersrf, frac_tersrf, 'CDN', zlev)
145
146    DO i=1, knon
147       z0m(i) = MAX(1.5e-05,SQRT(z0m(i)**2 + rugoro(i)**2))
148    END DO
149
150    ! compute the ratio to convert and print soil depths in meters (conv_ratio = (cond/cap)^0.5 and cap = I^2/cond)
151    DO j=1, nbtersrf
152       conv_ratio_tersrf(:,j) = hcond_tersrf(:,j)/inertie_tersrf(:,j)
153    ENDDO
154
155    !
156    ! *** Surface parameter aggregation ***
157    !
158    IF (iflag_hetero_surf == 1) THEN
159      !* Calcultaion of fluxes
160     
161      ! calculate total absorbed radiance at surface
162      radsol(:) = 0.0
163      radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
164
165      ! calculate constants (needeed for capsol and dif_grnd)
166      CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
167      IF (type_veget=='betaclim') THEN
168        CALL calbeta_clim(knon,jour,latitude(knindex(1:knon)),beta)
169      ENDIF
170
171      ! mean evapotranspiration coefficient
172      beta_eff = eff_surf_param(klon, nbtersrf, beta_tersrf, frac_tersrf, 'ARI')
173      beta = beta_eff
174
175      ! calculate temperature, heat capacity and conduction flux in soil
176      IF (soil_model) THEN
177         inertie_eff = eff_surf_param(klon, nbtersrf, inertie_tersrf, frac_tersrf, 'ARI')
178         conv_ratio_eff = eff_surf_param(klon, nbtersrf, conv_ratio_tersrf, frac_tersrf, 'ARI')
179         !
180         CALL soil_hetero(dtime, is_ter, knon, snow, tsurf, qsol,  &
181           & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux, &
182           & inertie_eff, conv_ratio_eff)
183         !
184         DO i=1, knon
185            cal(i) = RCPD / soilcap(i)
186            radsol(i) = radsol(i)  + soilflux(i)
187         END DO
188      ELSE
189         cal(:) = RCPD * capsol(:)
190         IF (klon_glo .EQ. 1) THEN
191           cal(:) = 0.
192         ENDIF
193      ENDIF
194
195      ! calculate fluxes
196      CALL calcul_fluxs(knon, is_ter, dtime, &
197           tsurf, p1lay, cal, beta, tq_cdrag, tq_cdrag, pref, &
198           precip_rain, precip_snow, snow, qsurf,  &
199           radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
200           1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
201           tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
202   
203     
204      !* Calculate snow height, run_off, age of snow
205
206      CALL fonte_neige( knon, is_ter, knindex, dtime, &
207           tsurf, precip_rain, precip_snow, &
208           snow, qsol, tsurf_new, evap, icesub)
209   
210      ! calculate the age of snow
211      CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
212   
213      WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
214   
215      DO i=1, knon
216         zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
217         alb_lim(i)  = alb_neig(i) *zfra(i) + alb_lim(i)*(1.0-zfra(i))
218      END DO
219
220     
221      !* Return albedo :
222      !    alb1_new and alb2_new are here given the same values
223     
224      alb1_new(:) = 0.0
225      alb2_new(:) = 0.0
226      alb1_new(1:knon) = alb_lim(1:knon)
227      alb2_new(1:knon) = alb_lim(1:knon)
228       
229      !* Send to coupler
230      !  The run-off from river and coast are not calculated in the bucket modele.
231      !  For testing purpose of the coupled modele we put the run-off to zero.
232      IF (type_ocean=='couple') THEN
233         dummy_riverflow(:)   = 0.0
234         dummy_coastalflow(:) = 0.0
235         CALL cpl_send_land_fields(itime, knon, knindex, &
236              dummy_riverflow, dummy_coastalflow)
237      ENDIF
238
239    !
240    ! ***  Flux aggregation ***
241    !
242    ELSE IF (iflag_hetero_surf == 2) THEN
243      ! initialize output tables
244      evap_tersrf(:,:)      = 0.
245      fluxsens_tersrf(:,:)  = 0.
246      fluxlat_tersrf(:,:)   = 0.
247      tsurf_new_tersrf(:,:) = 0.
248      dflux_s_tersrf(:,:)   = 0.
249      dflux_l_tersrf(:,:)   = 0.
250      radsol_tersrf(:,:)    = 0.
251      swnet_tersrf(:,:)     = 0.
252      lwnet_tersrf(:,:)     = 0.
253      ! hyp: surface emissivity = 1
254      emis_tersrf(:,:)      = 1.
255 
256      ! * calculate total absorbed radiance at surface
257      DO j=1, nbtersrf
258        ! SW
259        swnet_tersrf(klon_1D,j) = (1. - albedo_tersrf(klon_1D,j)) / (1. - alb_lim(klon_1D)) * swnet(klon_1D)
260        ! LW
261        ! first order
262        lwnet_tersrf(klon_1D,j) = lwnet(klon_1D) + 4. * emis_tersrf(klon_1D,j) * RSIGMA * tsurf(klon_1D)**3 * &
263                                  (tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))
264      ENDDO
265      ! LW second order corrections
266      !- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
267      IF (iflag_order2_sollw == 1) THEN
268        meansqT(:) = 0. ! as working buffer
269        !
270        DO j=1, nbtersrf
271          meansqT(klon_1D) = meansqT(klon_1D) + (tsurf_tersrf(klon_1D,j) - tsurf(klon_1D))**2 * frac_tersrf(klon_1D,j)
272        ENDDO
273        DO j=1, nbtersrf
274          lwnet_tersrf(klon_1D,j) = lwnet_tersrf(klon_1D,j) + 6. * RSIGMA * tsurf(klon_1D)**2 * (meansqT(klon_1D) - &
275                              (tsurf(klon_1D) - tsurf_tersrf(klon_1D,j))**2)
276        ENDDO
277      ENDIF
278      ! net radiation
279      radsol_tersrf(:,:) = swnet_tersrf(:,:) + lwnet_tersrf(:,:)
280 
281      ! * compute evapotranspiration coefficient
282      capsol(:) = 1.0/(2.5578E+06*0.15)
283      dif_grnd(:) = 0.
284
285      ! unused variables in cdrag routine
286      pblh(:) = 0.
287      ri_in(:) = 0.
288      sst(:) = 0.
289
290      ! Loop on sub-surfaces
291      DO j=1, nbtersrf
292        ! * drag coefficients
293        CALL cdrag(knon, is_ter, speed, temp_air, spechum, geop1, pref, pblh, &
294          tsurf_tersrf(:,j), qsurf_tersrf(:,j), z0m_tersrf(:,j), z0h_tersrf(:,j), ri_in, 0, &
295          cdragm_tersrf(:,j), cdragh_tersrf(:,j), zri_tersrf(:,j), plev, precip_rain, sst, p1lay)
296
297        ! * calculate temperature, heat capacity and conduction flux in soil
298        IF (soil_model) THEN
299          !
300          CALL soil_hetero(dtime, is_ter, knon, snow, tsurf_tersrf(:,j), qsol,  &
301            longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil_tersrf(:,:,j), soilcap, soilflux, &
302            inertie_tersrf(:,j), conv_ratio_tersrf(:,j))
303          !
304          cal(:) = RCPD / soilcap(:)
305          radsol_tersrf(:,j) = radsol_tersrf(:,j)  + soilflux(:)
306          !
307        ELSE
308          cal = RCPD * capsol
309          IF (klon_glo .EQ. 1) THEN
310            cal = 0.
311          ENDIF
312        ENDIF
313
314        ! * calcultaion of fluxes
315        CALL calcul_fluxs(knon, is_ter, dtime, &
316             tsurf_tersrf(klon_1D,j), p1lay, cal, beta_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), cdragh_tersrf(klon_1D,j), pref, &
317             precip_rain, precip_snow, snow, qsurf_tersrf(klon_1D,j),  &
318             radsol_tersrf(klon_1D,j), dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
319             1.,petAcoef, peqAcoef, petBcoef, peqBcoef, &
320             tsurf_new_tersrf(klon_1D,j), evap_tersrf(klon_1D,j), fluxlat_tersrf(klon_1D,j), fluxsens_tersrf(klon_1D,j), &
321             dflux_s_tersrf(klon_1D,j), dflux_l_tersrf(klon_1D,j))
322
323        ! if snow > 0
324        ! calculate snow height, run_off, age of snow
325        CALL fonte_neige( knon, is_ter, knindex, dtime, &
326             tsurf_tersrf(:,j), precip_rain, precip_snow, &
327             snow, qsol, tsurf_new_tersrf(:,j), evap_tersrf(:,j), icesub)
328
329      ENDDO ! loop on sub-surfaces
330
331      ! calculate the age of snow
332      CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
333      WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
334   
335      DO i=1, knon
336        zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
337        alb_lim(i)  = alb_neig(i) *zfra(i) + alb_lim(i)*(1.0-zfra(i))
338      END DO
339
340      ! return albedo :
341      !    alb1_new and alb2_new are here given the same values
342      alb1_new(:) = 0.0
343      alb2_new(:) = 0.0
344      alb1_new(1:knon) = alb_lim(1:knon)
345      alb2_new(1:knon) = alb_lim(1:knon)
346       
347      ! send to coupler
348      !  the run-off from river and coast are not calculated in the bucket modele.
349      !  for testing purpose of the coupled modele we put the run-off to zero.
350      IF (type_ocean=='couple') THEN
351         dummy_riverflow(:)   = 0.0
352         dummy_coastalflow(:) = 0.0
353         CALL cpl_send_land_fields(itime, knon, knindex, &
354              dummy_riverflow, dummy_coastalflow)
355      ENDIF
356
357      ! * average of fluxes and surface variables
358      qsurf = average_surf_var(klon, nbtersrf, qsurf_tersrf, frac_tersrf, 'ARI')
359      tsurf_new = average_surf_var(klon, nbtersrf, tsurf_new_tersrf, frac_tersrf, 'ARI')
360      evap = average_surf_var(klon, nbtersrf, evap_tersrf, frac_tersrf, 'ARI')
361      fluxlat = average_surf_var(klon, nbtersrf, fluxlat_tersrf, frac_tersrf, 'ARI')
362      fluxsens = average_surf_var(klon, nbtersrf, fluxsens_tersrf, frac_tersrf, 'ARI')
363      dflux_l = average_surf_var(klon, nbtersrf, dflux_l_tersrf, frac_tersrf, 'ARI')
364      dflux_s = average_surf_var(klon, nbtersrf, dflux_s_tersrf, frac_tersrf, 'ARI')
365      DO k=1, nsoilmx
366        tsoil(:,k) = average_surf_var(klon, nbtersrf, tsoil_tersrf(:,k,:), frac_tersrf, 'ARI')
367      ENDDO
368
369      ! order 2 correction to tsurf_new, for radiation computations (main atm effect of Ts)
370      IF (iflag_order2_sollw == 1) THEN
371        meansqT(:) = 0. ! as working buffer
372        DO j=1, nbtersrf
373          meansqT(klon_1D) = meansqT(klon_1D)+(tsurf_tersrf(klon_1D,j)-tsurf_new(klon_1D))**2 *frac_tersrf(klon_1D,j)
374        ENDDO
375        tsurf_new(:) = tsurf_new(:)+1.5*meansqT(:)/tsurf_new(:)
376      ENDIF   ! iflag_order2_sollw == 1
377
378    ENDIF ! iflag_hetero_surf
379!
380!* End
381!
382  END SUBROUTINE surf_land_bucket_hetero
383!
384!****************************************************************************************
385!
386END MODULE surf_land_bucket_hetero_mod
Note: See TracBrowser for help on using the repository browser.