source: LMDZ6/branches/Amaury_dev/libf/phylmd/surf_land_mod.F90 @ 5411

Last change on this file since 5411 was 5231, checked in by abarral, 3 months ago

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