source: LMDZ6/trunk/libf/phylmdiso/surf_land_mod.F90 @ 5943

Last change on this file since 5943 was 5943, checked in by yann meurdesoif, 3 weeks ago

Replace phylmdiso files
YM

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