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

Last change on this file since 5273 was 5273, checked in by abarral, 11 hours ago

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