source: LMDZ4/trunk/libf/phytherm/pbl_surface_mod.F90 @ 1080

Last change on this file since 1080 was 852, checked in by Laurent Fairhead, 17 years ago

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 48.4 KB
Line 
1!
2! $Header$
3!
4MODULE pbl_surface_mod
5!
6! Planetary Boundary Layer and Surface module
7!
8! This module manage the calculation of turbulent diffusion in the boundary layer
9! and all interactions towards the differents sub-surfaces.
10!
11!
12  USE dimphy
13  USE mod_phys_lmdz_para,  ONLY : mpi_size
14  USE ioipsl
15  USE surface_data,        ONLY : ocean, ok_veget
16  USE surf_land_mod,       ONLY : surf_land
17  USE surf_landice_mod,    ONLY : surf_landice
18  USE surf_ocean_mod,      ONLY : surf_ocean
19  USE surf_seaice_mod,     ONLY : surf_seaice
20  USE cpl_mod,             ONLY : gath2cpl
21  USE climb_hq_mod,        ONLY : climb_hq_down, climb_hq_up
22  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
23  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
24
25  IMPLICIT NONE
26
27! Declaration of variables saved in restart file
28  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: qsol
29  !$OMP THREADPRIVATE(qsol)
30  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder
31  !$OMP THREADPRIVATE(fder)
32  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow
33  !$OMP THREADPRIVATE(snow)
34  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf
35  !$OMP THREADPRIVATE(qsurf)
36  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: evap
37  !$OMP THREADPRIVATE(evap)
38  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: rugos
39  !$OMP THREADPRIVATE(rugos)
40  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: agesno
41  !$OMP THREADPRIVATE(agesno)
42  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil 
43  !$OMP THREADPRIVATE(ftsoil)
44
45CONTAINS
46!
47!****************************************************************************************
48!
49  SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,&
50       evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
51
52! This routine should be called after the restart file has been read.
53! This routine initialize the restart variables and does some validation tests
54! for the index of the different surfaces and tests the choice of type of ocean.
55
56    INCLUDE "indicesol.h"
57    INCLUDE "dimsoil.h"
58    INCLUDE "iniprint.h"
59 
60! Input variables
61!****************************************************************************************
62    REAL, DIMENSION(klon), INTENT(IN)                 :: qsol_rst
63    REAL, DIMENSION(klon), INTENT(IN)                 :: fder_rst
64    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: snow_rst
65    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: qsurf_rst
66    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: evap_rst
67    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: rugos_rst
68    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: agesno_rst
69    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
70
71 
72! Local variables
73!****************************************************************************************
74    INTEGER                       :: ierr
75    CHARACTER(len=80)             :: abort_message
76    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
77   
78
79!****************************************************************************************
80! Allocate and initialize module variables with fields read from restart file.
81!
82!****************************************************************************************   
83    ALLOCATE(qsol(klon), stat=ierr)
84    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
85
86    ALLOCATE(fder(klon), stat=ierr)
87    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
88
89    ALLOCATE(snow(klon,nbsrf), stat=ierr)
90    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
91
92    ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
93    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
94
95    ALLOCATE(evap(klon,nbsrf), stat=ierr)
96    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
97
98    ALLOCATE(rugos(klon,nbsrf), stat=ierr)
99    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
100
101    ALLOCATE(agesno(klon,nbsrf), stat=ierr)
102    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
103
104    ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
105    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
106
107
108    qsol(:)       = qsol_rst(:)
109    fder(:)       = fder_rst(:)
110    snow(:,:)     = snow_rst(:,:)
111    qsurf(:,:)    = qsurf_rst(:,:)
112    evap(:,:)     = evap_rst(:,:)
113    rugos(:,:)    = rugos_rst(:,:)
114    agesno(:,:)   = agesno_rst(:,:)
115    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
116
117
118!****************************************************************************************
119! Test for sub-surface indices
120!
121!****************************************************************************************
122    IF (is_ter /= 1) THEN
123      WRITE(lunout,*)" *** Warning ***"
124      WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
125      WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
126      abort_message="voir ci-dessus"
127      CALL abort_gcm(modname,abort_message,1)
128    ENDIF
129
130    IF ( is_oce > is_sic ) THEN
131      WRITE(lunout,*)' *** Warning ***'
132      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
133      WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
134      WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
135      abort_message='voir ci-dessus'
136      CALL abort_gcm(modname,abort_message,1)
137    ENDIF
138
139    IF ( is_lic > is_sic ) THEN
140      WRITE(lunout,*)' *** Warning ***'
141      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
142      WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
143      WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
144      abort_message='voir ci-dessus'
145      CALL abort_gcm(modname,abort_message,1)
146    ENDIF
147
148!****************************************************************************************
149! Validation of ocean mode
150!
151!****************************************************************************************
152
153    IF (ocean /= 'slab  ' .AND. ocean /= 'force ' .AND. ocean /= 'couple') THEN
154      WRITE(lunout,*)' *** Warning ***'
155      WRITE(lunout,*)'Option couplage pour l''ocean = ', ocean
156      abort_message='option pour l''ocean non valable'
157      CALL abort_gcm(modname,abort_message,1)
158    ENDIF
159
160!****************************************************************************************
161! Test of coherence between variable ok_veget and cpp key CPP_VEGET
162!
163!****************************************************************************************
164    IF (ok_veget) THEN
165#ifndef CPP_VEGET
166       abort_message='Pb de coherence: ok_veget = .true. mais CPP_VEGET = .false.'
167       CALL abort_gcm(modname,abort_message,1)
168#endif
169    ENDIF
170
171
172  END SUBROUTINE pbl_surface_init
173
174!****************************************************************************************
175
176
177  SUBROUTINE pbl_surface( &
178       dtime,     date0,     itap,     jour,          &
179       debut,     lafin,                              &
180       rlon,      rlat,      rugoro,   rmu0,          &
181       rain_f,    snow_f,    solsw_m,  sollw_m,       &
182       t,         q,         u,        v,             &
183       pplay,     paprs,     pctsrf,                  &
184       ts,        albe,      alblw,    u10m,   v10m,  &
185       sollwdown, cdragh,    cdragm,   zu1,    zv1,   &
186       albsol,    albsollw,  zxsens,   zxevap,        &
187       zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
188       d_t,       d_q,       d_u,      d_v,           &
189       zcoefh,    pctsrf_new,                         &
190       qsol_d,    zq2m,      s_pblh,   s_plcl,        &
191       s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,        &
192       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
193       zxrugs,    zu10m,     zv10m,    fder_print,    &
194       zxqsurf,   rh2m,      zxfluxu,  zxfluxv,       &
195       rugos_d,   agesno_d,  sollw,    solsw,         &
196       d_ts,      evap_d,    fluxlat,  t2m,           &
197       wfbils,    wfbilo,    flux_t,   flux_u, flux_v,&
198       dflux_t,   dflux_q,   zxsnow,                  &
199       zxfluxt,   zxfluxq,   q2m,      flux_q, tke    )
200!****************************************************************************************
201! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
202! Objet: interface de "couche limite" (diffusion verticale)
203!
204!AA REM:
205!AA-----
206!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
207!AA pour l'instant le calcul de la couche limite pour les traceurs
208!AA se fait avec cltrac et ne tient pas compte de la differentiation
209!AA des sous-fraction de sol.
210!AA REM bis :
211!AA----------
212!AA Pour pouvoir extraire les coefficient d'echanges et le vent
213!AA dans la premiere couche, 3 champs supplementaires ont ete crees
214!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
215!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
216!AA si les informations des subsurfaces doivent etre prises en compte
217!AA il faudra sortir ces memes champs en leur ajoutant une dimension,
218!AA c'est a dire nbsrf (nbre de subsurface).
219!
220! Arguments:
221!
222! dtime----input-R- interval du temps (secondes)
223! itap-----input-I- numero du pas de temps
224! date0----input-R- jour initial
225! t--------input-R- temperature (K)
226! q--------input-R- vapeur d'eau (kg/kg)
227! u--------input-R- vitesse u
228! v--------input-R- vitesse v
229! ts-------input-R- temperature du sol (en Kelvin)
230! paprs----input-R- pression a intercouche (Pa)
231! pplay----input-R- pression au milieu de couche (Pa)
232! rlat-----input-R- latitude en degree
233! rugos----input-R- longeur de rugosite (en m)
234!
235! d_t------output-R- le changement pour "t"
236! d_q------output-R- le changement pour "q"
237! d_u------output-R- le changement pour "u"
238! d_v------output-R- le changement pour "v"
239! d_ts-----output-R- le changement pour "ts"
240! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
241!                    (orientation positive vers le bas)
242! tke---input/output-R- tke (kg/m**2/s)
243! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
244! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
245! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
246! dflux_t--output-R- derive du flux sensible
247! dflux_q--output-R- derive du flux latent
248! zu1------output-R- le vent dans la premiere couche
249! zv1------output-R- le vent dans la premiere couche
250! trmb1----output-R- deep_cape
251! trmb2----output-R- inhibition
252! trmb3----output-R- Point Omega
253! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
254! plcl-----output-R- Niveau de condensation
255! pblh-----output-R- HCL
256! pblT-----output-R- T au nveau HCL
257!
258    INCLUDE "indicesol.h"
259    INCLUDE "dimsoil.h"
260    INCLUDE "YOMCST.h"
261    INCLUDE "iniprint.h"
262    INCLUDE "FCTTRE.h"
263    INCLUDE "clesphys.h"
264    INCLUDE "compbl.h"
265    INCLUDE "dimensions.h"
266    INCLUDE "YOETHF.h"
267    INCLUDE "temps.h"
268    INCLUDE "control.h"
269
270! Input variables
271!****************************************************************************************
272    REAL,                         INTENT(IN)        :: dtime
273    REAL,                         INTENT(IN)        :: date0
274    INTEGER,                      INTENT(IN)        :: itap
275    INTEGER,                      INTENT(IN)        :: jour    ! jour de l'annee en cours
276    LOGICAL,                      INTENT(IN)        :: debut, lafin
277    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon, rlat
278    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro
279    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosinus de l'angle solaire zenithal
280    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f, snow_f
281    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! mean value
282    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! mean value
283    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t, q
284    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u, v
285    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay
286    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs
287    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf
288
289! Input/Output variables
290!****************************************************************************************
291    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts
292    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: albe
293    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alblw
294    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m, v10m
295
296! Output variables
297!****************************************************************************************
298    REAL, DIMENSION(klon),        INTENT(OUT)       :: sollwdown
299    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh, cdragm
300    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1
301    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1
302    REAL, DIMENSION(klon),        INTENT(OUT)       :: albsol
303    REAL, DIMENSION(klon),        INTENT(OUT)       :: albsollw
304    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens, zxevap
305    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol
306    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat
307    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m
308    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
309    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t, d_q
310    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u, d_v
311    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefh
312    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pctsrf_new
313
314! Output only for diagnostics
315    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol_d
316    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m
317    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh
318    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl
319    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL
320    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL
321    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL
322    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT
323    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm
324    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1
325    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2
326    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3
327    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxrugs
328    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m
329    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m
330    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print
331    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf
332    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m
333    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu, zxfluxv
334    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: rugos_d
335    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: agesno_d
336    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw, solsw
337    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts
338    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: evap_d
339    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat
340    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m
341    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils, wfbilo
342    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t
343    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u, flux_v
344
345! Output not needed
346    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t, dflux_q
347    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow
348    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt, zxfluxq
349    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m
350    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q
351
352! Input/output
353    REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke
354
355
356! Local variables with attribute SAVE
357!****************************************************************************************
358    INTEGER, SAVE                            :: nhoridbg, nidbg
359!$OMP THREADPRIVATE(nhoridbg, nidbg)
360    LOGICAL, SAVE                            :: debugindex=.FALSE.
361!$OMP THREADPRIVATE(debugindex)
362    LOGICAL, SAVE                            :: first_call=.TRUE.
363!$OMP THREADPRIVATE(first_call)
364    CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
365!$OMP THREADPRIVATE(cl_surf)
366
367! Other local variables
368!****************************************************************************************
369    INTEGER                            :: i, k, nsrf
370    INTEGER                            :: knon, j
371    INTEGER                            :: idayref
372    INTEGER , DIMENSION(klon)          :: ni
373    REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
374    REAL                               :: amn, amx
375    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
376    REAL, DIMENSION(klon)              :: yts, yrugos, ypct, yz0_new
377    REAL, DIMENSION(klon)              :: yalb
378    REAL, DIMENSION(klon)              :: yalblw
379    REAL, DIMENSION(klon)              :: yu1, yv1
380    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
381    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
382    REAL, DIMENSION(klon)              :: ysollw, ysolsw, ysollwdown
383    REAL, DIMENSION(klon)              :: yfder
384    REAL, DIMENSION(klon)              :: yrads,yrugoro
385    REAL, DIMENSION(klon)              :: yfluxlat
386    REAL, DIMENSION(klon)              :: y_d_ts
387    REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
388    REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
389    REAL, DIMENSION(klon)              :: u1lay, v1lay
390    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
391    REAL, DIMENSION(klon)              :: yustar
392    REAL, DIMENSION(klon)              :: yu10mx
393    REAL, DIMENSION(klon)              :: yu10my
394    REAL, DIMENSION(klon)              :: ywindsp
395    REAL, DIMENSION(klon)              :: yt10m, yq10m
396    REAL, DIMENSION(klon)              :: ypblh
397    REAL, DIMENSION(klon)              :: ylcl
398    REAL, DIMENSION(klon)              :: ycapCL
399    REAL, DIMENSION(klon)              :: yoliqCL
400    REAL, DIMENSION(klon)              :: ycteiCL
401    REAL, DIMENSION(klon)              :: ypblT
402    REAL, DIMENSION(klon)              :: ytherm
403    REAL, DIMENSION(klon)              :: ytrmb1
404    REAL, DIMENSION(klon)              :: ytrmb2
405    REAL, DIMENSION(klon)              :: ytrmb3
406    REAL, DIMENSION(klon)              :: uzon, vmer
407    REAL, DIMENSION(klon)              :: tair1, qair1, tairsol
408    REAL, DIMENSION(klon)              :: psfce, patm
409    REAL, DIMENSION(klon)              :: qairsol, zgeo1
410    REAL, DIMENSION(klon)              :: rugo1
411    REAL, DIMENSION(klon)              :: yfluxsens, swdown
412    REAL, DIMENSION(klon)              :: petAcoef, peqAcoef, petBcoef, peqBcoef
413    REAL, DIMENSION(klon)              :: ypsref, epot_air
414    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb_new
415    REAL, DIMENSION(klon)              :: pctsrf_nsrf
416    REAL, DIMENSION(klon)              :: ztsol
417    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q
418    REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
419    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
420    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
421    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm
422    REAL, DIMENSION(klon,klev)         :: yu, yv
423    REAL, DIMENSION(klon,klev)         :: yt, yq
424    REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
425    REAL, DIMENSION(klon,klev)         :: delp
426    REAL, DIMENSION(klon,klev+1)       :: ypaprs
427    REAL, DIMENSION(klon,klev+1)       :: ytke
428    REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
429    REAL, DIMENSION(klon,nbsrf)        :: pctsrf_pot
430    CHARACTER(len=80)                  :: abort_message
431    CHARACTER(len=20)                  :: modname = 'pbl_surface'
432    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
433    LOGICAL, PARAMETER                 :: check=.FALSE.
434
435! For debugging with IOIPSL
436    INTEGER, DIMENSION(iim*(jjm+1))    :: ndexbg
437    REAL                               :: zjulian
438    REAL, DIMENSION(klon)              :: tabindx
439    REAL, DIMENSION(iim,jjm+1)         :: zx_lon, zx_lat
440    REAL, DIMENSION(iim,jjm+1)         :: debugtab
441
442
443    REAL, DIMENSION(klon,nbsrf)        :: pblh
444    REAL, DIMENSION(klon,nbsrf)        :: plcl
445    REAL, DIMENSION(klon,nbsrf)        :: capCL
446    REAL, DIMENSION(klon,nbsrf)        :: oliqCL
447    REAL, DIMENSION(klon,nbsrf)        :: cteiCL
448    REAL, DIMENSION(klon,nbsrf)        :: pblT
449    REAL, DIMENSION(klon,nbsrf)        :: therm
450    REAL, DIMENSION(klon,nbsrf)        :: trmb1
451    REAL, DIMENSION(klon,nbsrf)        :: trmb2
452    REAL, DIMENSION(klon,nbsrf)        :: trmb3
453    REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
454    REAL, DIMENSION(klon,nbsrf)        :: zx_qs1, zx_t1
455    REAL, DIMENSION(klon,nbsrf)        :: zdelta1, zcor1
456
457
458!jg+ temporary test
459    REAL, DIMENSION(klon,klev)         :: y_flux_u_old, y_flux_v_old
460    REAL, DIMENSION(klon,klev)         :: y_d_u_old, y_d_v_old
461!jg-
462   
463!****************************************************************************************
464! End of declarations
465!****************************************************************************************
466
467
468!****************************************************************************************
469! 1) Initialisation and validation tests
470!    Only done first time entering this subroutine
471!
472!****************************************************************************************
473
474
475    IF (first_call) THEN
476       first_call=.FALSE.
477     
478       ! Initilize debug IO
479       IF (debugindex .AND. mpi_size==1) THEN
480          ! initialize IOIPSL output
481          idayref = day_ini
482          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
483          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
484          DO i = 1, iim
485             zx_lon(i,1) = rlon(i+1)
486             zx_lon(i,jjm+1) = rlon(i+1)
487          ENDDO
488          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
489          CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
490               1,iim,1,jjm+1, &
491               itau_phy,zjulian,dtime,nhoridbg,nidbg)
492          ! no vertical axis
493          cl_surf(1)='ter'
494          cl_surf(2)='lic'
495          cl_surf(3)='oce'
496          cl_surf(4)='sic'
497          DO nsrf=1,nbsrf
498             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
499                  jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
500          END DO
501
502          CALL histend(nidbg)
503          CALL histsync(nidbg)
504
505       END IF
506       
507    ENDIF
508         
509!****************************************************************************************
510! 2) Initialization to zero
511!    Done for all local variables that will be compressed later
512!    and argument with INTENT(OUT)
513!****************************************************************************************
514    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
515    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0     ; zu1 = 0.0       
516    zv1 = 0.0     ; yqsurf = 0.0     ; yalb = 0.0      ; yalblw = 0.0   
517    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0   
518    ysollw = 0.0  ; ysollwdown = 0.0 ; yrugos = 0.0    ; yu1 = 0.0   
519    yv1 = 0.0     ; yrads = 0.0      ; ypaprs = 0.0    ; ypplay = 0.0
520    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
521    yq = 0.0      ; pctsrf_new = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
522    yrugoro = 0.0 ; yu10mx = 0.0     ; yu10my = 0.0    ; ywindsp = 0.0   
523    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
524    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
525    d_u = 0.0     ; d_v = 0.0        ; zcoefh = 0.0    ; yqsol = 0.0   
526    ytherm = 0.0  ; ytke=0.
527     
528    ytsoil = 999999.
529
530!****************************************************************************************
531! 3) - Calculate pressure thickness of each layer
532!    - Calculate the wind at first layer
533!
534!****************************************************************************************
535    DO k = 1, klev
536       DO i = 1, klon
537          delp(i,k) = paprs(i,k)-paprs(i,k+1)
538       ENDDO
539    ENDDO
540    DO i = 1, klon
541       zx_alf1 = 1.0
542       zx_alf2 = 1.0 - zx_alf1
543       u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
544       v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
545    ENDDO
546
547
548!****************************************************************************************
549! Test for rugos........ from physiq.. A la fin plutot???
550! Calcul de l'abedo moyen par maille
551!****************************************************************************************
552
553    zxrugs(:) = 0.0
554    DO nsrf = 1, nbsrf
555       DO i = 1, klon
556          rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
557          zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
558       ENDDO
559    ENDDO
560
561! Calcul de l'abedo moyen par maille
562    albsol(:)   = 0.0
563    albsollw(:) = 0.0
564    DO nsrf = 1, nbsrf
565       DO i = 1, klon
566          albsol(i)   = albsol(i)   + albe(i,nsrf)  * pctsrf(i,nsrf)
567          albsollw(i) = albsollw(i) + alblw(i,nsrf) * pctsrf(i,nsrf)
568       ENDDO
569    ENDDO
570
571
572
573! Calcule de ztsol (aussi fait dans physiq.F, pourrait etre un argument)
574    ztsol(:) = 0.0
575    DO nsrf = 1, nbsrf
576       DO i = 1, klon
577          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
578       ENDDO
579    ENDDO
580
581
582! Repartition du longwave par sous-surface linearisee
583    DO nsrf = 1, nbsrf
584       DO i = 1, klon
585          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
586          solsw(i,nsrf) = solsw_m(i)*(1.-albe(i,nsrf))/(1.-albsol(i))
587       ENDDO
588    ENDDO
589
590
591    DO i = 1, klon
592       sollwdown(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
593    ENDDO
594
595!****************************************************************************************
596! 4) Loop over different surfaces
597!
598! All points with a possibility of the current surface are used. This is done
599! to allow the sea-ice to appear or disappear. It is considered here that the
600! entier domaine of the ocean possibly can contain sea-ice.
601!
602!****************************************************************************************
603
604    pctsrf_pot = pctsrf
605    pctsrf_pot(:,is_oce) = 1. - zmasq(:)
606    pctsrf_pot(:,is_sic) = 1. - zmasq(:)
607     
608    loop_nbsrf: DO nsrf = 1, nbsrf
609
610! Search for index(ni) and size(knon) of domaine to treat
611       ni(:) = 0
612       knon  = 0
613       DO i = 1, klon
614          IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
615             knon = knon + 1
616             ni(knon) = i
617          ENDIF
618       ENDDO
619
620       ! write index, with IOIPSL
621       IF (debugindex .AND. mpi_size==1) THEN
622          tabindx(:)=0.
623          DO i=1,knon
624             tabindx(i)=FLOAT(i)
625          END DO
626          debugtab(:,:) = 0.
627          ndexbg(:) = 0
628          CALL gath2cpl(tabindx,debugtab,knon,ni)
629          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
630       ENDIF
631       
632!****************************************************************************************
633! 5) Compress variables
634!
635!****************************************************************************************
636
637       DO j = 1, knon
638          i = ni(j)
639          ypct(j) = pctsrf(i,nsrf)
640          yts(j) = ts(i,nsrf)
641          ysnow(j) = snow(i,nsrf)
642          yqsurf(j) = qsurf(i,nsrf)
643          yalb(j) = albe(i,nsrf)
644          yalblw(j) = alblw(i,nsrf)
645          yrain_f(j) = rain_f(i)
646          ysnow_f(j) = snow_f(i)
647          yagesno(j) = agesno(i,nsrf)
648          yfder(j) = fder(i)
649          ysolsw(j) = solsw(i,nsrf)
650          ysollw(j) = sollw(i,nsrf)
651          ysollwdown(j) = sollwdown(i)
652          yrugos(j) = rugos(i,nsrf)
653          yrugoro(j) = rugoro(i)
654          yu1(j) = u1lay(i)
655          yv1(j) = v1lay(i)
656          yrads(j) =  ysolsw(j)+ ysollw(j)
657          ypaprs(j,klev+1) = paprs(i,klev+1)
658          yu10mx(j) = u10m(i,nsrf)
659          yu10my(j) = v10m(i,nsrf)
660          ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )
661       END DO
662
663       DO k = 1, klev
664          DO j = 1, knon
665             i = ni(j)
666             ypaprs(j,k) = paprs(i,k)
667             ypplay(j,k) = pplay(i,k)
668             ydelp(j,k) = delp(i,k)
669             ytke(j,k)=tke(i,k,nsrf)
670             yu(j,k) = u(i,k)
671             yv(j,k) = v(i,k)
672             yt(j,k) = t(i,k)
673             yq(j,k) = q(i,k)
674          ENDDO
675       ENDDO
676       
677       DO k = 1, nsoilmx
678          DO j = 1, knon
679             i = ni(j)
680             ytsoil(j,k) = ftsoil(i,k,nsrf)
681          END DO
682       END DO
683       
684       ! qsol(water height in soil) only for bucket continental model
685       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
686          DO j = 1, knon
687             i = ni(j)
688             yqsol(j) = qsol(i)
689          END DO
690       ENDIF
691       
692!****************************************************************************************
693! 6) Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
694!    atmosphere and coefficients for turbulent diffusion at surface(Cdrag).
695!
696!****************************************************************************************
697
698       CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
699            ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf,  &
700            ycoefm, ycoefh,ytke)
701       
702!****************************************************************************************
703!
704! 8) "La descente" - "The downhill"
705
706!  climb_hq_down and climb_wind_down calculate the coefficients
707!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
708!  Only the coefficients at surface for H and Q are returned.
709!
710!****************************************************************************************
711
712! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
713       CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
714            ydelp, yt, yq, dtime, &
715            petAcoef, peqAcoef, petBcoef, peqBcoef)
716
717! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
718       CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv)
719     
720
721!****************************************************************************************
722! 9) Small calculations
723!
724!****************************************************************************************
725         
726       ypsref(:) = ypaprs(:,1) 
727       epot_air(:) = 0.0
728       epot_air(1:knon) = RCPD*yt(1:knon,1)*(ypsref(1:knon)/ypplay(1:knon,1))**RKAPPA
729
730       swdown(:) = 0.0
731       IF (nsrf .EQ. is_ter) THEN
732          swdown(1:knon) = ysolsw(1:knon)/(1-yalb(1:knon))
733       ELSE
734          swdown(1:knon) = ysolsw(1:knon)
735       ENDIF
736
737       ! constant CO2
738       r_co2_ppm(:) = co2_ppm
739
740!****************************************************************************************
741!
742! 10) Switch selon current surface
743!     It is necessary to start with the continental surfaces because the ocean
744!     needs their run-off.
745!
746!****************************************************************************************
747       SELECT CASE(nsrf)
748     
749       CASE(is_ter)
750          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
751               rlon, rlat, &
752               debut, lafin, ydelp(:,1), epot_air, r_co2_ppm, ysollwdown, ysolsw, swdown, &
753               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
754               petAcoef, peqAcoef, petBcoef, peqBcoef, &
755               ypsref, yu1, yv1, yrugoro, pctsrf, &
756               yrads, ysnow, yqsurf, yqsol, yagesno, &
757               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
758               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
759     
760       CASE(is_lic)
761          CALL surf_landice(itap, dtime, knon, ni, &
762               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
763               petAcoef, peqAcoef, petBcoef, peqBcoef, &
764               ypsref, yu1, yv1, yrugoro, pctsrf, &
765               yrads, ysnow, yqsurf, yqsol, yagesno, &
766               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
767               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
768         
769       CASE(is_oce)
770          CALL surf_ocean(rlon, rlat, ysollw, yalb, &
771               yrugos, ywindsp, rmu0, &
772               yfder, &
773               itap, dtime, jour, knon, ni, &
774               debut, swdown, &
775               ypplay(:,1), ycoefh(:,1), ycoefm(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
776               petAcoef, peqAcoef, petBcoef, peqBcoef, &
777               ypsref, yu1, yv1, yrugoro, pctsrf, &
778               yrads, ysnow, yqsurf, yagesno, &
779               yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
780               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
781         
782       CASE(is_sic)
783          CALL surf_seaice( &
784               rlon, rlat, ysollw, yalb, &
785               yfder, &
786               itap, dtime, jour, knon, ni, &
787               debut, lafin, swdown, &
788               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
789               petAcoef, peqAcoef, petBcoef, peqBcoef, &
790               ypsref, yu1, yv1, yrugoro, pctsrf, &
791               yrads, ysnow, yqsurf, yqsol, yagesno, &
792               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
793               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
794         
795
796       CASE DEFAULT
797          WRITE(lunout,*) 'Surface index = ', nsrf
798          abort_message = 'Surface index not valid'
799          CALL abort_gcm(modname,abort_message,1)
800       END SELECT
801
802!****************************************************************************************
803! Save the fraction of this sub-surface
804!
805!****************************************************************************************
806       pctsrf_new(:,nsrf) = pctsrf_nsrf(:)
807
808!****************************************************************************************
809! 11) - Calcul the increment of surface temperature
810!     - Update albedo
811!
812!****************************************************************************************
813       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
814 
815       yalb(1:knon) = yalb_new(1:knon)
816
817!****************************************************************************************
818!
819! 12) "La remontee" - "The uphill"
820!
821!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
822!  for X=H, Q, U and V, for all vertical levels.
823!
824!****************************************************************************************
825! H and Q
826       y_flux_t1(:) =  yfluxsens(:)
827       y_flux_q1(:) = -yevap(:)
828
829       CALL climb_hq_up(knon, dtime, yt, yq, &
830            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
831            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))   
832       
833! U and V
834       CALL climb_wind_up(knon, dtime, yu, yv, &
835            y_flux_u, y_flux_v, y_d_u, y_d_v)
836
837       DO j = 1, knon
838          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
839          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
840          yu1(j) = yu1(j) *  ypct(j)
841          yv1(j) = yv1(j) *  ypct(j)
842       ENDDO
843
844!****************************************************************************************
845! 13) Transform variables for output format :
846!     - Decompress
847!     - Multiply with pourcentage of current surface
848!     - Cumulate in global variable
849!
850!****************************************************************************************
851
852       tke(:,:,nsrf)=0.
853       DO k = 1, klev
854          DO j = 1, knon
855             i = ni(j)
856             ycoefh(j,k) = ycoefh(j,k) * ypct(j)
857             ycoefm(j,k) = ycoefm(j,k) * ypct(j)
858             y_d_t(j,k) = y_d_t(j,k) * ypct(j)
859             y_d_q(j,k) = y_d_q(j,k) * ypct(j)
860             y_d_u(j,k) = y_d_u(j,k) * ypct(j)
861             y_d_v(j,k) = y_d_v(j,k) * ypct(j)
862
863             flux_t(i,k,nsrf) = y_flux_t(j,k)
864             flux_q(i,k,nsrf) = y_flux_q(j,k)
865             flux_u(i,k,nsrf) = y_flux_u(j,k)
866             flux_v(i,k,nsrf) = y_flux_v(j,k)
867
868             tke(i,k,nsrf)=ytke(j,k)
869
870          ENDDO
871       ENDDO
872       
873       evap(:,nsrf) = - flux_q(:,1,nsrf)
874       
875       albe(:, nsrf) = 0.
876       alblw(:, nsrf) = 0.
877       snow(:, nsrf) = 0.
878       qsurf(:, nsrf) = 0.
879       rugos(:, nsrf) = 0.
880       fluxlat(:,nsrf) = 0.
881       DO j = 1, knon
882          i = ni(j)
883          d_ts(i,nsrf) = y_d_ts(j)
884          albe(i,nsrf) = yalb(j) 
885          alblw(i,nsrf) = yalblw(j)
886          snow(i,nsrf) = ysnow(j) 
887          qsurf(i,nsrf) = yqsurf(j)
888          rugos(i,nsrf) = yz0_new(j)
889          fluxlat(i,nsrf) = yfluxlat(j)
890          agesno(i,nsrf) = yagesno(j) 
891          cdragh(i) = cdragh(i) + ycoefh(j,1)
892          cdragm(i) = cdragm(i) + ycoefm(j,1)
893          dflux_t(i) = dflux_t(i) + y_dflux_t(j)
894          dflux_q(i) = dflux_q(i) + y_dflux_q(j)
895          zu1(i) = zu1(i) + yu1(j)
896          zv1(i) = zv1(i) + yv1(j)
897       END DO
898
899       IF ( nsrf .EQ. is_ter ) THEN
900          DO j = 1, knon
901             i = ni(j)
902             qsol(i) = yqsol(j)
903          END DO
904       END IF
905       
906       ftsoil(:,:,nsrf) = 0.
907       DO k = 1, nsoilmx
908          DO j = 1, knon
909             i = ni(j)
910             ftsoil(i, k, nsrf) = ytsoil(j,k)
911          END DO
912       END DO
913       
914       
915#ifdef CRAY
916       DO k = 1, klev
917          DO j = 1, knon
918             i = ni(j)
919#else
920       DO j = 1, knon
921          i = ni(j)
922          DO k = 1, klev
923#endif
924             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
925             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
926             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
927             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
928             zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
929#ifdef CRAY
930          END DO
931       END DO
932#else
933          END DO
934       END DO
935#endif
936
937!****************************************************************************************
938! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
939!     Call HBTM
940!
941!****************************************************************************************
942       t2m(:,nsrf)    = 0.
943       q2m(:,nsrf)    = 0.
944       u10m(:,nsrf)   = 0.
945       v10m(:,nsrf)   = 0.
946
947       pblh(:,nsrf)   = 0.        ! Hauteur de couche limite
948       plcl(:,nsrf)   = 0.        ! Niveau de condensation de la CLA
949       capCL(:,nsrf)  = 0.        ! CAPE de couche limite
950       oliqCL(:,nsrf) = 0.        ! eau_liqu integree de couche limite
951       cteiCL(:,nsrf) = 0.        ! cloud top instab. crit. couche limite
952       pblt(:,nsrf)   = 0.        ! T a la Hauteur de couche limite
953       therm(:,nsrf)  = 0.
954       trmb1(:,nsrf)  = 0.        ! deep_cape
955       trmb2(:,nsrf)  = 0.        ! inhibition
956       trmb3(:,nsrf)  = 0.        ! Point Omega
957
958#undef T2m     
959#define T2m     
960#ifdef T2m
961! diagnostic t,q a 2m et u, v a 10m
962
963       DO j=1, knon
964          i = ni(j)
965          uzon(j) = yu(j,1) + y_d_u(j,1)
966          vmer(j) = yv(j,1) + y_d_v(j,1)
967          tair1(j) = yt(j,1) + y_d_t(j,1)
968          qair1(j) = yq(j,1) + y_d_q(j,1)
969          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
970               * (ypaprs(j,1)-ypplay(j,1))
971          tairsol(j) = yts(j) + y_d_ts(j)
972          rugo1(j) = yrugos(j)
973          IF(nsrf.EQ.is_oce) THEN
974             rugo1(j) = rugos(i,nsrf)
975          ENDIF
976          psfce(j)=ypaprs(j,1)
977          patm(j)=ypplay(j,1)
978          qairsol(j) = yqsurf(j)
979       END DO
980       
981
982! Calculate the temperature et relative humidity at 2m and the wind at 10m
983       CALL stdlevvar(klon, knon, nsrf, zxli, &
984            uzon, vmer, tair1, qair1, zgeo1, &
985            tairsol, qairsol, rugo1, psfce, patm, &
986            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
987
988       DO j=1, knon
989          i = ni(j)
990          t2m(i,nsrf)=yt2m(j)
991         
992          q2m(i,nsrf)=yq2m(j)
993
994! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
995          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
996          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
997
998       END DO
999
1000
1001       CALL HBTM(knon, ypaprs, ypplay, &
1002            yt2m,yt10m,yq2m,yq10m,yustar, &
1003            y_flux_t,y_flux_q,yu,yv,yt,yq, &
1004            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
1005            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1006       
1007       DO j=1, knon
1008          i = ni(j)
1009          pblh(i,nsrf)   = ypblh(j)
1010          plcl(i,nsrf)   = ylcl(j)
1011          capCL(i,nsrf)  = ycapCL(j)
1012          oliqCL(i,nsrf) = yoliqCL(j)
1013          cteiCL(i,nsrf) = ycteiCL(j)
1014          pblT(i,nsrf)   = ypblT(j)
1015          therm(i,nsrf)  = ytherm(j)
1016          trmb1(i,nsrf)  = ytrmb1(j)
1017          trmb2(i,nsrf)  = ytrmb2(j)
1018          trmb3(i,nsrf)  = ytrmb3(j)
1019       END DO
1020       
1021#else
1022! not defined T2m
1023! No calculation
1024! Set output variables to zero
1025       DO j = 1, knon
1026          i = ni(j)
1027          pblh(i,nsrf)   = 0.
1028          plcl(i,nsrf)   = 0.
1029          capCL(i,nsrf)  = 0.
1030          oliqCL(i,nsrf) = 0.
1031          cteiCL(i,nsrf) = 0.
1032          pblT(i,nsrf)   = 0.
1033          therm(i,nsrf)  = 0.
1034          trmb1(i,nsrf)  = 0.
1035          trmb2(i,nsrf)  = 0.
1036          trmb3(i,nsrf)  = 0.
1037       END DO
1038       DO j = 1, knon
1039          i = ni(j)
1040          t2m(i,nsrf)=0.
1041          q2m(i,nsrf)=0.
1042          u10m(i,nsrf)=0.
1043          v10m(i,nsrf)=0.
1044       END DO
1045#endif
1046
1047!****************************************************************************************
1048! 15) End of loop over different surfaces
1049!
1050!****************************************************************************************
1051    END DO loop_nbsrf
1052
1053!****************************************************************************************
1054! 16) Calculate the mean value over all sub-surfaces for som variables
1055!
1056!     NB!!! jg : Pour garder la convergence numerique j'utilise pctsrf_new comme c'etait
1057!     fait dans physiq.F mais ca devrait plutot etre pctsrf???!!!!! A verifier!
1058!****************************************************************************************
1059   
1060    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1061    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1062    DO nsrf = 1, nbsrf
1063       DO k = 1, klev
1064          DO i = 1, klon
1065             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf_new(i,nsrf)
1066             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf_new(i,nsrf)
1067             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf_new(i,nsrf)
1068             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf_new(i,nsrf)
1069          END DO
1070       END DO
1071    END DO
1072
1073    DO i = 1, klon
1074       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1075       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
1076       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1077    ENDDO
1078   
1079
1080    DO i = 1, klon
1081       IF ( ABS( pctsrf_new(i, is_ter) + pctsrf_new(i, is_lic) + &
1082            pctsrf_new(i, is_oce) + pctsrf_new(i, is_sic)  - 1.) .GT. EPSFRA) &
1083            THEN
1084          WRITE(*,*) 'physiq : pb sous surface au point ', i, &
1085               pctsrf_new(i, 1 : nbsrf)
1086       ENDIF
1087    ENDDO
1088
1089!
1090! Incrementer la temperature du sol
1091!
1092    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
1093    zt2m(:) = 0.0    ; zq2m(:) = 0.0
1094    zu10m(:) = 0.0   ; zv10m(:) = 0.0
1095    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0
1096    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1097    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1098    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1099    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1100   
1101   
1102    DO nsrf = 1, nbsrf
1103       DO i = 1, klon         
1104          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1105         
1106          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1107               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf_new(i,nsrf)
1108          wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1109               pctsrf_new(i,nsrf)
1110
1111          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf_new(i,nsrf)
1112          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf_new(i,nsrf)
1113         
1114          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf_new(i,nsrf)
1115          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf_new(i,nsrf)
1116          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf_new(i,nsrf)
1117          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf_new(i,nsrf)
1118
1119          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf_new(i,nsrf)
1120          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf_new(i,nsrf)
1121          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf_new(i,nsrf)
1122          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf_new(i,nsrf)
1123          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf_new(i,nsrf)
1124          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf_new(i,nsrf)
1125          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf_new(i,nsrf)
1126          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf_new(i,nsrf)
1127          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf_new(i,nsrf)
1128          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf_new(i,nsrf)
1129       END DO
1130    END DO
1131
1132    IF (check) THEN
1133       amn=MIN(ts(1,is_ter),1000.)
1134       amx=MAX(ts(1,is_ter),-1000.)
1135       DO i=2, klon
1136          amn=MIN(ts(i,is_ter),amn)
1137          amx=MAX(ts(i,is_ter),amx)
1138       ENDDO
1139       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1140    ENDIF
1141!
1142! If a sub-surface does not exsist for a grid point, the mean value for all
1143! sub-surfaces is distributed.
1144!
1145    DO nsrf = 1, nbsrf
1146       DO i = 1, klon
1147          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1148             ts(i,nsrf)     = zxtsol(i)
1149             t2m(i,nsrf)    = zt2m(i)
1150             q2m(i,nsrf)    = zq2m(i)
1151             u10m(i,nsrf)   = zu10m(i)
1152             v10m(i,nsrf)   = zv10m(i)
1153
1154! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1155             pblh(i,nsrf)   = s_pblh(i)
1156             plcl(i,nsrf)   = s_plcl(i)
1157             capCL(i,nsrf)  = s_capCL(i)
1158             oliqCL(i,nsrf) = s_oliqCL(i)
1159             cteiCL(i,nsrf) = s_cteiCL(i)
1160             pblT(i,nsrf)   = s_pblT(i)
1161             therm(i,nsrf)  = s_therm(i)
1162             trmb1(i,nsrf)  = s_trmb1(i)
1163             trmb2(i,nsrf)  = s_trmb2(i)
1164             trmb3(i,nsrf)  = s_trmb3(i)
1165          ENDIF
1166       ENDDO
1167    ENDDO
1168
1169
1170    DO i = 1, klon
1171       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3
1172    ENDDO
1173   
1174    zxqsurf(:) = 0.0
1175    zxsnow(:)  = 0.0
1176    DO nsrf = 1, nbsrf
1177       DO i = 1, klon
1178          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf_new(i,nsrf)
1179          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf_new(i,nsrf)
1180       END DO
1181    END DO
1182
1183!
1184!IM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
1185!IM ajout dependance type surface
1186    rh2m(:)   = 0.0
1187    qsat2m(:) = 0.0
1188
1189    DO i = 1, klon
1190       DO nsrf=1, nbsrf
1191          zx_t1(i,nsrf) = t2m(i,nsrf)
1192          IF (thermcep) THEN
1193             zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
1194             zx_qs1(i,nsrf)  = r2es * &
1195                  FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
1196             zx_qs1(i,nsrf)  = MIN(0.5,zx_qs1(i,nsrf))
1197             zcor1(i,nsrf)   = 1./(1.-retv*zx_qs1(i,nsrf))
1198             zx_qs1(i,nsrf)  = zx_qs1(i,nsrf)*zcor1(i,nsrf)
1199          END IF
1200          zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
1201          zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
1202          rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf_new(i,nsrf)
1203          qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf_new(i,nsrf)
1204       END DO
1205    END DO
1206
1207! Some of the module declared variables are returned for printing in physiq.F
1208    qsol_d(:)     = qsol(:)
1209    evap_d(:,:)   = evap(:,:)
1210    rugos_d(:,:)  = rugos(:,:)
1211    agesno_d(:,:) = agesno(:,:)
1212
1213
1214  END SUBROUTINE pbl_surface
1215!
1216!****************************************************************************************
1217!
1218  SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1219       evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1220
1221    INCLUDE "indicesol.h"
1222    INCLUDE "dimsoil.h"
1223
1224! Ouput variables
1225!****************************************************************************************
1226    REAL, DIMENSION(klon), INTENT(OUT)                 :: qsol_rst
1227    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
1228    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
1229    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
1230    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: evap_rst
1231    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: rugos_rst
1232    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: agesno_rst
1233    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1234
1235 
1236!****************************************************************************************
1237! Return module variables for writing to restart file
1238!
1239!****************************************************************************************   
1240    qsol_rst(:)       = qsol(:)
1241    fder_rst(:)       = fder(:)
1242    snow_rst(:,:)     = snow(:,:)
1243    qsurf_rst(:,:)    = qsurf(:,:)
1244    evap_rst(:,:)     = evap(:,:)
1245    rugos_rst(:,:)    = rugos(:,:)
1246    agesno_rst(:,:)   = agesno(:,:)
1247    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1248
1249!****************************************************************************************
1250! Deallocate module variables
1251!
1252!****************************************************************************************
1253    DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1254
1255  END SUBROUTINE pbl_surface_final
1256
1257!****************************************************************************************
1258
1259
1260END MODULE pbl_surface_mod
Note: See TracBrowser for help on using the repository browser.