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

Last change on this file since 4523 was 4523, checked in by evignon, 17 months ago

merge de la branche blowing snow vers la trunk
premiere tentative
Etienne
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    ! See comments in each module surf_land_orchidee_xxx for compatiblity with ORCHIDEE
35#ifdef ORCHIDEE_NOOPENMP
36    ! Compilation with cpp key ORCHIDEE NOOPENMP
37    USE surf_land_orchidee_noopenmp_mod
38#else
39#if ORCHIDEE_NOZ0H
40    ! Compilation with cpp key ORCHIDEE NOZ0H
41    USE surf_land_orchidee_noz0h_mod
42#else
43#if ORCHIDEE_NOFREIN
44    ! Compilation with cpp key ORCHIDEE_NOFREIN
45    USE surf_land_orchidee_nofrein_mod
46#else
47#if ORCHIDEE_NOUNSTRUCT
48    ! Compilation with cpp key ORCHIDEE_NOUNSTRUCT
49    USE surf_land_orchidee_nounstruct_mod
50#else
51#if ORCHIDEE_NOLIC
52    ! Compilation with cpp key ORCHIDEE_NOLIC
53    USE surf_land_orchidee_nolic_mod
54#else
55    ! Default version#else
56    USE surf_land_orchidee_mod
57#endif
58#endif
59#endif
60#endif
61#endif
62
63    USE surf_land_bucket_mod
64    USE calcul_fluxs_mod
65    USE indice_sol_mod
66#ifdef ISO
67    use infotrac_phy, ONLY: ntiso,niso
68    use isotopes_mod, ONLY: nudge_qsol, iso_eau
69#ifdef ISOVERIF
70    use isotopes_verif_mod
71#endif
72#endif
73   
74    USE print_control_mod, ONLY: lunout
75
76    INCLUDE "dimsoil.h"
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
112! In/Output variables
113!****************************************************************************************
114    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
115    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
116    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
117#ifdef ISO
118    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsnow, xtsol
119#endif
120
121! Output variables
122!****************************************************************************************
123    REAL, DIMENSION(klon), INTENT(OUT)       :: z0m, z0h
124!albedo SB >>>
125!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new ! albdeo for shortwave interval 1(visible)
126!    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new ! albedo for shortwave interval 2(near infrared)
127    REAL, DIMENSION(6), INTENT(IN) :: SFRWL
128    REAL, DIMENSION(klon,nsw), INTENT(OUT)       :: alb_dir_new,alb_dif_new
129!albedo SB <<<
130    REAL, DIMENSION(klon), INTENT(OUT)       :: evap
131    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat, fluxbs
132    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
133    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
134    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l     
135    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1  ! flux for U and V at first model level
136    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: veget,lai
137    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: height
138#ifdef ISO
139    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
140    REAL, DIMENSION(klon), INTENT(OUT)      :: h1
141    REAL, DIMENSION(niso,klon), INTENT(OUT)      :: xtrunoff_diag
142    REAL, DIMENSION(klon), INTENT(OUT)      :: runoff_diag
143    REAL, DIMENSION(niso,klon), INTENT(IN)        :: Rland_ice
144#endif
145
146! Local variables
147!****************************************************************************************
148    REAL, DIMENSION(klon) :: p1lay_tmp
149    REAL, DIMENSION(klon) :: pref_tmp
150    REAL, DIMENSION(klon) :: swdown     ! downwelling shortwave radiation at land surface
151    REAL, DIMENSION(klon) :: epot_air           ! potential air temperature
152    REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
153    REAL, DIMENSION(klon) :: u0, v0     ! surface speed
154    REAL, DIMENSION(klon) :: precip_totsnow     ! total solid precip
155    INTEGER               :: i
156
157!albedo SB >>>
158    REAL, DIMENSION(klon)      :: alb1_new,alb2_new
159!albedo SB <<<
160
161#ifdef ISO       
162      real, parameter :: t_coup = 273.15
163      real, dimension(klon) :: fqfonte_diag
164      real, dimension(klon) ::  snow_evap_diag
165      real, dimension(klon) ::  fqcalving_diag
166      integer :: ixt
167#endif 
168
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_gcm('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,fluxbs, &             
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    ENDIF ! ok_veget
304
305        ! blowing snow not treated yet over land
306        fluxbs(:)=0.
307!****************************************************************************************
308! Calculation for all land models
309! - Flux calculation at first modele level for U and V
310!****************************************************************************************
311! Suppose zero surface speed
312    u0(:)=0.0
313    v0(:)=0.0
314    CALL calcul_flux_wind(knon, dtime, &
315         u0, v0, u1, v1, gustiness, cdragm, &
316         AcoefU, AcoefV, BcoefU, BcoefV, &
317         p1lay, temp_air, &
318         flux_u1, flux_v1)
319
320#ifdef ISO
321#ifdef ISOVERIF
322!     write(*,*) 'surf_land 237: sortie'   
323     do i=1,knon
324        if (iso_eau.gt.0) then
325             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
326    &            'surf_land 241',errmax,errmaxrel)
327        endif !if (iso_eau.gt.0) then     
328      enddo !do i=1,knon 
329#endif
330#endif
331
332!albedo SB >>>
333     SELECT CASE(NSW)
334     CASE(2)
335       alb_dir_new(1:knon,1)=alb1_new(1:knon)
336       alb_dir_new(1:knon,2)=alb2_new(1:knon)
337     CASE(4)
338       alb_dir_new(1:knon,1)=alb1_new(1:knon)
339       alb_dir_new(1:knon,2)=alb2_new(1:knon)
340       alb_dir_new(1:knon,3)=alb2_new(1:knon)
341       alb_dir_new(1:knon,4)=alb2_new(1:knon)
342     CASE(6)
343       alb_dir_new(1:knon,1)=alb1_new(1:knon)
344       alb_dir_new(1:knon,2)=alb1_new(1:knon)
345       alb_dir_new(1:knon,3)=alb1_new(1:knon)
346       alb_dir_new(1:knon,4)=alb2_new(1:knon)
347       alb_dir_new(1:knon,5)=alb2_new(1:knon)
348       alb_dir_new(1:knon,6)=alb2_new(1:knon)
349     END SELECT
350
351     alb_dif_new=alb_dir_new
352!albedo SB <<<
353   
354  END SUBROUTINE surf_land
355
356
357#ifdef ISO
358    SUBROUTINE surf_land_nudge_qsol(knon,rlat,rlon,qsol,xtsol,knindex)
359
360    USE dimphy   
361    USE infotrac_phy, ONLY: niso
362    use isotopes_mod, ONLY: region_nudge_qsol   
363    INTEGER, INTENT(IN)                     ::  knon         
364    REAL, DIMENSION(klon), INTENT(IN)       :: rlon, rlat
365    REAL, DIMENSION(klon), INTENT(INOUT)    :: qsol
366    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex   
367    REAL, DIMENSION(niso,klon), INTENT(INOUT)    :: xtsol
368    real lat_min_nudge_qsol,lat_max_nudge_qsol
369    real lon_min_nudge_qsol,lon_max_nudge_qsol
370    integer i,ixt
371    real qsol_new
372
373  if (region_nudge_qsol.eq.1) then
374        ! Aamzonie du Sud
375        lat_min_nudge_qsol=-15.0
376        lat_max_nudge_qsol=-5.0
377        lon_min_nudge_qsol=-70.0
378        lon_max_nudge_qsol=-50.0
379  else if (region_nudge_qsol.eq.2) then
380        ! Aamzonie du Nord
381        lat_min_nudge_qsol=-5.0
382        lat_max_nudge_qsol=5.0
383        lon_min_nudge_qsol=-70.0
384        lon_max_nudge_qsol=-50.0
385  else
386        write(*,*) 'surf_land 298: cas pas prevu'
387        write(*,*) 'region_nudge_qsol=',region_nudge_qsol
388        stop
389  endif
390
391!  write(*,*) 'surf_land 314: knon=',knon
392!  write(*,*) 'rlat=',rlat
393!  write(*,*) 'rlon=',rlon
394!  write(*,*) 'region_nudge_qsol=',region_nudge_qsol
395
396    do i=1,knon
397     if ((rlat(knindex(i)).ge.lat_min_nudge_qsol).and. &
398  &    (rlat(knindex(i)).le.lat_max_nudge_qsol).and. &
399  &    (rlon(knindex(i)).ge.lon_min_nudge_qsol).and. &
400  &    (rlon(knindex(i)).le.lon_max_nudge_qsol)) then
401!        write(*,*) 'surf_land 324: bon domaine: rlat,rlon,qsol=', &
402!  &             rlat(knindex(i)),rlon(knindex(i)),qsol(knindex(i))
403        qsol_new=qsol(i)
404        if (region_nudge_qsol.eq.1) then       
405           qsol_new=max(qsol(i),50.0)   
406        else if (region_nudge_qsol.eq.2) then     
407           qsol_new=max(qsol(i),120.0)
408        else !if (region_nudge_qsol.eq.1) then
409          write(*,*) 'surf_land 317: cas pas prevu'
410          write(*,*) 'region_nudge_qsol=',region_nudge_qsol
411          stop
412        endif !if (region_nudge_qsol.eq.1) then
413        if (qsol(i).gt.0.0) then
414          do ixt=1,niso
415           xtsol(ixt,i)=xtsol(ixt,i)*qsol_new/qsol(i)
416          enddo
417        else !if (qsol(i).gt.0.0) then
418          do ixt=1,niso
419           xtsol(ixt,i)=0.0
420          enddo
421        endif !if (qsol(i).gt.0.0) then
422        qsol(i)=qsol_new
423        write(*,*) 'surf_land 346: qsol_new=',qsol(i)     
424     endif ! if ((rlat(i).ge.lat_min_nudge_qsol).and.
425  enddo !do i=1,knon
426
427  END SUBROUTINE surf_land_nudge_qsol
428#endif
429
430!
431!****************************************************************************************
432! 
433END MODULE surf_land_mod
434!
435!****************************************************************************************
436! 
Note: See TracBrowser for help on using the repository browser.