Context Navigation

pbl_surface_mod.F90 @ 785

Last change on this file since 785 was 781, checked in by Laurent Fairhead, 18 years ago
Routines pour la nouvelle interface avec les surfaces de Josefine JG LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 49.2 KB

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

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

Download in other formats:

Original Format