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

Last change on this file since 4009 was 3927, checked in by Laurent Fairhead, 3 years ago

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