source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/surf_land_mod.F90

Last change on this file was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

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