Context Navigation

pbl_surface_mod.F90 @ 862

Last change on this file since 862 was 852, checked in by Laurent Fairhead, 17 years ago
Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007 LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 48.4 KB

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

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

Download in other formats:

Original Format