Context Navigation

pbl_surface_mod.F90 @ 1080

Last change on this file since 1080 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

Rev	Line
[814]	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
[815]	15	USE surface_data, ONLY : ocean, ok_veget
[814]	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, &
[852]	199	zxfluxt, zxfluxq, q2m, flux_q, tke )
[814]	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)
[852]	242	! tke---input/output-R- tke (kg/m**2/s)
[814]	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
[852]	352	! Input/output
	353	REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke
[814]	354
[852]	355
[814]	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
[852]	427	REAL, DIMENSION(klon,klev+1) :: ytke
[814]	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
[852]	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
[814]	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
[852]	526	ytherm = 0.0 ; ytke=0.
[814]	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)
[852]	669	ytke(j,k)=tke(i,k,nsrf)
[814]	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, &
[852]	700	ycoefm, ycoefh,ytke)
[814]	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]	852	tke(:,:,nsrf)=0.
[814]	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)
[852]	867
	868	tke(i,k,nsrf)=ytke(j,k)
	869
[814]	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