source: LMDZ4/trunk/libf/phy_IPCC_AR4/pbl_surface_mod.F90 @ 868

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

Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4
par la version de la physique avec thermique. On garde le repertoire phylmd
pour un petit moment pour que les utilisateurs ne soient pas trop perdus ...
phy_IPCC_AR4 = phylmd
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.0 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 )
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! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
243! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
244! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
245! dflux_t--output-R- derive du flux sensible
246! dflux_q--output-R- derive du flux latent
247! zu1------output-R- le vent dans la premiere couche
248! zv1------output-R- le vent dans la premiere couche
249! trmb1----output-R- deep_cape
250! trmb2----output-R- inhibition
251! trmb3----output-R- Point Omega
252! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
253! plcl-----output-R- Niveau de condensation
254! pblh-----output-R- HCL
255! pblT-----output-R- T au nveau HCL
256!
257    INCLUDE "indicesol.h"
258    INCLUDE "dimsoil.h"
259    INCLUDE "YOMCST.h"
260    INCLUDE "iniprint.h"
261    INCLUDE "FCTTRE.h"
262    INCLUDE "clesphys.h"
263    INCLUDE "compbl.h"
264    INCLUDE "dimensions.h"
265    INCLUDE "YOETHF.h"
266    INCLUDE "temps.h"
267    INCLUDE "control.h"
268
269! Input variables
270!****************************************************************************************
271    REAL,                         INTENT(IN)        :: dtime
272    REAL,                         INTENT(IN)        :: date0
273    INTEGER,                      INTENT(IN)        :: itap
274    INTEGER,                      INTENT(IN)        :: jour    ! jour de l'annee en cours
275    LOGICAL,                      INTENT(IN)        :: debut, lafin
276    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon, rlat
277    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro
278    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosinus de l'angle solaire zenithal
279    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f, snow_f
280    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! mean value
281    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! mean value
282    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t, q
283    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u, v
284    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay
285    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs
286    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf
287
288! Input/Output variables
289!****************************************************************************************
290    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts
291    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: albe
292    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alblw
293    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m, v10m
294
295! Output variables
296!****************************************************************************************
297    REAL, DIMENSION(klon),        INTENT(OUT)       :: sollwdown
298    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh, cdragm
299    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1
300    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1
301    REAL, DIMENSION(klon),        INTENT(OUT)       :: albsol
302    REAL, DIMENSION(klon),        INTENT(OUT)       :: albsollw
303    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens, zxevap
304    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol
305    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat
306    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m
307    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
308    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t, d_q
309    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u, d_v
310    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefh
311    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pctsrf_new
312
313! Output only for diagnostics
314    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol_d
315    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m
316    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh
317    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl
318    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL
319    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL
320    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL
321    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT
322    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm
323    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1
324    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2
325    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3
326    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxrugs
327    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m
328    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m
329    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print
330    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf
331    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m
332    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu, zxfluxv
333    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: rugos_d
334    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: agesno_d
335    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw, solsw
336    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts
337    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: evap_d
338    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat
339    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m
340    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils, wfbilo
341    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t
342    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u, flux_v
343
344! Output not needed
345    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t, dflux_q
346    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow
347    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt, zxfluxq
348    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m
349    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q
350
351
352! Local variables with attribute SAVE
353!****************************************************************************************
354    INTEGER, SAVE                            :: nhoridbg, nidbg
355!$OMP THREADPRIVATE(nhoridbg, nidbg)
356    LOGICAL, SAVE                            :: debugindex=.FALSE.
357!$OMP THREADPRIVATE(debugindex)
358    LOGICAL, SAVE                            :: first_call=.TRUE.
359!$OMP THREADPRIVATE(first_call)
360    CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
361!$OMP THREADPRIVATE(cl_surf)
362
363! Other local variables
364!****************************************************************************************
365    INTEGER                            :: i, k, nsrf
366    INTEGER                            :: knon, j
367    INTEGER                            :: idayref
368    INTEGER , DIMENSION(klon)          :: ni
369    REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
370    REAL                               :: amn, amx
371    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
372    REAL, DIMENSION(klon)              :: yts, yrugos, ypct, yz0_new
373    REAL, DIMENSION(klon)              :: yalb
374    REAL, DIMENSION(klon)              :: yalblw
375    REAL, DIMENSION(klon)              :: yu1, yv1
376    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
377    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
378    REAL, DIMENSION(klon)              :: ysollw, ysolsw, ysollwdown
379    REAL, DIMENSION(klon)              :: yfder
380    REAL, DIMENSION(klon)              :: yrads,yrugoro
381    REAL, DIMENSION(klon)              :: yfluxlat
382    REAL, DIMENSION(klon)              :: y_d_ts
383    REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
384    REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
385    REAL, DIMENSION(klon)              :: u1lay, v1lay
386    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
387    REAL, DIMENSION(klon)              :: yustar
388    REAL, DIMENSION(klon)              :: yu10mx
389    REAL, DIMENSION(klon)              :: yu10my
390    REAL, DIMENSION(klon)              :: ywindsp
391    REAL, DIMENSION(klon)              :: yt10m, yq10m
392    REAL, DIMENSION(klon)              :: ypblh
393    REAL, DIMENSION(klon)              :: ylcl
394    REAL, DIMENSION(klon)              :: ycapCL
395    REAL, DIMENSION(klon)              :: yoliqCL
396    REAL, DIMENSION(klon)              :: ycteiCL
397    REAL, DIMENSION(klon)              :: ypblT
398    REAL, DIMENSION(klon)              :: ytherm
399    REAL, DIMENSION(klon)              :: ytrmb1
400    REAL, DIMENSION(klon)              :: ytrmb2
401    REAL, DIMENSION(klon)              :: ytrmb3
402    REAL, DIMENSION(klon)              :: uzon, vmer
403    REAL, DIMENSION(klon)              :: tair1, qair1, tairsol
404    REAL, DIMENSION(klon)              :: psfce, patm
405    REAL, DIMENSION(klon)              :: qairsol, zgeo1
406    REAL, DIMENSION(klon)              :: rugo1
407    REAL, DIMENSION(klon)              :: yfluxsens, swdown
408    REAL, DIMENSION(klon)              :: petAcoef, peqAcoef, petBcoef, peqBcoef
409    REAL, DIMENSION(klon)              :: ypsref, epot_air
410    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb_new
411    REAL, DIMENSION(klon)              :: pctsrf_nsrf
412    REAL, DIMENSION(klon)              :: ztsol
413    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q
414    REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
415    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
416    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
417    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm
418    REAL, DIMENSION(klon,klev)         :: yu, yv
419    REAL, DIMENSION(klon,klev)         :: yt, yq
420    REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
421    REAL, DIMENSION(klon,klev)         :: delp
422    REAL, DIMENSION(klon,klev+1)       :: ypaprs
423    REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
424    REAL, DIMENSION(klon,nbsrf)        :: pctsrf_pot
425    CHARACTER(len=80)                  :: abort_message
426    CHARACTER(len=20)                  :: modname = 'pbl_surface'
427    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
428    LOGICAL, PARAMETER                 :: check=.FALSE.
429
430! For debugging with IOIPSL
431    INTEGER, DIMENSION(iim*(jjm+1))    :: ndexbg
432    REAL                               :: zjulian
433    REAL, DIMENSION(klon)              :: tabindx
434    REAL, DIMENSION(iim,jjm+1)         :: zx_lon, zx_lat
435    REAL, DIMENSION(iim,jjm+1)         :: debugtab
436
437
438    REAL, DIMENSION(klon,nbsrf)        :: pblh
439    REAL, DIMENSION(klon,nbsrf)        :: plcl
440    REAL, DIMENSION(klon,nbsrf)        :: capCL
441    REAL, DIMENSION(klon,nbsrf)        :: oliqCL
442    REAL, DIMENSION(klon,nbsrf)        :: cteiCL
443    REAL, DIMENSION(klon,nbsrf)        :: pblT
444    REAL, DIMENSION(klon,nbsrf)        :: therm
445    REAL, DIMENSION(klon,nbsrf)        :: trmb1
446    REAL, DIMENSION(klon,nbsrf)        :: trmb2
447    REAL, DIMENSION(klon,nbsrf)        :: trmb3
448    REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
449    REAL, DIMENSION(klon,nbsrf)        :: zx_qs1, zx_t1
450    REAL, DIMENSION(klon,nbsrf)        :: zdelta1, zcor1
451
452
453!****************************************************************************************
454! End of declarations
455!****************************************************************************************
456
457
458!****************************************************************************************
459! 1) Initialisation and validation tests
460!    Only done first time entering this subroutine
461!
462!****************************************************************************************
463
464
465    IF (first_call) THEN
466       first_call=.FALSE.
467     
468       ! Initilize debug IO
469       IF (debugindex .AND. mpi_size==1) THEN
470          ! initialize IOIPSL output
471          idayref = day_ini
472          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
473          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
474          DO i = 1, iim
475             zx_lon(i,1) = rlon(i+1)
476             zx_lon(i,jjm+1) = rlon(i+1)
477          ENDDO
478          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
479          CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
480               1,iim,1,jjm+1, &
481               itau_phy,zjulian,dtime,nhoridbg,nidbg)
482          ! no vertical axis
483          cl_surf(1)='ter'
484          cl_surf(2)='lic'
485          cl_surf(3)='oce'
486          cl_surf(4)='sic'
487          DO nsrf=1,nbsrf
488             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
489                  jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
490          END DO
491
492          CALL histend(nidbg)
493          CALL histsync(nidbg)
494
495       END IF
496       
497    ENDIF
498         
499!****************************************************************************************
500! 2) Initialization to zero
501!    Done for all local variables that will be compressed later
502!    and argument with INTENT(OUT)
503!****************************************************************************************
504    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
505    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0     ; zu1 = 0.0       
506    zv1 = 0.0     ; yqsurf = 0.0     ; yalb = 0.0      ; yalblw = 0.0   
507    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0   
508    ysollw = 0.0  ; ysollwdown = 0.0 ; yrugos = 0.0    ; yu1 = 0.0   
509    yv1 = 0.0     ; yrads = 0.0      ; ypaprs = 0.0    ; ypplay = 0.0
510    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
511    yq = 0.0      ; pctsrf_new = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
512    yrugoro = 0.0 ; yu10mx = 0.0     ; yu10my = 0.0    ; ywindsp = 0.0   
513    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
514    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
515    d_u = 0.0     ; d_v = 0.0        ; zcoefh = 0.0    ; yqsol = 0.0   
516    ytherm = 0.0
517     
518    ytsoil = 999999.
519
520!****************************************************************************************
521! 3) - Calculate pressure thickness of each layer
522!    - Calculate the wind at first layer
523!
524!****************************************************************************************
525    DO k = 1, klev
526       DO i = 1, klon
527          delp(i,k) = paprs(i,k)-paprs(i,k+1)
528       ENDDO
529    ENDDO
530    DO i = 1, klon
531       zx_alf1 = 1.0
532       zx_alf2 = 1.0 - zx_alf1
533       u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
534       v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
535    ENDDO
536
537
538!****************************************************************************************
539! Test for rugos........ from physiq.. A la fin plutot???
540! Calcul de l'abedo moyen par maille
541!****************************************************************************************
542
543    zxrugs(:) = 0.0
544    DO nsrf = 1, nbsrf
545       DO i = 1, klon
546          rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
547          zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
548       ENDDO
549    ENDDO
550
551! Calcul de l'abedo moyen par maille
552    albsol(:)   = 0.0
553    albsollw(:) = 0.0
554    DO nsrf = 1, nbsrf
555       DO i = 1, klon
556          albsol(i)   = albsol(i)   + albe(i,nsrf)  * pctsrf(i,nsrf)
557          albsollw(i) = albsollw(i) + alblw(i,nsrf) * pctsrf(i,nsrf)
558       ENDDO
559    ENDDO
560
561
562
563! Calcule de ztsol (aussi fait dans physiq.F, pourrait etre un argument)
564    ztsol(:) = 0.0
565    DO nsrf = 1, nbsrf
566       DO i = 1, klon
567          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
568       ENDDO
569    ENDDO
570
571
572! Repartition du longwave par sous-surface linearisee
573    DO nsrf = 1, nbsrf
574       DO i = 1, klon
575          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
576          solsw(i,nsrf) = solsw_m(i)*(1.-albe(i,nsrf))/(1.-albsol(i))
577       ENDDO
578    ENDDO
579
580
581    DO i = 1, klon
582       sollwdown(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
583    ENDDO
584
585!****************************************************************************************
586! 4) Loop over different surfaces
587!
588! All points with a possibility of the current surface are used. This is done
589! to allow the sea-ice to appear or disappear. It is considered here that the
590! entier domaine of the ocean possibly can contain sea-ice.
591!
592!****************************************************************************************
593
594    pctsrf_pot = pctsrf
595    pctsrf_pot(:,is_oce) = 1. - zmasq(:)
596    pctsrf_pot(:,is_sic) = 1. - zmasq(:)
597     
598    loop_nbsrf: DO nsrf = 1, nbsrf
599
600! Search for index(ni) and size(knon) of domaine to treat
601       ni(:) = 0
602       knon  = 0
603       DO i = 1, klon
604          IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
605             knon = knon + 1
606             ni(knon) = i
607          ENDIF
608       ENDDO
609
610       ! write index, with IOIPSL
611       IF (debugindex .AND. mpi_size==1) THEN
612          tabindx(:)=0.
613          DO i=1,knon
614             tabindx(i)=FLOAT(i)
615          END DO
616          debugtab(:,:) = 0.
617          ndexbg(:) = 0
618          CALL gath2cpl(tabindx,debugtab,knon,ni)
619          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
620       ENDIF
621       
622!****************************************************************************************
623! 5) Compress variables
624!
625!****************************************************************************************
626
627       DO j = 1, knon
628          i = ni(j)
629          ypct(j) = pctsrf(i,nsrf)
630          yts(j) = ts(i,nsrf)
631          ysnow(j) = snow(i,nsrf)
632          yqsurf(j) = qsurf(i,nsrf)
633          yalb(j) = albe(i,nsrf)
634          yalblw(j) = alblw(i,nsrf)
635          yrain_f(j) = rain_f(i)
636          ysnow_f(j) = snow_f(i)
637          yagesno(j) = agesno(i,nsrf)
638          yfder(j) = fder(i)
639          ysolsw(j) = solsw(i,nsrf)
640          ysollw(j) = sollw(i,nsrf)
641          ysollwdown(j) = sollwdown(i)
642          yrugos(j) = rugos(i,nsrf)
643          yrugoro(j) = rugoro(i)
644          yu1(j) = u1lay(i)
645          yv1(j) = v1lay(i)
646          yrads(j) =  ysolsw(j)+ ysollw(j)
647          ypaprs(j,klev+1) = paprs(i,klev+1)
648          yu10mx(j) = u10m(i,nsrf)
649          yu10my(j) = v10m(i,nsrf)
650          ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )
651       END DO
652
653       DO k = 1, klev
654          DO j = 1, knon
655             i = ni(j)
656             ypaprs(j,k) = paprs(i,k)
657             ypplay(j,k) = pplay(i,k)
658             ydelp(j,k) = delp(i,k)
659             yu(j,k) = u(i,k)
660             yv(j,k) = v(i,k)
661             yt(j,k) = t(i,k)
662             yq(j,k) = q(i,k)
663          ENDDO
664       ENDDO
665       
666       DO k = 1, nsoilmx
667          DO j = 1, knon
668             i = ni(j)
669             ytsoil(j,k) = ftsoil(i,k,nsrf)
670          END DO
671       END DO
672       
673       ! qsol(water height in soil) only for bucket continental model
674       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
675          DO j = 1, knon
676             i = ni(j)
677             yqsol(j) = qsol(i)
678          END DO
679       ENDIF
680       
681!****************************************************************************************
682! 6) Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
683!    atmosphere and coefficients for turbulent diffusion at surface(Cdrag).
684!
685!****************************************************************************************
686
687       CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
688            ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf,  &
689            ycoefm, ycoefh)
690       
691!****************************************************************************************
692!
693! 8) "La descente" - "The downhill"
694
695!  climb_hq_down and climb_wind_down calculate the coefficients
696!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
697!  Only the coefficients at surface for H and Q are returned.
698!
699!****************************************************************************************
700
701! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
702       CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
703            ydelp, yt, yq, dtime, &
704            petAcoef, peqAcoef, petBcoef, peqBcoef)
705
706! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
707       CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv)
708     
709
710!****************************************************************************************
711! 9) Small calculations
712!
713!****************************************************************************************
714         
715       ypsref(:) = ypaprs(:,1) 
716       epot_air(:) = 0.0
717       epot_air(1:knon) = RCPD*yt(1:knon,1)*(ypsref(1:knon)/ypplay(1:knon,1))**RKAPPA
718
719       swdown(:) = 0.0
720       IF (nsrf .EQ. is_ter) THEN
721          swdown(1:knon) = ysolsw(1:knon)/(1-yalb(1:knon))
722       ELSE
723          swdown(1:knon) = ysolsw(1:knon)
724       ENDIF
725
726       ! constant CO2
727       r_co2_ppm(:) = co2_ppm
728
729!****************************************************************************************
730!
731! 10) Switch selon current surface
732!     It is necessary to start with the continental surfaces because the ocean
733!     needs their run-off.
734!
735!****************************************************************************************
736       SELECT CASE(nsrf)
737     
738       CASE(is_ter)
739          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
740               rlon, rlat, &
741               debut, lafin, ydelp(:,1), epot_air, r_co2_ppm, ysollwdown, ysolsw, swdown, &
742               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
743               petAcoef, peqAcoef, petBcoef, peqBcoef, &
744               ypsref, yu1, yv1, yrugoro, pctsrf, &
745               yrads, ysnow, yqsurf, yqsol, yagesno, &
746               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
747               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
748     
749       CASE(is_lic)
750          CALL surf_landice(itap, dtime, knon, ni, &
751               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
752               petAcoef, peqAcoef, petBcoef, peqBcoef, &
753               ypsref, yu1, yv1, yrugoro, pctsrf, &
754               yrads, ysnow, yqsurf, yqsol, yagesno, &
755               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
756               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
757         
758       CASE(is_oce)
759          CALL surf_ocean(rlon, rlat, ysollw, yalb, &
760               yrugos, ywindsp, rmu0, &
761               yfder, &
762               itap, dtime, jour, knon, ni, &
763               debut, swdown, &
764               ypplay(:,1), ycoefh(:,1), ycoefm(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
765               petAcoef, peqAcoef, petBcoef, peqBcoef, &
766               ypsref, yu1, yv1, yrugoro, pctsrf, &
767               yrads, ysnow, yqsurf, yagesno, &
768               yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
769               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
770         
771       CASE(is_sic)
772          CALL surf_seaice( &
773               rlon, rlat, ysollw, yalb, &
774               yfder, &
775               itap, dtime, jour, knon, ni, &
776               debut, lafin, swdown, &
777               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
778               petAcoef, peqAcoef, petBcoef, peqBcoef, &
779               ypsref, yu1, yv1, yrugoro, pctsrf, &
780               yrads, ysnow, yqsurf, yqsol, yagesno, &
781               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
782               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
783         
784
785       CASE DEFAULT
786          WRITE(lunout,*) 'Surface index = ', nsrf
787          abort_message = 'Surface index not valid'
788          CALL abort_gcm(modname,abort_message,1)
789       END SELECT
790
791!****************************************************************************************
792! Save the fraction of this sub-surface
793!
794!****************************************************************************************
795       pctsrf_new(:,nsrf) = pctsrf_nsrf(:)
796
797!****************************************************************************************
798! 11) - Calcul the increment of surface temperature
799!     - Update albedo
800!
801!****************************************************************************************
802       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
803 
804       yalb(1:knon) = yalb_new(1:knon)
805
806!****************************************************************************************
807!
808! 12) "La remontee" - "The uphill"
809!
810!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
811!  for X=H, Q, U and V, for all vertical levels.
812!
813!****************************************************************************************
814! H and Q
815       y_flux_t1(:) =  yfluxsens(:)
816       y_flux_q1(:) = -yevap(:)
817
818       CALL climb_hq_up(knon, dtime, yt, yq, &
819            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
820            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))   
821       
822! U and V
823       CALL climb_wind_up(knon, dtime, yu, yv, &
824            y_flux_u, y_flux_v, y_d_u, y_d_v)
825
826       DO j = 1, knon
827          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
828          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
829          yu1(j) = yu1(j) *  ypct(j)
830          yv1(j) = yv1(j) *  ypct(j)
831       ENDDO
832
833!****************************************************************************************
834! 13) Transform variables for output format :
835!     - Decompress
836!     - Multiply with pourcentage of current surface
837!     - Cumulate in global variable
838!
839!****************************************************************************************
840
841       DO k = 1, klev
842          DO j = 1, knon
843             i = ni(j)
844             ycoefh(j,k) = ycoefh(j,k) * ypct(j)
845             ycoefm(j,k) = ycoefm(j,k) * ypct(j)
846             y_d_t(j,k) = y_d_t(j,k) * ypct(j)
847             y_d_q(j,k) = y_d_q(j,k) * ypct(j)
848             y_d_u(j,k) = y_d_u(j,k) * ypct(j)
849             y_d_v(j,k) = y_d_v(j,k) * ypct(j)
850
851             flux_t(i,k,nsrf) = y_flux_t(j,k)
852             flux_q(i,k,nsrf) = y_flux_q(j,k)
853             flux_u(i,k,nsrf) = y_flux_u(j,k)
854             flux_v(i,k,nsrf) = y_flux_v(j,k)
855          ENDDO
856       ENDDO
857       
858       evap(:,nsrf) = - flux_q(:,1,nsrf)
859       
860       albe(:, nsrf) = 0.
861       alblw(:, nsrf) = 0.
862       snow(:, nsrf) = 0.
863       qsurf(:, nsrf) = 0.
864       rugos(:, nsrf) = 0.
865       fluxlat(:,nsrf) = 0.
866       DO j = 1, knon
867          i = ni(j)
868          d_ts(i,nsrf) = y_d_ts(j)
869          albe(i,nsrf) = yalb(j) 
870          alblw(i,nsrf) = yalblw(j)
871          snow(i,nsrf) = ysnow(j) 
872          qsurf(i,nsrf) = yqsurf(j)
873          rugos(i,nsrf) = yz0_new(j)
874          fluxlat(i,nsrf) = yfluxlat(j)
875          agesno(i,nsrf) = yagesno(j) 
876          cdragh(i) = cdragh(i) + ycoefh(j,1)
877          cdragm(i) = cdragm(i) + ycoefm(j,1)
878          dflux_t(i) = dflux_t(i) + y_dflux_t(j)
879          dflux_q(i) = dflux_q(i) + y_dflux_q(j)
880          zu1(i) = zu1(i) + yu1(j)
881          zv1(i) = zv1(i) + yv1(j)
882       END DO
883
884       IF ( nsrf .EQ. is_ter ) THEN
885          DO j = 1, knon
886             i = ni(j)
887             qsol(i) = yqsol(j)
888          END DO
889       END IF
890       
891       ftsoil(:,:,nsrf) = 0.
892       DO k = 1, nsoilmx
893          DO j = 1, knon
894             i = ni(j)
895             ftsoil(i, k, nsrf) = ytsoil(j,k)
896          END DO
897       END DO
898       
899       
900#ifdef CRAY
901       DO k = 1, klev
902          DO j = 1, knon
903             i = ni(j)
904#else
905       DO j = 1, knon
906          i = ni(j)
907          DO k = 1, klev
908#endif
909             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
910             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
911             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
912             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
913             zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
914#ifdef CRAY
915          END DO
916       END DO
917#else
918          END DO
919       END DO
920#endif
921
922!****************************************************************************************
923! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
924!     Call HBTM
925!
926!****************************************************************************************
927       t2m(:,nsrf)    = 0.
928       q2m(:,nsrf)    = 0.
929       u10m(:,nsrf)   = 0.
930       v10m(:,nsrf)   = 0.
931
932       pblh(:,nsrf)   = 0.        ! Hauteur de couche limite
933       plcl(:,nsrf)   = 0.        ! Niveau de condensation de la CLA
934       capCL(:,nsrf)  = 0.        ! CAPE de couche limite
935       oliqCL(:,nsrf) = 0.        ! eau_liqu integree de couche limite
936       cteiCL(:,nsrf) = 0.        ! cloud top instab. crit. couche limite
937       pblt(:,nsrf)   = 0.        ! T a la Hauteur de couche limite
938       therm(:,nsrf)  = 0.
939       trmb1(:,nsrf)  = 0.        ! deep_cape
940       trmb2(:,nsrf)  = 0.        ! inhibition
941       trmb3(:,nsrf)  = 0.        ! Point Omega
942
943#undef T2m     
944#define T2m     
945#ifdef T2m
946! diagnostic t,q a 2m et u, v a 10m
947
948       DO j=1, knon
949          i = ni(j)
950          uzon(j) = yu(j,1) + y_d_u(j,1)
951          vmer(j) = yv(j,1) + y_d_v(j,1)
952          tair1(j) = yt(j,1) + y_d_t(j,1)
953          qair1(j) = yq(j,1) + y_d_q(j,1)
954          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
955               * (ypaprs(j,1)-ypplay(j,1))
956          tairsol(j) = yts(j) + y_d_ts(j)
957          rugo1(j) = yrugos(j)
958          IF(nsrf.EQ.is_oce) THEN
959             rugo1(j) = rugos(i,nsrf)
960          ENDIF
961          psfce(j)=ypaprs(j,1)
962          patm(j)=ypplay(j,1)
963          qairsol(j) = yqsurf(j)
964       END DO
965       
966
967! Calculate the temperature et relative humidity at 2m and the wind at 10m
968       CALL stdlevvar(klon, knon, nsrf, zxli, &
969            uzon, vmer, tair1, qair1, zgeo1, &
970            tairsol, qairsol, rugo1, psfce, patm, &
971            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
972
973       DO j=1, knon
974          i = ni(j)
975          t2m(i,nsrf)=yt2m(j)
976         
977          q2m(i,nsrf)=yq2m(j)
978
979! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
980          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
981          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
982
983       END DO
984
985
986       CALL HBTM(knon, ypaprs, ypplay, &
987            yt2m,yt10m,yq2m,yq10m,yustar, &
988            y_flux_t,y_flux_q,yu,yv,yt,yq, &
989            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
990            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
991       
992       DO j=1, knon
993          i = ni(j)
994          pblh(i,nsrf)   = ypblh(j)
995          plcl(i,nsrf)   = ylcl(j)
996          capCL(i,nsrf)  = ycapCL(j)
997          oliqCL(i,nsrf) = yoliqCL(j)
998          cteiCL(i,nsrf) = ycteiCL(j)
999          pblT(i,nsrf)   = ypblT(j)
1000          therm(i,nsrf)  = ytherm(j)
1001          trmb1(i,nsrf)  = ytrmb1(j)
1002          trmb2(i,nsrf)  = ytrmb2(j)
1003          trmb3(i,nsrf)  = ytrmb3(j)
1004       END DO
1005       
1006#else
1007! not defined T2m
1008! No calculation
1009! Set output variables to zero
1010       DO j = 1, knon
1011          i = ni(j)
1012          pblh(i,nsrf)   = 0.
1013          plcl(i,nsrf)   = 0.
1014          capCL(i,nsrf)  = 0.
1015          oliqCL(i,nsrf) = 0.
1016          cteiCL(i,nsrf) = 0.
1017          pblT(i,nsrf)   = 0.
1018          therm(i,nsrf)  = 0.
1019          trmb1(i,nsrf)  = 0.
1020          trmb2(i,nsrf)  = 0.
1021          trmb3(i,nsrf)  = 0.
1022       END DO
1023       DO j = 1, knon
1024          i = ni(j)
1025          t2m(i,nsrf)=0.
1026          q2m(i,nsrf)=0.
1027          u10m(i,nsrf)=0.
1028          v10m(i,nsrf)=0.
1029       END DO
1030#endif
1031
1032!****************************************************************************************
1033! 15) End of loop over different surfaces
1034!
1035!****************************************************************************************
1036    END DO loop_nbsrf
1037
1038!****************************************************************************************
1039! 16) Calculate the mean value over all sub-surfaces for som variables
1040!
1041!     NB!!! jg : Pour garder la convergence numerique j'utilise pctsrf_new comme c'etait
1042!     fait dans physiq.F mais ca devrait plutot etre pctsrf???!!!!! A verifier!
1043!****************************************************************************************
1044   
1045    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1046    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1047    DO nsrf = 1, nbsrf
1048       DO k = 1, klev
1049          DO i = 1, klon
1050             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf_new(i,nsrf)
1051             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf_new(i,nsrf)
1052             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf_new(i,nsrf)
1053             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf_new(i,nsrf)
1054          END DO
1055       END DO
1056    END DO
1057
1058    DO i = 1, klon
1059       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1060       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
1061       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1062    ENDDO
1063   
1064
1065    DO i = 1, klon
1066       IF ( ABS( pctsrf_new(i, is_ter) + pctsrf_new(i, is_lic) + &
1067            pctsrf_new(i, is_oce) + pctsrf_new(i, is_sic)  - 1.) .GT. EPSFRA) &
1068            THEN
1069          WRITE(*,*) 'physiq : pb sous surface au point ', i, &
1070               pctsrf_new(i, 1 : nbsrf)
1071       ENDIF
1072    ENDDO
1073
1074!
1075! Incrementer la temperature du sol
1076!
1077    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
1078    zt2m(:) = 0.0    ; zq2m(:) = 0.0
1079    zu10m(:) = 0.0   ; zv10m(:) = 0.0
1080    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0
1081    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1082    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1083    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1084    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1085   
1086   
1087    DO nsrf = 1, nbsrf
1088       DO i = 1, klon         
1089          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1090         
1091          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1092               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf_new(i,nsrf)
1093          wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1094               pctsrf_new(i,nsrf)
1095
1096          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf_new(i,nsrf)
1097          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf_new(i,nsrf)
1098         
1099          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf_new(i,nsrf)
1100          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf_new(i,nsrf)
1101          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf_new(i,nsrf)
1102          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf_new(i,nsrf)
1103
1104          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf_new(i,nsrf)
1105          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf_new(i,nsrf)
1106          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf_new(i,nsrf)
1107          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf_new(i,nsrf)
1108          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf_new(i,nsrf)
1109          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf_new(i,nsrf)
1110          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf_new(i,nsrf)
1111          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf_new(i,nsrf)
1112          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf_new(i,nsrf)
1113          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf_new(i,nsrf)
1114       END DO
1115    END DO
1116
1117    IF (check) THEN
1118       amn=MIN(ts(1,is_ter),1000.)
1119       amx=MAX(ts(1,is_ter),-1000.)
1120       DO i=2, klon
1121          amn=MIN(ts(i,is_ter),amn)
1122          amx=MAX(ts(i,is_ter),amx)
1123       ENDDO
1124       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1125    ENDIF
1126!
1127! If a sub-surface does not exsist for a grid point, the mean value for all
1128! sub-surfaces is distributed.
1129!
1130    DO nsrf = 1, nbsrf
1131       DO i = 1, klon
1132          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1133             ts(i,nsrf)     = zxtsol(i)
1134             t2m(i,nsrf)    = zt2m(i)
1135             q2m(i,nsrf)    = zq2m(i)
1136             u10m(i,nsrf)   = zu10m(i)
1137             v10m(i,nsrf)   = zv10m(i)
1138
1139! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1140             pblh(i,nsrf)   = s_pblh(i)
1141             plcl(i,nsrf)   = s_plcl(i)
1142             capCL(i,nsrf)  = s_capCL(i)
1143             oliqCL(i,nsrf) = s_oliqCL(i)
1144             cteiCL(i,nsrf) = s_cteiCL(i)
1145             pblT(i,nsrf)   = s_pblT(i)
1146             therm(i,nsrf)  = s_therm(i)
1147             trmb1(i,nsrf)  = s_trmb1(i)
1148             trmb2(i,nsrf)  = s_trmb2(i)
1149             trmb3(i,nsrf)  = s_trmb3(i)
1150          ENDIF
1151       ENDDO
1152    ENDDO
1153
1154
1155    DO i = 1, klon
1156       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3
1157    ENDDO
1158   
1159    zxqsurf(:) = 0.0
1160    zxsnow(:)  = 0.0
1161    DO nsrf = 1, nbsrf
1162       DO i = 1, klon
1163          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf_new(i,nsrf)
1164          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf_new(i,nsrf)
1165       END DO
1166    END DO
1167
1168!
1169!IM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
1170!IM ajout dependance type surface
1171    rh2m(:)   = 0.0
1172    qsat2m(:) = 0.0
1173
1174    DO i = 1, klon
1175       DO nsrf=1, nbsrf
1176          zx_t1(i,nsrf) = t2m(i,nsrf)
1177          IF (thermcep) THEN
1178             zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
1179             zx_qs1(i,nsrf)  = r2es * &
1180                  FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
1181             zx_qs1(i,nsrf)  = MIN(0.5,zx_qs1(i,nsrf))
1182             zcor1(i,nsrf)   = 1./(1.-retv*zx_qs1(i,nsrf))
1183             zx_qs1(i,nsrf)  = zx_qs1(i,nsrf)*zcor1(i,nsrf)
1184          END IF
1185          zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
1186          zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
1187          rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf_new(i,nsrf)
1188          qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf_new(i,nsrf)
1189       END DO
1190    END DO
1191
1192! Some of the module declared variables are returned for printing in physiq.F
1193    qsol_d(:)     = qsol(:)
1194    evap_d(:,:)   = evap(:,:)
1195    rugos_d(:,:)  = rugos(:,:)
1196    agesno_d(:,:) = agesno(:,:)
1197
1198
1199  END SUBROUTINE pbl_surface
1200!
1201!****************************************************************************************
1202!
1203  SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1204       evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1205
1206    INCLUDE "indicesol.h"
1207    INCLUDE "dimsoil.h"
1208
1209! Ouput variables
1210!****************************************************************************************
1211    REAL, DIMENSION(klon), INTENT(OUT)                 :: qsol_rst
1212    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
1213    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
1214    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
1215    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: evap_rst
1216    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: rugos_rst
1217    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: agesno_rst
1218    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1219
1220 
1221!****************************************************************************************
1222! Return module variables for writing to restart file
1223!
1224!****************************************************************************************   
1225    qsol_rst(:)       = qsol(:)
1226    fder_rst(:)       = fder(:)
1227    snow_rst(:,:)     = snow(:,:)
1228    qsurf_rst(:,:)    = qsurf(:,:)
1229    evap_rst(:,:)     = evap(:,:)
1230    rugos_rst(:,:)    = rugos(:,:)
1231    agesno_rst(:,:)   = agesno(:,:)
1232    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1233
1234!****************************************************************************************
1235! Deallocate module variables
1236!
1237!****************************************************************************************
1238    DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1239
1240  END SUBROUTINE pbl_surface_final
1241
1242!****************************************************************************************
1243
1244
1245END MODULE pbl_surface_mod
Note: See TracBrowser for help on using the repository browser.