source: LMDZ6/trunk/libf/phylmd/surf_land_mod.F90 @ 5715

Last change on this file since 5715 was 5627, checked in by amaison, 2 months ago

Representation of heterogeneous continental subsurfaces with parameter or flux aggregation in the simplified surface model (bucket) for 1D case studies.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1!
2MODULE surf_land_mod
3
4  IMPLICIT NONE
5
6CONTAINS
7!
8!****************************************************************************************
9! 
10  SUBROUTINE surf_land(itime, dtime, date0, jour, knon, knindex, &
11       rlon, rlat, yrmu0, &
12       debut, lafin, zlev, ccanopy, swnet, lwnet, albedo, &
13       tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
14       AcoefH, AcoefQ, BcoefH, BcoefQ, &
15       AcoefU, AcoefV, BcoefU, BcoefV, &
16       pref, u1, v1, gustiness, rugoro, pctsrf, &
17       lwdown_m, q2m, t2m, &
18       snow, qsol, agesno, tsoil, &
19       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, fluxbs, &   
20       qsurf, tsurf_new, dflux_s, dflux_l, &
21       flux_u1, flux_v1 , &
22       veget,lai,height, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
23       cdragm_tersrf, cdragh_tersrf, &
24       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &
25#ifdef ISO
26       ,xtprecip_rain, xtprecip_snow,xtspechum, &
27       xtsnow, xtsol,xtevap,h1, &
28       runoff_diag,xtrunoff_diag,Rland_ice &
29#endif               
30               )
31
32    USE dimphy
33    USE surface_data, ONLY    : ok_veget
34    USE carbon_cycle_mod
35
36
37    ! See comments in each module surf_land_orchidee_xxx for compatiblity with ORCHIDEE
38#ifdef ORCHIDEE_NOOPENMP
39    ! Compilation with cpp key ORCHIDEE NOOPENMP
40    USE surf_land_orchidee_noopenmp_mod
41#else
42#if ORCHIDEE_NOZ0H
43    ! Compilation with cpp key ORCHIDEE NOZ0H
44    USE surf_land_orchidee_noz0h_mod
45#else
46#if ORCHIDEE_NOFREIN
47    ! Compilation with cpp key ORCHIDEE_NOFREIN
48    USE surf_land_orchidee_nofrein_mod
49#else
50#if ORCHIDEE_NOUNSTRUCT
51    ! Compilation with cpp key ORCHIDEE_NOUNSTRUCT
52    USE surf_land_orchidee_nounstruct_mod
53#else
54#if ORCHIDEE_NOLIC
55    ! Compilation with cpp key ORCHIDEE_NOLIC
56    USE surf_land_orchidee_nolic_mod
57#else
58    ! Default version
59    USE surf_land_orchidee_mod
60#endif
61#endif
62#endif
63#endif
64#endif
65   
66    USE surf_land_bucket_mod
67    USE surf_land_bucket_hetero_mod
68    USE calcul_fluxs_mod
69    USE indice_sol_mod
70#ifdef ISO
71    use infotrac_phy, ONLY: ntiso,niso
72    use isotopes_mod, ONLY: nudge_qsol, iso_eau
73#ifdef ISOVERIF
74    use isotopes_verif_mod
75#endif
76#endif
77
78USE dimpft_mod_h
79        USE clesphys_mod_h
80    USE yomcst_mod_h
81USE print_control_mod, ONLY: lunout
82    USE dimsoil_mod_h, ONLY: nsoilmx
83    USE compbl_mod_h
84
85! Input variables 
86!****************************************************************************************
87    INTEGER, INTENT(IN)                     :: itime, jour, knon
88    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
89    REAL, INTENT(IN)                        :: date0
90    REAL, DIMENSION(klon), INTENT(IN)       :: rlon, rlat
91    REAL, DIMENSION(klon), INTENT(IN)       :: yrmu0  ! cosine of solar zenith angle
92    LOGICAL, INTENT(IN)                     :: debut, lafin
93    REAL, INTENT(IN)                        :: dtime
94    REAL, DIMENSION(klon), INTENT(IN)       :: ccanopy
95    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
96    REAL, DIMENSION(klon), INTENT(IN)       :: albedo  ! albedo for whole short-wave interval
97    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
98    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
99    REAL, DIMENSION(klon), INTENT(IN)       :: cdragh, cdragm
100    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow, precip_bs
101    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
102    REAL, DIMENSION(klon), INTENT(IN)       :: AcoefH, AcoefQ, BcoefH, BcoefQ
103    REAL, DIMENSION(klon), INTENT(IN)       :: AcoefU, AcoefV, BcoefU, BcoefV
104    REAL, DIMENSION(klon), INTENT(IN)       :: pref   ! pressure reference
105    REAL, DIMENSION(klon), INTENT(IN)       :: u1, v1, gustiness
106    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
107    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
108    REAL, DIMENSION(klon), INTENT(IN)       :: lwdown_m  ! downwelling longwave radiation at mean surface
109                                                         ! corresponds to previous sollwdown
110    REAL, DIMENSION(klon), INTENT(IN)       :: q2m, t2m
111    REAL, DIMENSION(klon, nbtersrf), INTENT(IN) :: tsurf_tersrf
112#ifdef ISO
113    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
114    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtspechum
115#endif
116! In/Output variables
117!****************************************************************************************
118    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
119    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
120    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
121    REAL, DIMENSION(klon), INTENT(INOUT)          :: zlev
122    REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf
123#ifdef ISO
124    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsnow, xtsol
125#endif
126
127! Output variables
128!****************************************************************************************
129    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
130!albedo SB >>>
131!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new ! albdeo for shortwave interval 1(visible)
132!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new ! albedo for shortwave interval 2(near infrared)
133    REAL, DIMENSION(6), INTENT(IN) :: SFRWL
134    REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
135!albedo SB <<<
136    REAL, DIMENSION(klon), INTENT(OUT)       :: evap
137    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat, fluxbs
138    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
139    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
140    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l     
141    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1  ! flux for U and V at first model level
142    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: veget,lai
143    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: height
144! AM
145    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: tsurf_new_tersrf
146    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: qsurf_tersrf
147    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragm_tersrf
148    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: cdragh_tersrf
149    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: swnet_tersrf
150    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: lwnet_tersrf
151    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxsens_tersrf
152    REAL, DIMENSION(klon, nbtersrf), INTENT(OUT) :: fluxlat_tersrf
153#ifdef ISO
154    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
155    REAL, DIMENSION(klon), INTENT(OUT)      :: h1
156    REAL, DIMENSION(niso,klon), INTENT(OUT)      :: xtrunoff_diag
157    REAL, DIMENSION(klon), INTENT(OUT)      :: runoff_diag
158    REAL, DIMENSION(niso,klon), INTENT(IN)        :: Rland_ice
159#endif
160
161! Local variables
162!****************************************************************************************
163    REAL, DIMENSION(klon) :: p1lay_tmp
164    REAL, DIMENSION(klon) :: pref_tmp
165    REAL, DIMENSION(klon) :: swdown     ! downwelling shortwave radiation at land surface
166    REAL, DIMENSION(klon) :: epot_air           ! potential air temperature
167    REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
168    REAL, DIMENSION(klon) :: u0, v0     ! surface speed
169    REAL, DIMENSION(klon) :: precip_totsnow     ! total solid precip
170    INTEGER               :: i,j
171    CHARACTER (len = 20)  :: modname = 'surf_land'
172    CHARACTER (len = 100) :: abort_message
173
174!albedo SB >>>
175    REAL, DIMENSION(klon)      :: alb1_new,alb2_new
176!albedo SB <<<
177
178#ifdef ISO       
179      real, parameter :: t_coup = 273.15
180      real, dimension(klon) :: fqfonte_diag
181      real, dimension(klon) :: snow_evap_diag
182      real, dimension(klon) :: fqcalving_diag
183      integer :: ixt
184#endif
185!****************************************************************************************
186!Total solid precip
187
188IF (ok_bs) THEN
189precip_totsnow(:)=precip_snow(:)+precip_bs(:)
190ELSE
191precip_totsnow(:)=precip_snow(:)
192ENDIF
193!****************************************************************************************
194#ifdef ISO
195#ifdef ISOVERIF
196!        write(*,*) 'surf_land_mod 162'
197        do i=1,knon
198          if (iso_eau.gt.0) then
199            call iso_verif_egalite_choix(precip_snow(i), &
200     &          xtprecip_snow(iso_eau,i),'surf_land_mod 129', &
201     &          errmax,errmaxrel)
202            call iso_verif_egalite_choix(qsol(i), &
203     &          xtsol(iso_eau,i),'surf_land_mod 139', &
204     &          errmax,errmaxrel)
205          endif 
206        enddo
207#endif
208#ifdef ISOVERIF
209!       write(*,*) 'surf_land 169: ok_veget=',ok_veget
210        do i=1,knon
211         do ixt=1,ntiso
212           call iso_verif_noNaN(xtprecip_snow(ixt,i),'surf_land 146')
213         enddo
214        enddo
215#endif
216#endif
217
218
219!****************************************************************************************
220! Choice between call to vegetation model (ok_veget=true) or simple calculation below
221!
222!****************************************************************************************
223   IF (ok_veget) THEN
224!****************************************************************************************
225!  Call model sechiba in model ORCHIDEE
226!
227!****************************************************************************************
228       p1lay_tmp(:)      = 0.0
229       pref_tmp(:)       = 0.0
230       p1lay_tmp(1:knon) = p1lay(1:knon)/100.
231       pref_tmp(1:knon)  = pref(1:knon)/100.
232!
233!* Calculate incoming flux for SW and LW interval: swdown
234!
235       swdown(:) = 0.0
236       DO i = 1, knon
237          swdown(i) = swnet(i)/(1-albedo(i))
238       END DO
239!
240!* Calculate potential air temperature
241!
242       epot_air(:) = 0.0
243       DO i = 1, knon
244          epot_air(i) = RCPD*temp_air(i)*(pref(i)/p1lay(i))**RKAPPA
245       END DO
246
247#ifdef ISO
248      CALL abort_physic('surf_land_mod 220','isos pas prevus dans orchidee',1)
249#endif
250       ! temporary for keeping same results using lwdown_m instead of lwdown
251       CALL surf_land_orchidee(itime, dtime, date0, knon, &
252            knindex, rlon, rlat, yrmu0, pctsrf, &
253            debut, lafin, &
254            zlev,  u1, v1, gustiness, temp_air, spechum, epot_air, ccanopy, &
255            cdragh, AcoefH, AcoefQ, BcoefH, BcoefQ, &
256            precip_rain, precip_totsnow, lwdown_m, swnet, swdown, &
257            pref_tmp, q2m, t2m, &
258            evap, fluxsens, fluxlat,  &             
259            tsol_rad, tsurf_new, alb1_new, alb2_new, &
260            emis_new, z0m, z0h, qsurf, &
261            veget, lai, height &
262!#ifdef ISO
263!            , xtprecip_rain, xtprecip_snow, xtspechum, xtevap &
264!#endif
265            )                 
266
267#ifdef ISO
268#ifdef ISOVERIF
269     write(*,*) 'surf_land 193: apres surf_land_orchidee'   
270     do i=1,knon
271        if (iso_eau.gt.0) then
272             call iso_verif_egalite_choix(xtevap(iso_eau,i),evap(i), &
273    &            'surf_land 197',errmax,errmaxrel)
274        endif !if (iso_eau.gt.0) then     
275      enddo !do i=1,knon 
276#endif
277#endif
278! 
279!* Add contribution of relief to surface roughness
280! 
281       DO i=1,knon
282          z0m(i) = MAX(1.5e-05,SQRT(z0m(i)**2 + rugoro(i)**2))
283       ENDDO
284
285    ELSE  ! not ok_veget
286!****************************************************************************************
287! No extern vegetation model choosen, call simple bucket calculations instead.
288!
289!****************************************************************************************
290#ifdef ISO
291#ifdef ISOVERIF
292!       write(*,*) 'surf_land 247'
293        call iso_verif_egalite_vect1D( &
294     &           xtsnow,snow,'surf_land_mod 207',niso,klon)
295#endif
296#endif
297
298#ifdef ISO
299        if (nudge_qsol.eq.1) then
300          call surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)
301        endif
302        !write(*,*) 'surf_land 258'
303#endif
304      IF (iflag_hetero_surf .GT. 0) THEN
305        IF (klon .EQ. 1) THEN
306          !
307          CALL surf_land_bucket_hetero(itime, jour, knon, knindex, debut, dtime,&
308              tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
309              spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, pref, zlev, &
310              u1, v1, gustiness, rugoro, swnet, lwnet, &
311              snow, qsol, agesno, tsoil, &
312              qsurf, z0m, z0h, alb1_new, alb2_new, evap, &
313              fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l, &
314              tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
315              cdragm_tersrf, cdragh_tersrf, &
316              swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)
317        ELSE
318          abort_message = 'Heterogeneous continental subsurfaces (iflag_hetero_surf > 0) are only compatible in 1D cases.'
319          CALL abort_physic(modname,abort_message,1)
320        ENDIF
321      !
322      ELSE
323       CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
324            tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
325            spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, pref, &
326            u1, v1, gustiness, rugoro, swnet, lwnet, &
327            snow, qsol, agesno, tsoil, &
328            qsurf, z0m, alb1_new, alb2_new, evap, &
329            fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l &
330#ifdef ISO
331            ,xtprecip_rain, xtprecip_snow,xtspechum, &
332            xtsnow, xtsol,xtevap,h1, &
333     &      runoff_diag, xtrunoff_diag,Rland_ice &
334#endif           
335     &       )
336        z0h(1:knon)=z0m(1:knon) ! En attendant mieux
337
338      ENDIF ! iflag_hetero_surf
339
340    ENDIF ! ok_veget
341
342        ! blowing snow not treated yet over land
343        fluxbs(:)=0.
344
345
346!****************************************************************************************
347! Calculation for all land models
348! - Flux calculation at first modele level for U and V
349!****************************************************************************************
350! Suppose zero surface speed
351    u0(:)=0.0
352    v0(:)=0.0
353    CALL calcul_flux_wind(knon, dtime, &
354         u0, v0, u1, v1, gustiness, cdragm, &
355         AcoefU, AcoefV, BcoefU, BcoefV, &
356         p1lay, temp_air, &
357         flux_u1, flux_v1)
358
359#ifdef ISO
360#ifdef ISOVERIF
361!     write(*,*) 'surf_land 237: sortie'   
362      DO i=1,knon
363        IF (iso_eau >= 0) THEN
364             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
365    &            'surf_land 241',errmax,errmaxrel)
366        ENDIF !if (iso_eau.gt.0) then     
367      ENDDO !do i=1,knon 
368#endif
369#endif
370
371!albedo SB >>>
372     SELECT CASE(NSW)
373     CASE(2)
374       alb_dir_new(1:knon,1)=alb1_new(1:knon)
375       alb_dir_new(1:knon,2)=alb2_new(1:knon)
376     CASE(4)
377       alb_dir_new(1:knon,1)=alb1_new(1:knon)
378       alb_dir_new(1:knon,2)=alb2_new(1:knon)
379       alb_dir_new(1:knon,3)=alb2_new(1:knon)
380       alb_dir_new(1:knon,4)=alb2_new(1:knon)
381     CASE(6)
382       alb_dir_new(1:knon,1)=alb1_new(1:knon)
383       alb_dir_new(1:knon,2)=alb1_new(1:knon)
384       alb_dir_new(1:knon,3)=alb1_new(1:knon)
385       alb_dir_new(1:knon,4)=alb2_new(1:knon)
386       alb_dir_new(1:knon,5)=alb2_new(1:knon)
387       alb_dir_new(1:knon,6)=alb2_new(1:knon)
388     END SELECT
389
390     alb_dif_new=alb_dir_new
391!albedo SB <<<
392   
393  END SUBROUTINE surf_land
394
395
396#ifdef ISO
397  SUBROUTINE surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)
398
399    USE dimphy   
400    USE infotrac_phy, ONLY: niso
401    USE isotopes_mod, ONLY: region_nudge_qsol   
402    INTEGER, INTENT(IN)                       :: knon         
403    REAL, DIMENSION(klon), INTENT(IN)         :: rlon, rlat
404    REAL, DIMENSION(klon), INTENT(INOUT)      :: qsol
405    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex   
406    REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsol
407    REAL :: lat_min_nudge_qsol,lat_max_nudge_qsol
408    REAL :: lon_min_nudge_qsol,lon_max_nudge_qsol
409    INTEGER :: i,ixt
410    REAL :: qsol_new
411
412    IF (region_nudge_qsol == 1) THEN
413        ! Aamzonie du Sud
414        lat_min_nudge_qsol=-15.0
415        lat_max_nudge_qsol=-5.0
416        lon_min_nudge_qsol=-70.0
417        lon_max_nudge_qsol=-50.0
418    ELSE IF (region_nudge_qsol == 2) THEN
419        ! Aamzonie du Nord
420        lat_min_nudge_qsol=-5.0
421        lat_max_nudge_qsol=5.0
422        lon_min_nudge_qsol=-70.0
423        lon_max_nudge_qsol=-50.0
424    ELSE
425        WRITE(*,*) 'surf_land 298: cas pas prevu'
426        WRITE(*,*) 'region_nudge_qsol=',region_nudge_qsol
427        stop
428    ENDIF
429
430!    write(*,*) 'surf_land 314: knon=',knon
431!    write(*,*) 'rlat=',rlat
432!    write(*,*) 'rlon=',rlon
433!    write(*,*) 'region_nudge_qsol=',region_nudge_qsol
434
435    DO i=1,knon
436      IF ((rlat(knindex(i)) >= lat_min_nudge_qsol).and. &
437  &       (rlat(knindex(i)) <= lat_max_nudge_qsol).and. &
438  &       (rlon(knindex(i)) >= lon_min_nudge_qsol).and. &
439  &       (rlon(knindex(i)) <= lon_max_nudge_qsol)) THEN
440!        write(*,*) 'surf_land 324: bon domaine: rlat,rlon,qsol=', &
441!  &             rlat(knindex(i)),rlon(knindex(i)),qsol(knindex(i))
442        qsol_new=qsol(i)
443        IF (region_nudge_qsol == 1) THEN   
444           qsol_new=max(qsol(i),50.0)   
445        ELSE IF (region_nudge_qsol == 2) THEN     
446           qsol_new=max(qsol(i),120.0)
447        ELSE !if (region_nudge_qsol.eq.1) then
448           WRITE(*,*) 'surf_land 317: cas pas prevu'
449           WRITE(*,*) 'region_nudge_qsol=',region_nudge_qsol
450           STOP
451        ENDIF !if (region_nudge_qsol.eq.1) then
452        IF (qsol(i) > 0.0) THEN
453           DO ixt=1,niso
454              xtsol(ixt,i)=xtsol(ixt,i)*qsol_new/qsol(i)
455           ENDDO
456        ELSE !IF (qsol(i) > 0.0) THEN
457           DO ixt=1,niso
458             xtsol(ixt,i)=0.0
459           ENDDO
460        ENDIF !IF (qsol(i) > 0.0) THEN
461        qsol(i)=qsol_new
462        WRITE(*,*) 'surf_land 346: qsol_new=',qsol(i)     
463     ENDIF ! if ((rlat(i).ge.lat_min_nudge_qsol).and.
464  ENDDO !DO i=1,knon
465
466  END SUBROUTINE surf_land_nudge_qsol
467#endif
468
469!
470!****************************************************************************************
471! 
472END MODULE surf_land_mod
473!
474!****************************************************************************************
475! 
Note: See TracBrowser for help on using the repository browser.