Context Navigation

pbl_surface_mod.F90 @ 850

Last change on this file since 850 was 803, checked in by Laurent Fairhead, 17 years ago
On elimine l'ancien calcul de clvent JG 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	!
4	MODULE 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
45	CONTAINS
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/m2/s (W/m2)
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.0RSIGMAztsol(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) + RSIGMAztsol(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) = RCPDyt(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.0RSIGMAzxtsol(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
1245	END MODULE pbl_surface_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format