Context Navigation

pbl_surface_mod.F90 @ 3962

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

correction du bug ISCCP (n'ecrire ptop que quand ISCCP est appelé)
petite inversion de boucle dans isccp_cloud_types.F pour aller + vite
"CFisation" d'un certain nombre d'unités pour les hist*
les suggestions de JL pour rugoro

SD
LF

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 48.0 KB

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

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

Download in other formats:

Original Format