Context Navigation

pbl_surface_mod.F90 @ 789

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

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

Download in other formats:

Original Format