Context Navigation

pbl_surface_mod.F90 @ 882

Last change on this file since 882 was 882, checked in by Laurent Fairhead, 17 years ago
Modifs pour intégrer le 1D FH 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.9 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	! Declarations specifiques pour le 1D. A reprendre
465	REAL :: fsens,flat
466	LOGICAL ok_flux_surf
467	data ok_flux_surf/.false./
468	common /flux_arp/fsens,flat,ok_flux_surf
469
470	!****************************************************************************************
471	! End of declarations
472	!****************************************************************************************
473
474
475	!****************************************************************************************
476	! 1) Initialisation and validation tests
477	! Only done first time entering this subroutine
478	!
479	!****************************************************************************************
480
481
482	IF (first_call) THEN
483	first_call=.FALSE.
484
485	! Initilize debug IO
486	IF (debugindex .AND. mpi_size==1) THEN
487	! initialize IOIPSL output
488	idayref = day_ini
489	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
490	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
491	DO i = 1, iim
492	zx_lon(i,1) = rlon(i+1)
493	zx_lon(i,jjm+1) = rlon(i+1)
494	ENDDO
495	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
496	CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
497	1,iim,1,jjm+1, &
498	itau_phy,zjulian,dtime,nhoridbg,nidbg)
499	! no vertical axis
500	cl_surf(1)='ter'
501	cl_surf(2)='lic'
502	cl_surf(3)='oce'
503	cl_surf(4)='sic'
504	DO nsrf=1,nbsrf
505	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
506	jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
507	END DO
508
509	CALL histend(nidbg)
510	CALL histsync(nidbg)
511
512	END IF
513
514	ENDIF
515
516	!****************************************************************************************
517	! 2) Initialization to zero
518	! Done for all local variables that will be compressed later
519	! and argument with INTENT(OUT)
520	!****************************************************************************************
521	cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
522	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0 ; zu1 = 0.0
523	zv1 = 0.0 ; yqsurf = 0.0 ; yalb = 0.0 ; yalblw = 0.0
524	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
525	ysollw = 0.0 ; ysollwdown = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
526	yv1 = 0.0 ; yrads = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
527	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
528	yq = 0.0 ; pctsrf_new = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
529	yrugoro = 0.0 ; yu10mx = 0.0 ; yu10my = 0.0 ; ywindsp = 0.0
530	d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
531	flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
532	d_u = 0.0 ; d_v = 0.0 ; zcoefh = 0.0 ; yqsol = 0.0
533	ytherm = 0.0 ; ytke=0.
534
535	ytsoil = 999999.
536
537	!****************************************************************************************
538	! 3) - Calculate pressure thickness of each layer
539	! - Calculate the wind at first layer
540	!
541	!****************************************************************************************
542	DO k = 1, klev
543	DO i = 1, klon
544	delp(i,k) = paprs(i,k)-paprs(i,k+1)
545	ENDDO
546	ENDDO
547	DO i = 1, klon
548	zx_alf1 = 1.0
549	zx_alf2 = 1.0 - zx_alf1
550	u1lay(i) = u(i,1)zx_alf1 + u(i,2)zx_alf2
551	v1lay(i) = v(i,1)zx_alf1 + v(i,2)zx_alf2
552	ENDDO
553
554
555	!****************************************************************************************
556	! Test for rugos........ from physiq.. A la fin plutot???
557	! Calcul de l'abedo moyen par maille
558	!****************************************************************************************
559
560	zxrugs(:) = 0.0
561	DO nsrf = 1, nbsrf
562	DO i = 1, klon
563	rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
564	zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
565	ENDDO
566	ENDDO
567
568	! Calcul de l'abedo moyen par maille
569	albsol(:) = 0.0
570	albsollw(:) = 0.0
571	DO nsrf = 1, nbsrf
572	DO i = 1, klon
573	albsol(i) = albsol(i) + albe(i,nsrf) * pctsrf(i,nsrf)
574	albsollw(i) = albsollw(i) + alblw(i,nsrf) * pctsrf(i,nsrf)
575	ENDDO
576	ENDDO
577
578
579
580	! Calcule de ztsol (aussi fait dans physiq.F, pourrait etre un argument)
581	ztsol(:) = 0.0
582	DO nsrf = 1, nbsrf
583	DO i = 1, klon
584	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
585	ENDDO
586	ENDDO
587
588
589	! Repartition du longwave par sous-surface linearisee
590	DO nsrf = 1, nbsrf
591	DO i = 1, klon
592	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
593	solsw(i,nsrf) = solsw_m(i)*(1.-albe(i,nsrf))/(1.-albsol(i))
594	ENDDO
595	ENDDO
596
597
598	DO i = 1, klon
599	sollwdown(i) = sollw_m(i) + RSIGMAztsol(i)*4
600	ENDDO
601
602	!****************************************************************************************
603	! 4) Loop over different surfaces
604	!
605	! All points with a possibility of the current surface are used. This is done
606	! to allow the sea-ice to appear or disappear. It is considered here that the
607	! entier domaine of the ocean possibly can contain sea-ice.
608	!
609	!****************************************************************************************
610
611	pctsrf_pot = pctsrf
612	pctsrf_pot(:,is_oce) = 1. - zmasq(:)
613	pctsrf_pot(:,is_sic) = 1. - zmasq(:)
614
615	loop_nbsrf: DO nsrf = 1, nbsrf
616
617	! Search for index(ni) and size(knon) of domaine to treat
618	ni(:) = 0
619	knon = 0
620	DO i = 1, klon
621	IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
622	knon = knon + 1
623	ni(knon) = i
624	ENDIF
625	ENDDO
626
627	! write index, with IOIPSL
628	IF (debugindex .AND. mpi_size==1) THEN
629	tabindx(:)=0.
630	DO i=1,knon
631	tabindx(i)=FLOAT(i)
632	END DO
633	debugtab(:,:) = 0.
634	ndexbg(:) = 0
635	CALL gath2cpl(tabindx,debugtab,knon,ni)
636	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
637	ENDIF
638
639	!****************************************************************************************
640	! 5) Compress variables
641	!
642	!****************************************************************************************
643
644	DO j = 1, knon
645	i = ni(j)
646	ypct(j) = pctsrf(i,nsrf)
647	yts(j) = ts(i,nsrf)
648	ysnow(j) = snow(i,nsrf)
649	yqsurf(j) = qsurf(i,nsrf)
650	yalb(j) = albe(i,nsrf)
651	yalblw(j) = alblw(i,nsrf)
652	yrain_f(j) = rain_f(i)
653	ysnow_f(j) = snow_f(i)
654	yagesno(j) = agesno(i,nsrf)
655	yfder(j) = fder(i)
656	ysolsw(j) = solsw(i,nsrf)
657	ysollw(j) = sollw(i,nsrf)
658	ysollwdown(j) = sollwdown(i)
659	yrugos(j) = rugos(i,nsrf)
660	yrugoro(j) = rugoro(i)
661	yu1(j) = u1lay(i)
662	yv1(j) = v1lay(i)
663	yrads(j) = ysolsw(j)+ ysollw(j)
664	ypaprs(j,klev+1) = paprs(i,klev+1)
665	yu10mx(j) = u10m(i,nsrf)
666	yu10my(j) = v10m(i,nsrf)
667	ywindsp(j) = SQRT(yu10mx(j)yu10mx(j) + yu10my(j)yu10my(j) )
668	END DO
669
670	DO k = 1, klev
671	DO j = 1, knon
672	i = ni(j)
673	ypaprs(j,k) = paprs(i,k)
674	ypplay(j,k) = pplay(i,k)
675	ydelp(j,k) = delp(i,k)
676	ytke(j,k)=tke(i,k,nsrf)
677	yu(j,k) = u(i,k)
678	yv(j,k) = v(i,k)
679	yt(j,k) = t(i,k)
680	yq(j,k) = q(i,k)
681	ENDDO
682	ENDDO
683
684	DO k = 1, nsoilmx
685	DO j = 1, knon
686	i = ni(j)
687	ytsoil(j,k) = ftsoil(i,k,nsrf)
688	END DO
689	END DO
690
691	! qsol(water height in soil) only for bucket continental model
692	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
693	DO j = 1, knon
694	i = ni(j)
695	yqsol(j) = qsol(i)
696	END DO
697	ENDIF
698
699	!****************************************************************************************
700	! 6) Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
701	! atmosphere and coefficients for turbulent diffusion at surface(Cdrag).
702	!
703	!****************************************************************************************
704
705	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
706	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, &
707	ycoefm, ycoefh,ytke)
708
709	!****************************************************************************************
710	!
711	! 8) "La descente" - "The downhill"
712	!
713	! climb_hq_down and climb_wind_down calculate the coefficients
714	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
715	! Only the coefficients at surface for H and Q are returned.
716	!
717	!****************************************************************************************
718
719	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
720	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
721	ydelp, yt, yq, dtime, &
722	petAcoef, peqAcoef, petBcoef, peqBcoef)
723
724	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
725	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv)
726
727
728	!****************************************************************************************
729	! 9) Small calculations
730	!
731	!****************************************************************************************
732
733	ypsref(:) = ypaprs(:,1)
734	epot_air(:) = 0.0
735	epot_air(1:knon) = RCPDyt(1:knon,1)(ypsref(1:knon)/ypplay(1:knon,1))**RKAPPA
736
737	swdown(:) = 0.0
738	IF (nsrf .EQ. is_ter) THEN
739	swdown(1:knon) = ysolsw(1:knon)/(1-yalb(1:knon))
740	ELSE
741	swdown(1:knon) = ysolsw(1:knon)
742	ENDIF
743
744	! constant CO2
745	r_co2_ppm(:) = co2_ppm
746
747	!****************************************************************************************
748	!
749	! 10) Switch selon current surface
750	! It is necessary to start with the continental surfaces because the ocean
751	! needs their run-off.
752	!
753	!****************************************************************************************
754	SELECT CASE(nsrf)
755
756	CASE(is_ter)
757	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
758	rlon, rlat, &
759	debut, lafin, ydelp(:,1), epot_air, r_co2_ppm, ysollwdown, ysolsw, swdown, &
760	yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
761	petAcoef, peqAcoef, petBcoef, peqBcoef, &
762	ypsref, yu1, yv1, yrugoro, pctsrf, &
763	yrads, ysnow, yqsurf, yqsol, yagesno, &
764	ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
765	ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
766
767	CASE(is_lic)
768	CALL surf_landice(itap, dtime, knon, ni, &
769	yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
770	petAcoef, peqAcoef, petBcoef, peqBcoef, &
771	ypsref, yu1, yv1, yrugoro, pctsrf, &
772	yrads, ysnow, yqsurf, yqsol, yagesno, &
773	ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
774	ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
775
776	CASE(is_oce)
777	CALL surf_ocean(rlon, rlat, ysollw, yalb, &
778	yrugos, ywindsp, rmu0, &
779	yfder, &
780	itap, dtime, jour, knon, ni, &
781	debut, swdown, &
782	ypplay(:,1), ycoefh(:,1), ycoefm(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
783	petAcoef, peqAcoef, petBcoef, peqBcoef, &
784	ypsref, yu1, yv1, yrugoro, pctsrf, &
785	yrads, ysnow, yqsurf, yagesno, &
786	yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
787	ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
788
789	CASE(is_sic)
790	CALL surf_seaice( &
791	rlon, rlat, ysollw, yalb, &
792	yfder, &
793	itap, dtime, jour, knon, ni, &
794	debut, lafin, swdown, &
795	yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
796	petAcoef, peqAcoef, petBcoef, peqBcoef, &
797	ypsref, yu1, yv1, yrugoro, pctsrf, &
798	yrads, ysnow, yqsurf, yqsol, yagesno, &
799	ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
800	ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
801
802
803	CASE DEFAULT
804	WRITE(lunout,*) 'Surface index = ', nsrf
805	abort_message = 'Surface index not valid'
806	CALL abort_gcm(modname,abort_message,1)
807	END SELECT
808
809	!****************************************************************************************
810	! Save the fraction of this sub-surface
811	!
812	!****************************************************************************************
813	pctsrf_new(:,nsrf) = pctsrf_nsrf(:)
814
815	!****************************************************************************************
816	! 11) - Calcul the increment of surface temperature
817	! - Update albedo
818	!
819	!****************************************************************************************
820	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
821
822	yalb(1:knon) = yalb_new(1:knon)
823
824	!****************************************************************************************
825	!
826	! 12) "La remontee" - "The uphill"
827	!
828	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
829	! for X=H, Q, U and V, for all vertical levels.
830	!
831	!****************************************************************************************
832	! H and Q
833	print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
834	print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
835	if (ok_flux_surf) then
836	y_flux_t1(:) = fsens
837	y_flux_q1(:) = flat/RLVTT
838	yfluxlat(:) = flat
839	else
840	y_flux_t1(:) = yfluxsens(:)
841	y_flux_q1(:) = -yevap(:)
842	endif
843
844	CALL climb_hq_up(knon, dtime, yt, yq, &
845	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
846	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
847
848	! U and V
849	CALL climb_wind_up(knon, dtime, yu, yv, &
850	y_flux_u, y_flux_v, y_d_u, y_d_v)
851
852	DO j = 1, knon
853	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
854	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
855	yu1(j) = yu1(j) * ypct(j)
856	yv1(j) = yv1(j) * ypct(j)
857	ENDDO
858
859	!****************************************************************************************
860	! 13) Transform variables for output format :
861	! - Decompress
862	! - Multiply with pourcentage of current surface
863	! - Cumulate in global variable
864	!
865	!****************************************************************************************
866
867	tke(:,:,nsrf)=0.
868	DO k = 1, klev
869	DO j = 1, knon
870	i = ni(j)
871	ycoefh(j,k) = ycoefh(j,k) * ypct(j)
872	ycoefm(j,k) = ycoefm(j,k) * ypct(j)
873	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
874	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
875	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
876	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
877
878	flux_t(i,k,nsrf) = y_flux_t(j,k)
879	flux_q(i,k,nsrf) = y_flux_q(j,k)
880	flux_u(i,k,nsrf) = y_flux_u(j,k)
881	flux_v(i,k,nsrf) = y_flux_v(j,k)
882
883	tke(i,k,nsrf)=ytke(j,k)
884
885	ENDDO
886	ENDDO
887
888	evap(:,nsrf) = - flux_q(:,1,nsrf)
889
890	albe(:, nsrf) = 0.
891	alblw(:, nsrf) = 0.
892	snow(:, nsrf) = 0.
893	qsurf(:, nsrf) = 0.
894	rugos(:, nsrf) = 0.
895	fluxlat(:,nsrf) = 0.
896	DO j = 1, knon
897	i = ni(j)
898	d_ts(i,nsrf) = y_d_ts(j)
899	albe(i,nsrf) = yalb(j)
900	alblw(i,nsrf) = yalblw(j)
901	snow(i,nsrf) = ysnow(j)
902	qsurf(i,nsrf) = yqsurf(j)
903	rugos(i,nsrf) = yz0_new(j)
904	fluxlat(i,nsrf) = yfluxlat(j)
905	agesno(i,nsrf) = yagesno(j)
906	cdragh(i) = cdragh(i) + ycoefh(j,1)
907	cdragm(i) = cdragm(i) + ycoefm(j,1)
908	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
909	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
910	zu1(i) = zu1(i) + yu1(j)
911	zv1(i) = zv1(i) + yv1(j)
912	END DO
913
914	IF ( nsrf .EQ. is_ter ) THEN
915	DO j = 1, knon
916	i = ni(j)
917	qsol(i) = yqsol(j)
918	END DO
919	END IF
920
921	ftsoil(:,:,nsrf) = 0.
922	DO k = 1, nsoilmx
923	DO j = 1, knon
924	i = ni(j)
925	ftsoil(i, k, nsrf) = ytsoil(j,k)
926	END DO
927	END DO
928
929
930	#ifdef CRAY
931	DO k = 1, klev
932	DO j = 1, knon
933	i = ni(j)
934	#else
935	DO j = 1, knon
936	i = ni(j)
937	DO k = 1, klev
938	#endif
939	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
940	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
941	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
942	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
943	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
944	#ifdef CRAY
945	END DO
946	END DO
947	#else
948	END DO
949	END DO
950	#endif
951
952	!****************************************************************************************
953	! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
954	! Call HBTM
955	!
956	!****************************************************************************************
957	t2m(:,nsrf) = 0.
958	q2m(:,nsrf) = 0.
959	u10m(:,nsrf) = 0.
960	v10m(:,nsrf) = 0.
961
962	pblh(:,nsrf) = 0. ! Hauteur de couche limite
963	plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA
964	capCL(:,nsrf) = 0. ! CAPE de couche limite
965	oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite
966	cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite
967	pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite
968	therm(:,nsrf) = 0.
969	trmb1(:,nsrf) = 0. ! deep_cape
970	trmb2(:,nsrf) = 0. ! inhibition
971	trmb3(:,nsrf) = 0. ! Point Omega
972
973	#undef T2m
974	#define T2m
975	#ifdef T2m
976	! diagnostic t,q a 2m et u, v a 10m
977
978	DO j=1, knon
979	i = ni(j)
980	uzon(j) = yu(j,1) + y_d_u(j,1)
981	vmer(j) = yv(j,1) + y_d_v(j,1)
982	tair1(j) = yt(j,1) + y_d_t(j,1)
983	qair1(j) = yq(j,1) + y_d_q(j,1)
984	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
985	* (ypaprs(j,1)-ypplay(j,1))
986	tairsol(j) = yts(j) + y_d_ts(j)
987	rugo1(j) = yrugos(j)
988	IF(nsrf.EQ.is_oce) THEN
989	rugo1(j) = rugos(i,nsrf)
990	ENDIF
991	psfce(j)=ypaprs(j,1)
992	patm(j)=ypplay(j,1)
993	qairsol(j) = yqsurf(j)
994	END DO
995
996
997	! Calculate the temperature et relative humidity at 2m and the wind at 10m
998	CALL stdlevvar(klon, knon, nsrf, zxli, &
999	uzon, vmer, tair1, qair1, zgeo1, &
1000	tairsol, qairsol, rugo1, psfce, patm, &
1001	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1002
1003	DO j=1, knon
1004	i = ni(j)
1005	t2m(i,nsrf)=yt2m(j)
1006
1007	q2m(i,nsrf)=yq2m(j)
1008
1009	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
1010	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
1011	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
1012
1013	END DO
1014
1015
1016	CALL HBTM(knon, ypaprs, ypplay, &
1017	yt2m,yt10m,yq2m,yq10m,yustar, &
1018	y_flux_t,y_flux_q,yu,yv,yt,yq, &
1019	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
1020	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1021
1022	DO j=1, knon
1023	i = ni(j)
1024	pblh(i,nsrf) = ypblh(j)
1025	plcl(i,nsrf) = ylcl(j)
1026	capCL(i,nsrf) = ycapCL(j)
1027	oliqCL(i,nsrf) = yoliqCL(j)
1028	cteiCL(i,nsrf) = ycteiCL(j)
1029	pblT(i,nsrf) = ypblT(j)
1030	therm(i,nsrf) = ytherm(j)
1031	trmb1(i,nsrf) = ytrmb1(j)
1032	trmb2(i,nsrf) = ytrmb2(j)
1033	trmb3(i,nsrf) = ytrmb3(j)
1034	END DO
1035
1036	#else
1037	! not defined T2m
1038	! No calculation
1039	! Set output variables to zero
1040	DO j = 1, knon
1041	i = ni(j)
1042	pblh(i,nsrf) = 0.
1043	plcl(i,nsrf) = 0.
1044	capCL(i,nsrf) = 0.
1045	oliqCL(i,nsrf) = 0.
1046	cteiCL(i,nsrf) = 0.
1047	pblT(i,nsrf) = 0.
1048	therm(i,nsrf) = 0.
1049	trmb1(i,nsrf) = 0.
1050	trmb2(i,nsrf) = 0.
1051	trmb3(i,nsrf) = 0.
1052	END DO
1053	DO j = 1, knon
1054	i = ni(j)
1055	t2m(i,nsrf)=0.
1056	q2m(i,nsrf)=0.
1057	u10m(i,nsrf)=0.
1058	v10m(i,nsrf)=0.
1059	END DO
1060	#endif
1061
1062	!****************************************************************************************
1063	! 15) End of loop over different surfaces
1064	!
1065	!****************************************************************************************
1066	END DO loop_nbsrf
1067
1068	!****************************************************************************************
1069	! 16) Calculate the mean value over all sub-surfaces for som variables
1070	!
1071	! NB!!! jg : Pour garder la convergence numerique j'utilise pctsrf_new comme c'etait
1072	! fait dans physiq.F mais ca devrait plutot etre pctsrf???!!!!! A verifier!
1073	!****************************************************************************************
1074
1075	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1076	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1077	DO nsrf = 1, nbsrf
1078	DO k = 1, klev
1079	DO i = 1, klon
1080	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf_new(i,nsrf)
1081	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf_new(i,nsrf)
1082	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf_new(i,nsrf)
1083	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf_new(i,nsrf)
1084	END DO
1085	END DO
1086	END DO
1087
1088	DO i = 1, klon
1089	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1090	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1091	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1092	ENDDO
1093
1094
1095	DO i = 1, klon
1096	IF ( ABS( pctsrf_new(i, is_ter) + pctsrf_new(i, is_lic) + &
1097	pctsrf_new(i, is_oce) + pctsrf_new(i, is_sic) - 1.) .GT. EPSFRA) &
1098	THEN
1099	WRITE(,) 'physiq : pb sous surface au point ', i, &
1100	pctsrf_new(i, 1 : nbsrf)
1101	ENDIF
1102	ENDDO
1103
1104	!
1105	! Incrementer la temperature du sol
1106	!
1107	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
1108	zt2m(:) = 0.0 ; zq2m(:) = 0.0
1109	zu10m(:) = 0.0 ; zv10m(:) = 0.0
1110	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
1111	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1112	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1113	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1114	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1115
1116
1117	DO nsrf = 1, nbsrf
1118	DO i = 1, klon
1119	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1120
1121	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1122	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf_new(i,nsrf)
1123	wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1124	pctsrf_new(i,nsrf)
1125
1126	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf_new(i,nsrf)
1127	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf_new(i,nsrf)
1128
1129	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf_new(i,nsrf)
1130	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf_new(i,nsrf)
1131	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf_new(i,nsrf)
1132	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf_new(i,nsrf)
1133
1134	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf_new(i,nsrf)
1135	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf_new(i,nsrf)
1136	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf_new(i,nsrf)
1137	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf_new(i,nsrf)
1138	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf_new(i,nsrf)
1139	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf_new(i,nsrf)
1140	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf_new(i,nsrf)
1141	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf_new(i,nsrf)
1142	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf_new(i,nsrf)
1143	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf_new(i,nsrf)
1144	END DO
1145	END DO
1146
1147	IF (check) THEN
1148	amn=MIN(ts(1,is_ter),1000.)
1149	amx=MAX(ts(1,is_ter),-1000.)
1150	DO i=2, klon
1151	amn=MIN(ts(i,is_ter),amn)
1152	amx=MAX(ts(i,is_ter),amx)
1153	ENDDO
1154	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1155	ENDIF
1156	!
1157	! If a sub-surface does not exsist for a grid point, the mean value for all
1158	! sub-surfaces is distributed.
1159	!
1160	DO nsrf = 1, nbsrf
1161	DO i = 1, klon
1162	IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1163	ts(i,nsrf) = zxtsol(i)
1164	t2m(i,nsrf) = zt2m(i)
1165	q2m(i,nsrf) = zq2m(i)
1166	u10m(i,nsrf) = zu10m(i)
1167	v10m(i,nsrf) = zv10m(i)
1168
1169	! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1170	pblh(i,nsrf) = s_pblh(i)
1171	plcl(i,nsrf) = s_plcl(i)
1172	capCL(i,nsrf) = s_capCL(i)
1173	oliqCL(i,nsrf) = s_oliqCL(i)
1174	cteiCL(i,nsrf) = s_cteiCL(i)
1175	pblT(i,nsrf) = s_pblT(i)
1176	therm(i,nsrf) = s_therm(i)
1177	trmb1(i,nsrf) = s_trmb1(i)
1178	trmb2(i,nsrf) = s_trmb2(i)
1179	trmb3(i,nsrf) = s_trmb3(i)
1180	ENDIF
1181	ENDDO
1182	ENDDO
1183
1184
1185	DO i = 1, klon
1186	fder(i) = - 4.0RSIGMAzxtsol(i)**3
1187	ENDDO
1188
1189	zxqsurf(:) = 0.0
1190	zxsnow(:) = 0.0
1191	DO nsrf = 1, nbsrf
1192	DO i = 1, klon
1193	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf_new(i,nsrf)
1194	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf_new(i,nsrf)
1195	END DO
1196	END DO
1197
1198	!
1199	!IM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
1200	!IM ajout dependance type surface
1201	rh2m(:) = 0.0
1202	qsat2m(:) = 0.0
1203
1204	DO i = 1, klon
1205	DO nsrf=1, nbsrf
1206	zx_t1(i,nsrf) = t2m(i,nsrf)
1207	IF (thermcep) THEN
1208	zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
1209	zx_qs1(i,nsrf) = r2es * &
1210	FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
1211	zx_qs1(i,nsrf) = MIN(0.5,zx_qs1(i,nsrf))
1212	zcor1(i,nsrf) = 1./(1.-retv*zx_qs1(i,nsrf))
1213	zx_qs1(i,nsrf) = zx_qs1(i,nsrf)*zcor1(i,nsrf)
1214	END IF
1215	zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
1216	zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
1217	rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf_new(i,nsrf)
1218	qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf_new(i,nsrf)
1219	END DO
1220	END DO
1221
1222	! Some of the module declared variables are returned for printing in physiq.F
1223	qsol_d(:) = qsol(:)
1224	evap_d(:,:) = evap(:,:)
1225	rugos_d(:,:) = rugos(:,:)
1226	agesno_d(:,:) = agesno(:,:)
1227
1228
1229	END SUBROUTINE pbl_surface
1230	!
1231	!****************************************************************************************
1232	!
1233	SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1234	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1235
1236	INCLUDE "indicesol.h"
1237	INCLUDE "dimsoil.h"
1238
1239	! Ouput variables
1240	!****************************************************************************************
1241	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
1242	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
1243	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
1244	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
1245	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
1246	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
1247	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
1248	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1249
1250
1251	!****************************************************************************************
1252	! Return module variables for writing to restart file
1253	!
1254	!****************************************************************************************
1255	qsol_rst(:) = qsol(:)
1256	fder_rst(:) = fder(:)
1257	snow_rst(:,:) = snow(:,:)
1258	qsurf_rst(:,:) = qsurf(:,:)
1259	evap_rst(:,:) = evap(:,:)
1260	rugos_rst(:,:) = rugos(:,:)
1261	agesno_rst(:,:) = agesno(:,:)
1262	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1263
1264	!****************************************************************************************
1265	! Deallocate module variables
1266	!
1267	!****************************************************************************************
1268	DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1269
1270	END SUBROUTINE pbl_surface_final
1271	!
1272	!****************************************************************************************
1273	!
1274
1275	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format