Context Navigation

module_lmdz_phys.F @ 13

Last change on this file since 13 was 6, checked in by lfita, 11 years ago

Fixing problem with the solar radiation:

1.- Correct values for the julian and the hour (0,1) to the LMDZ physiq

Adding NaN checking for 'PSFC' and 'T_2'

File size: 88.8 KB

Rev	Line
[1]	1	MODULE module_lmdz_phys
	2	!
	3	! Module to include LMDZ (http://lmdz.lmd.jussieu.fr/) physical packages in WRF
	4	!
	5	! L. Fita, Laboratoire Meterologie Dynamique, IPSL, UPMC, CNRS, Jussieu, Paris, France
	6	! July 2013
	7	!
	8	CONTAINS
	9
	10	SUBROUTINE call_lmdz_phys(wrf_grid, wrf_xtime, wrf_restart_alarm, &
	11	& wrf_Lon, wrf_Lat, &
	12	& wrf_T, wrf_U, wrf_V, &
	13	& wrf_perP, wrf_baseP, &
	14	& wbdyw, wrf_isrestart, &
	15	& wrf_islowbdyin &
	16	! dimension arguments
	17	& ,wids,wide, wjds,wjde, wkds,wkde &
	18	& ,wims,wime, wjms,wjme, wkms,wkme &
	19	& ,wips,wipe, wjps,wjpe, wkps,wkpe &
	20	& ,wi_start,wi_end &
	21	& ,wj_start,wj_end &
	22	& ,wkts, wkte &
	23	& ,wnum_tiles &
	24	& ,wnum3dm, wparfirstscal, &
	25	& wnx, wny, wnz, &
	26	& wjulday, wgmt, wtime_step, &
	27	& wrf_fulleta, wrf_halfeta, wrf_dfulleta, &
	28	& wrf_fullpres, wrf_pergeopot, wrf_basegeopot, &
	29	& wrf_MOIST, wrf_w, wrf_ptop, &
	30	& wrf_permass, wrf_basemass, &
	31	& wrf_mut, wrf_muu, wrf_muv, &
	32	& wrf_Utend, wrf_Vtend, wrf_Ttend, &
	33	& wrf_MOISTtend, wrf_psfctend, &
	34	& wrf_qvid, wrf_qcid, wrf_qrid, &
	35	& wrf_qsid, wrf_qiid, wrf_qhid, wrf_qgid, &
	36	& wrf_mapft, wrf_mapfu, wrf_mapfv, wrf_dx, wrf_dy, &
	37	& wrf_dbg, landcat, soilcat, &
	38	& wrf_l_pbl, wrf_l_con, wrf_l_thermals, wrf_l_wake, &
	39	& wrf_nsoillayers, &
	40	& icheck_p, jcheck_p, kcheck_p &
	41	& )
	42
	43	! WRF modules
	44	USE module_model_constants
	45	USE module_domain_type
	46
	47	! LMDZ modules
	48	USE infotrac, ONLY: infotrac_init, type_trac
	49	USE comgeomphy
	50	USE control_mod
	51	USE indice_sol_mod
	52	USE phys_state_var_mod
	53	! USE fonte_neige_mod, ONLY: fonte_neige_init, fonte_neige_final
	54	USE fonte_neige_mod
	55	USE iostart, ONLY: get_var
	56	USE surface_data
	57	USE pbl_surface_mod
	58	USE mod_grid_phy_lmdz, only : klon_glo
	59
	60	! WRF+LMDZ
	61
	62	! Used to initialize variables in an old fashion way... Now using 'lmdz_wrf_variables_mod'
	63	USE wrf_lmdz_mod
	64	USE lmdz_wrf_variables_mod
	65	USE NOread_limit_sub_variables
	66	USE output_lmdz_NOmodule
	67
	68	IMPLICIT NONE
	69
	70	INCLUDE "../lmdz/dimensions.h"
	71	INCLUDE "../lmdz/comvert.h"
	72	INCLUDE "../lmdz/temps.h"
	73	INCLUDE "../lmdz/dimsoil.h"
	74
	75	TYPE(domain) , TARGET :: wrf_grid
	76
	77	REAL, INTENT(IN) :: wrf_xtime, wrf_Ptop, &
	78	wgmt, wrf_dx, wrf_dy
	79	LOGICAL, INTENT(IN) :: wrf_restart_alarm, &
	80	wrf_isrestart, wrf_islowbdyin
	81
	82	INTEGER, INTENT(IN) :: wnx, wny, wnz, wbdyw
	83
	84	INTEGER, INTENT(IN) :: &
	85	& wids,wide, wjds,wjde, wkds,wkde, &
	86	wims,wime, wjms,wjme, wkms,wkme, &
	87	wips,wipe, wjps,wjpe, wkps,wkpe, &
	88	wkts,wkte, wnum_tiles
	89
	90	INTEGER, DIMENSION(wnum_tiles), INTENT(IN) :: &
	91	& wi_start,wi_end,wj_start,wj_end
	92	INTEGER, INTENT(IN) :: wparfirstscal, wnum3dm
[6]	93	REAL, INTENT(IN) :: wjulday
[1]	94	REAL, INTENT(IN) :: wtime_step
	95	INTEGER, INTENT(IN) :: wrf_qvid, wrf_qcid, &
	96	& wrf_qrid, wrf_qsid, wrf_qiid, wrf_qhid, wrf_qgid, wrf_dbg
	97
	98	REAL, INTENT(IN), DIMENSION(wkms:wkme) :: wrf_fulleta
	99	REAL, INTENT(IN), DIMENSION(wkms:wkme) :: wrf_halfeta, wrf_dfulleta
	100
	101	REAL, INTENT(IN), DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: &
	102	& wrf_T, wrf_perP, wrf_baseP, wrf_U, wrf_V
	103
	104	REAL, INTENT(OUT), DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: &
	105	& wrf_Ttend, wrf_Utend, wrf_Vtend
	106
	107	REAL, INTENT(IN), DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: &
	108	& wrf_fullpres, wrf_pergeopot, wrf_basegeopot, wrf_w
	109
	110	REAL, INTENT(IN), DIMENSION(wims:wime,wjms:wjme) :: wrf_Lon, wrf_Lat, &
	111	& wrf_basemass, wrf_permass, wrf_mapft, wrf_mapfu, wrf_mapfv, wrf_mut, &
	112	& wrf_muu, wrf_muv
	113
	114	REAL, INTENT(OUT), DIMENSION(wims:wime,wjms:wjme) :: wrf_psfctend
	115
	116	REAL, INTENT(IN), DIMENSION(wims:wime,wkms:wkme,wjms:wjme,1:wnum3dm) :: &
	117	& wrf_MOIST
	118
	119	REAL, INTENT(OUT), DIMENSION(wims:wime,wkms:wkme,wjms:wjme,1:wnum3dm) :: &
	120	& wrf_MOISTtend
	121
	122	CHARACTER(LEN=256), INTENT(IN) :: landcat, soilcat
	123	INTEGER, INTENT(IN) :: wrf_l_pbl, wrf_l_con, &
	124	wrf_l_thermals, wrf_l_wake
	125	INTEGER, INTENT(IN) :: wrf_nsoillayers
	126	INTEGER, INTENT(IN) :: icheck_p, jcheck_p, &
	127	kcheck_p
	128	! Local
	129	INTEGER :: ip, jp, kp, lp
	130
	131	REAL, DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: wrf_Temp, &
	132	& wrf_Temptend, wrf_P, wrfMASSValues
	133	REAL, DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: wrfUdestagg, wrfVdestagg,&
	134	& wrfUtenddestagg, wrfVtenddestagg
	135	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw), wkms:wkme-1) :: &
	136	& TValues, TtendValues, HalfPressValues, WFluxMassValues
	137	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw), wkms:wkme-1) :: & !! DE-STAGGED !!
	138	& UValues, UtendValues
	139	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw), wkms:wkme-1) :: & !! DE-STAGGED !!
	140	& VValues, VtendValues
	141	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw), wkms:wkme) :: & !! Z-STAGG !!
	142	& FullPressValues, GeopotValues
	143	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw)) :: &
	144	& sfcGeopotValues
	145	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: MixingRatioValues, &
	146	& MixingRatiotendValues
	147
	148	REAL, DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: wrfWValues, &
	149	& wrfPERMASSValues, wrfZFluxMassValues
	150	REAL, DIMENSION(wims:wime,wkms:wkme,wjms:wjme) :: wrf_geopot, w3Dmat
	151	REAL, DIMENSION(wims:wime,wjms:wjme) :: wrf_sfcgeopot
	152
	153	LOGICAL :: lmdz_debut, lmdz_lafin
	154
	155	INTEGER :: k, ix, iy, iz, iq, ixy, &
	156	& lmdzp, ilmdz, jlmdz, ir, ixx, iyy
	157	INTEGER :: dimx, dimy, dimz, dimxy
	158	INTEGER :: ddimx, ddimy, ddimxy, ddimxy2
	159	INTEGER :: dim_klon_glo, dim_ngridmx
	160	INTEGER :: im2, jm2, km2
	161	INTEGER :: lmdzmixingratios
	162	CHARACTER (LEN=1024) :: message, wrf_err_message
	163
	164	INTEGER, PARAMETER :: Nphykeys=20
	165	REAL, DIMENSION(Nphykeys) :: phykeys
	166
	167	REAL, DIMENSION(2,(wime-wims-2wbdyw)(wjme-wjms-2*wbdyw),wkms:wkme-1) :: &
	168	& lmdz_dudyn, l3Dmat
	169	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw),3) :: lmdz_PVtheta
	170	REAL, DIMENSION(wkms:wkme-1) :: lmdzoutP
	171	REAL, DIMENSION((wime-wims-2wbdyw),(wjme-wjms-2wbdyw)) :: wrf_area
	172	REAL, DIMENSION((wime-wims-2wbdyw),(wjme-wjms-2wbdyw)) :: wrf_cu
	173	REAL, DIMENSION((wime-wims-2wbdyw),(wjme-wjms-2wbdyw)) :: wrf_cv
	174	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw)) :: AreaValues, &
	175	LonValues, LatValues, PsfctendValues
	176	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw)) :: CUValues, CVValues
	177	REAL, DIMENSION((wime-wims-2wbdyw)(wjme-wjms-2*wbdyw)) :: runoff_lic
	178	LOGICAL :: existsfile
	179	INTEGER :: ierr
	180	CHARACTER(LEN=50) :: errmsg, wrnmsg, subname
	181	CHARACTER(LEN=100) :: defilen, defvarn, defstr
	182	INTEGER :: funit, defint
	183	LOGICAL :: is_used, defbool
	184	REAL, DIMENSION(20) :: lmdz_clesphy0
	185	INTEGER :: iflag_phys
	186	REAL :: defreal,ksoilmax
	187	REAL, DIMENSION(4) :: ksoils
	188	REAL :: Pi
	189	CHARACTER(LEN=50) :: LMDZvarmethod, Spt, Spt3D
	190	CHARACTER(LEN=4) :: ipS, jpS, kpS
	191	CHARACTER(LEN=8) :: lpS
	192	INTEGER :: nlanduse
	193	REAL, DIMENSION(11) :: ORCHIDEElevels
	194	INTEGER :: Nsoillayers,begdataset, &
	195	Nvars,Nrs,Ncs,Nwrfsoiltypes
	196	REAL :: soilTmass
	197	CHARACTER(LEN=50) :: datasetname
	198	CHARACTER(LEN=100) :: tablename
	199	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: soilvals
	200	CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:) :: soilvalns
	201	INTEGER :: lpbl, lcon, lth, lwake
	202	TYPE(WRFU_Alarm), TARGET, DIMENSION(MAX_WRF_ALARMS) :: wrf_alarms
	203
	204	!!!!!!! Variables
	205	!
	206	! wrf_xtime: WRF actual timestep (in minutes since simulation started)
	207	! wrf_restart_alarm: WRF value indicating whether this is a restart time-step
	208	! wrf_T: WRF 3D temperature
	209	! wrf_perP: WRF 3D perturbation of the pressure
	210	! wrf_baseP: WRF 3D base pressure
	211	! wn[x/y/z]: WRF dimensions at each axis
	212	! wbdyw: number of grid points in the sponge zone
	213	! wids: start for i in domain
	214	! wide: end for i in domain
	215	! wjds: start for j in domain
	216	! wjde: end for j in domain
	217	! wkds: start for k in domain
	218	! wkde: end for k in domain
	219	! wims: start for i in memory
	220	! wime: end for i in memory
	221	! wjms: start for j in memory
	222	! wjme: end for j in memory
	223	! wips: start for i in patch
	224	! wipe: end for i in patch
	225	! wjps: start for j in patch
	226	! wjpe: end for j in patch
	227	! wkms: start for k in memory
	228	! wkme: end for k in memory
	229	! wi_start: start indices for i in tile
	230	! wi_end: end indices for i in tile
	231	! wj_start: start indices for j in tile
	232	! wj_end: end indices for j in tile
	233	! wkts: start for k in tile
	234	! wkte: end for k in tile
	235	! wnum_tiles: number of tiles
	236	! wjulday: WRF julian day
[6]	237	! wgmt: WRF gmt hour (0, 1)
[1]	238	! wtime_step: WRF time-step (in seconds)
	239	! wrf_dbg: level of debug from WRF
	240	! wrf_isrestart: flag indicating whether simulation is beggining from a restart or
	241	! not (zero time since simulation started )
	242	!
	243	! lmdz_debut: flag indicating if this is the first step for LMDZ
	244	! lmdz_lafin: flag indicating if this is the last step for LMDZ
	245	! wrf_fulleta: WRF full mass eta levels (w)
	246	! wrf_dfulleta: WRF full mass high between eta levels (w)
	247	! wrf_halfeta: WRF half mass eta levels (t, u, v, h)
	248	! wrf_Ptop: WRF top pressure
	249	! wrf_Lon: WRF longitudes
	250	! wrf_Lat: WRF latitudes
	251	! wrf_T: WRF potential temperature
	252	! wrf_Ttend: WRF potential temperature tendency
	253	! wrf_Temp: WRF temperature
	254	! wrf_Temptend: WRF temperature tendency
	255	! wrf_perP: WRF perturbation of pressure
	256	! wrf_baseP: WRF base of pressure
	257	! wrf_P: WRF pressure
	258	! wrf_U: WRF x-wind speed
	259	! wrf_Utend: WRF x-wind speed tendency
	260	! wrf_V: WRF y-wind speed
	261	! wrf_Vtend: WRF y-wind speed tendency
	262	! wrf_fullpres: WRF pressure in full eta levels
	263	! wrf_pergeopot: WRF perturbation of geopotential
	264	! wrf_basegeopot: WRF base of geopotential
	265	! wrf_geopot: WRF geopotential
	266	! wrf_psfctend: WRF surface pressure tendency
	267	! wrf_sfcgeopot: WRF surface geopotential
	268	! wrf_MOIST: WRF MOIST ARRAY (water spieces and tracers)
	269	! wrf_MOISTtend: WRF MOIST tendency ARRAY (water spieces and tracers)
	270	! wrf_qvid: WRF water vapor index of wrf_MOIST
	271	! wrf_qcid: WRF cloud liquid water index of wrf_MOIST
	272	! wrf_qrid: WRF rain liquid water index of wrf_MOIST
	273	! wnum3dm: WRF number of water species
	274	! wrf_w: WRF vertical velocity
	275	! wrf_permass: WRF perturbation dry air mass in column
	276	! wrf_basemass: WRF base dry air mass in column
	277	! wrf_mut: WRF total dry air mass column on mass points
	278	! wrf_muu: WRF total dry air mass column on u-wind grid points
	279	! wrf_muv: WRF total dry air mass column on v-wind gid points
	280	! wrfWValues: staggered WRF vertical velocity (in m s-1)
	281	! wrfPERMASSValues: de-staggered WRF perturbation of dry mass column (in Pa)
	282	! wrfMASSValues: de-staggered WRF perturbation of dry mass column (in kg)
	283	! wrfZFluxMassValues: WRF vertical flux mass values (kg s-1)
	284	! wrf_mapf[t/u/v]: WRF map factor for Temperature, u-wind and v-wind grid points
	285	! wrf_d[x/y]: length of the grid spacing in the x and y direction (in m)
	286	! wrf_area: WRF area of the grid cells
	287	! wrf_c[u/v]: WRF covariant values for wind components
	288	! landcat: land category used in the simulation
	289	! nlanduse: number of soil types in the WRF landuse variable
	290	!
	291	! LonValues: 2D-LMDZ array of WRF longitudes (in rad)
	292	! LatValues: 2D-LMDZ array of WRF latitudes (in rad)
	293	! TValues: 2D-LMDZ array of WRF temperatures (in K)
	294	! TtendValues: 2D-LMDZ tendency array of WRF temperatures (in K)
	295	! HalfPressValues: 2D-LMDZ array of WRF pressure at half-mass levels (in Pa)
	296	! UValues: 2D-LMDZ array of WRF x-wind speed (in m s-1)
	297	! UtendValues: 2D-LMDZ array of WRF tendency x-wind speed (in m s-2)
	298	! VValues: 2D-LMDZ array of WRF y-wind speed (in m s-1)
	299	! VtendValues: 2D-LMDZ array of WRF tendency y-wind speed (in m s-2)
	300	! PsfctendValues: 2D-LMDZ array of surface pressure tendency y-wind speed (in Pa s-1)
	301	! FullPressValues: 2D-LMDZ array of WRF pressure at full-mass levels (in Pa)
	302	! GeopotValues: 2D-LMDZ array of WRF geopotential (in m2 s-2)
	303	! sfcGeopotValues: LMDZ array of WRF sfc geopotential (in m2 s-2)
	304	! lmdzmixingratios: number of 'mixing ratios' to pass to the LMDZ physics
	305	! MixingRatioValues: 2D-LMDZ array of Mixing Ratios (in kg kg-1)
	306	! MixingRatiotendValues: 2D-LMDZ array of tendency Mixing Ratios (in kg kg-1 s-1)
	307	! WFluxMassValues: 2Dvertical mass flux values (in kg m s-1)
	308	! phykeys: LMDZ obsolete control I/O vector
	309	! lmdz_dudyn: LMDZ obsolete vector
	310	! lmdz_PVtheta: LMDZ obsolete PV values at given theta levels
	311	! lmdzoutP: LMDZ output pressure values (using same eta levels, but imposing Psfc=1000 hPa)
	312	! preff: refernce pressure for lmdzouP
	313	! AreaValues: LMDZ area (m2)
	314	! CUArea: LMDZ cu coeff. (u_covariant = cu * u)
	315	! CVArea: LMDZ cv coeff. (v_covariant = cv * v)
	316	!
	317	!!!!!!! Subroutines/Functions
	318	! WRFmoist_LMDZmoist: Subroutine to transform from WRF moisture array to the LMDZ moisture array
	319	! eta_to_pressure: Subroutine to transform from eta levels to pressure levels
	320	! vect_mat: Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
	321	! mat_vect: Subroutine to transform a 3D matrix to a LMDZ physic 1D vector
	322	! temp_to_theta: Subroutine to transform from temperatures to WRF potential temperatures
	323	! theta_to_temp: Subroutine to transform from WRF potential temperatures to temperature
	324	! physiq: Subroutine with the LMDZ physics
	325
	326	! CAPITAL variables stand for the name of the variables inside the subroutines/functions
	327
	328	!!!!!!! !!!!!! !!!!! !!!! !!! !! !
	329	WRITE(message, *)' starting: module_lmdz_phys.F'
	330	CALL wrf_debug(75, message)
	331
	332	PRINT *,' --------------------------------------------------------------------------------------------------'
	333
	334	! Basic definitions
	335	!!
	336	errmsg='ERROR -- error -- ERROR -- error'
	337	wrnmsg='WARNING -- warning -- WARNING -- warning'
	338	subname='module_lmdz_phys'
	339
	340	Pi = 2.*ASIN(1.)
	341
	342	type_trac = 'lmdz'
	343
	344	dimx=wime-wims
	345	dimy=wjme-wjms
	346	dimz=wkme-wkms
	347	dimxy=dimx*dimy
	348	dim_klon_glo = dimx*dimy
	349	dim_ngridmx = dimx*dimy
	350	ddimx = dimx-2*wbdyw
	351	ddimy = dimy-2*wbdyw
	352	ddimxy=ddimx*ddimy
	353	ddimxy2=ddimx*(ddimy/2-1)+ddimx/2
	354
	355	! Point to retrieve values for debugging
	356	! ip=ddimx/2
	357	! jp=ddimy/2
	358	! kp=dimz/2
	359	! ip=5
	360	! jp=11
	361	ip=icheck_p
	362	jp=jcheck_p
	363	kp=kcheck_p
	364
	365	lp=ddimx*(jp-1)+ip
	366	WRITE(ipS,'(I4)')ip
	367	WRITE(jpS,'(I4)')jp
	368	WRITE(kpS,'(I4)')kp
	369	WRITE(lpS,'(I8)')lp
	370	Spt = TRIM(ipS) // ', '//TRIM(jpS)
	371	Spt3D = TRIM(ipS) // ', '//TRIM(jpS)//', '//TRIM(kpS)
	372	PRINT *,' checking point ip: ',ip,' jp: ',jp,' kp: ',kp,' lp: ',lp
	373
	374	PRINT *,' Lluis in module_lmdz_phy itime_step= ', wrf_grid%itimestep, &
	375	' wrf_xtime: ',wrf_xtime,' wrf_rst: ', wrf_isrestart,' wrf_lwbdy: ', &
	376	wrf_islowbdyin
	377
	378	phykeys=0.
	379	lmdz_dudyn=0.
	380	lmdz_PVtheta=0.
	381
	382	! Checking of first/last time-step for LMDZ
	383	!!
	384	lmdz_debut = .FALSE.
	385	lmdz_lafin = .FALSE.
	386
	387	DO iy=1,ddimy
	388	DO ix=1,ddimx
	389	IF (wrf_Lon(ix,iy) < 0.) THEN
	390	LonValues(ddimx(iy-1)+ix) = (360. + wrf_Lon(ix,iy))Pi/180.
	391	ELSE
	392	LonValues(ddimx(iy-1)+ix) = wrf_Lon(ix,iy)Pi/180.
	393	END IF
	394	LatValues(ddimx(iy-1)+ix) = wrf_Lat(ix,iy)Pi/180.
	395	PsfctendValues(ddimx*(iy-1)+ix) = wrf_psfctend(ix,iy)
	396	END DO
	397	END DO
	398
	399	! Reshape can not be used? Different initial/ending indices !
	400	! LonValues = RESHAPE(wrf_Lon, (/ dimx*dimy/), ORDER=(/1/))
	401	! LatValues = RESHAPE(wrf_Lat, (/ dimx*dimy/), ORDER=(/1/))
	402
	403	firststep: IF ( wrf_xtime == 0 .OR. wrf_isrestart) THEN
	404	lmdz_debut = .TRUE.
	405
	406	DO funit=10,100
	407	INQUIRE(UNIT=funit, OPENED=is_used)
	408	IF (.NOT. is_used) EXIT
	409	END DO
	410	lmdz_clesphy0 = 0.
	411
	412	CALL conf_gcm(funit, .TRUE., lmdz_clesphy0)
	413	! L. Fita, wee need to provide the value
	414	annee_ref = anneeref
	415	day_ref = dayref
	416
	417	! Uploading LMDZ general configuration
	418	! used instead call conf_gcm, since variables are not found in memory (or at least with a value)
	419	defilen = 'gcm.def'
	420	defvarn = 'iflag_phys'
	421	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	422	iflag_phys = defint
	423
	424	! Is not working it has an space...
	425	! defvarn = 'type_trac'
	426	! CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	427	! type_trac = TRIM(defstr)
	428
	429	! used instead call conf_gcm, since variables are not found in memory (or at least with a value)
	430	defilen = 'config.def'
	431	defvarn = 'type_ocean'
	432	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	433	type_ocean= TRIM(defstr)
	434
	435	! used instead call conf_gcm, since variables are not found in memory (or at least with a value)
	436	! Looking for consistency for some namelist variables and *.def LMDZ values
	437	PRINT ,' Checking consistency between namelist and LMDZ .def values'
	438	defilen = 'physiq.def'
	439	defvarn = 'iflag_pbl'
	440	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	441	lpbl = defint
	442	defvarn = 'iflag_con'
	443	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	444	lcon = defint
	445	defvarn = 'iflag_thermals'
	446	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	447	lth = defint
	448	defvarn = 'iflag_wake'
	449	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	450	lwake = defint
	451
	452	IF (wrf_l_pbl /= lpbl) THEN
	453	PRINT *,TRIM(errmsg)
	454	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ pbl value: ', &
	455	wrf_l_pbl,' differs from iflag_pbl= ',lpbl ,' in physiq.def !!!!'
	456	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
	457	CALL wrf_error_fatal(TRIM(wrf_err_message))
	458	END IF
	459	IF (wrf_l_con /= lcon) THEN
	460	PRINT *,TRIM(errmsg)
	461	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ convection value: ', &
	462	wrf_l_con,' differs from iflag_con= ',lcon ,' in physiq.def !!!!'
	463	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
	464	CALL wrf_error_fatal(TRIM(wrf_err_message))
	465	END IF
	466	IF (wrf_l_thermals /= lth) THEN
	467	PRINT *,TRIM(errmsg)
	468	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ thermals value: ', &
	469	wrf_l_thermals,' differs from iflag_thermals= ',lth ,' in physiq.def !!!!'
	470	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
	471	CALL wrf_error_fatal(TRIM(wrf_err_message))
	472	END IF
	473	IF (wrf_l_wake /= lwake) THEN
	474	PRINT *,TRIM(errmsg)
	475	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ wakes value: ', &
	476	wrf_l_wake,' differs from iflag_wake= ',lwake ,' in physiq.def !!!!'
	477	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
	478	CALL wrf_error_fatal(TRIM(wrf_err_message))
	479	END IF
	480
	481	! Initialization of the variables
	482	!!
	483
	484	INQUIRE(FILE='traceur.def', EXIST=existsfile)
	485
	486	! L. Fita, LMD. August, this would be much more nice if it is in config_flags structure....
	487	IF (existsfile) THEN
	488	DO funit=10,100
	489	INQUIRE(UNIT=funit, OPENED=is_used)
	490	IF (.NOT. is_used) EXIT
	491	END DO
	492	OPEN(UNIT=funit, FILE='traceur.def', STATUS='OLD')
	493	READ(funit,*)lmdzmixingratios
	494	CLOSE(UNIT=funit)
	495	ELSE
	496	WRITE(message, *)errmsg
	497	CALL wrf_debug(75, message)
	498	WRITE(wrf_err_message,*)' required file traceur.def not found !!'
	499	CALL wrf_error_fatal(TRIM(wrf_err_message))
	500	END IF
	501
	502	CALL infotrac_init
	503
	504	! Initializing LMDZ matrices...
	505
	506	CALL Init_Phys_lmdz(1,ddimxy,dimz,1,(/ddimxy/))
	507	CALL InitComgeomphy()
	508
	509	CALL phys_state_var_init(0)
	510	!! CALL neige_initialize()
	511	!! CALL limit_initialize()
	512
	513	!! CALL pbl_initialize()
	514
	515	PRINT *,'Lluis after initializing all LMDZ!'
	516
	517	rlon=LonValues
	518	rlat=LatValues
	519
	520	xtime0: IF (wrf_xtime == 0.) THEN
	521
	522	nlanduse=UBOUND(wrf_grid%landusef, 2)
	523	CALL WRFlanduse_LMDZKsoil(landcat,wims,wime,wjms,wjme,ddimx,ddimy,nlanduse, &
	524	wrf_grid%landusef,wrf_grid%landmask,wrf_grid%lter,wrf_grid%llic, &
	525	wrf_grid%loce,wrf_grid%lsic)
	526
	527	! LANDMASK in LMDZ as fractions
	528	wrf_grid%lmsk = wrf_grid%lter + wrf_grid%llic
	529	! Pseudo WRF-initialization
	530
	531	! LANDMASK in LMDZ as 1/0
	532	! wrf_grid%lmsk=0.
	533	! DO ix=1,ddimx
	534	! DO iy=1,ddimy
	535	! IF (wrf_grid%lter(ix,iy) >= 0.5) THEN
	536	! wrf_grid%lmsk(ix,iy)=1.
	537	! wrf_grid%lter(ix,iy)=1.
	538	! wrf_grid%llic(ix,iy)=0.
	539	! wrf_grid%loce(ix,iy)=0.
	540	! wrf_grid%lsic(ix,iy)=0.
	541	! ELSE
	542	! wrf_grid%lmsk(ix,iy)=0.
	543	! wrf_grid%lter(ix,iy)=0.
	544	! wrf_grid%llic(ix,iy)=0.
	545	! wrf_grid%loce(ix,iy)=1.
	546	! wrf_grid%lsic(ix,iy)=0.
	547	! END IF
	548	! END DO
	549	! END DO
	550	! End of Pseudo WRF-initialization
	551	!!!!!
	552
	553	! Value is given at the first time-step, not any more, since it will be continuously
	554	! used by LMDZ
	555	! l. Fita, LMD. January 2014
	556	! This should be in this way, but, because we can not use any scheme DZS=0., so
	557	! direct met_em values (Registry modified) are used. soil layers are in [cm]!
	558	! DO iz=1,SIZE(wrf_grid%soil_layers(ip,:,jp))
	559	! wrf_grid%lwsol = wrf_grid%lwsol + wrf_grid%sm(:,iz,:)* &
	560	! wrf_grid%soil_layers(:,iz,:)/100.
	561	! END DO
	562	DO iz=1,4
	563	wrf_grid%ltksoil(:,iz,:)=wrf_grid%st(:,1,:)
	564	END DO
	565	!
	566	! Only using SST on that ocean grid points
	567	DO ix=1,ddimx
	568	DO iy=1,ddimy
	569	IF (wrf_grid%loce(ix,iy) > 0.) THEN
	570	! Full SST grid points must have a SST value from the wrfinput (WRF 1/0 landmask)
	571	! but that LMDZ points with fractions of land/mask could not. Using the closest one
	572	! if it is an isolated grid point, using TSK instead
	573	IF (wrf_grid%sst(ix,iy) < 200.15) THEN
	574	wrf_grid%ltksoil(ix,3,iy) = -9.
	575	IF ( (ix > 1).AND.(ix < ddimx) .AND. (iy > 1).AND.(iy < ddimy) ) THEN
	576	DO ixx=-1,1
	577	DO iyy=-1,1
	578	! PRINT *,ixx,iyy,wrf_grid%sst(ix+ixx,iy+iyy)
	579	IF (wrf_grid%sst(ix+ixx,iy+iyy) > 200.15) THEN
	580	wrf_grid%ltksoil(ix,3,iy) = wrf_grid%sst(ix+ixx,iy+iyy)
	581	EXIT
	582	END IF
	583	END DO
	584	IF ( wrf_grid%ltksoil(ix,3,iy) > 200.15) EXIT
	585	END DO
	586	IF ( wrf_grid%ltksoil(ix,3,iy) == -9.) wrf_grid%ltksoil(ix,3,iy) = &
	587	wrf_grid%tsk(ix,iy)
	588	ELSE
	589	wrf_grid%ltksoil(ix,3,iy) = wrf_grid%tsk(ix,iy)
	590	END IF
	591	ELSE
	592	wrf_grid%ltksoil(ix,3,iy) = wrf_grid%sst(ix,iy)
	593	END IF
	594	IF ( (wrf_grid%ltksoil(ix,3,iy) < 200.15) .OR. &
	595	(wrf_grid%ltksoil(ix,3,iy) > 370.) ) THEN
	596	PRINT *,TRIM(errmsg)
	597	WRITE(wrf_err_message,*) ' ' // TRIM(subname) // &
	598	' SST inconsistency at point: ', ix, ', ', iy,' with loce: ', &
	599	wrf_grid%loce(ix,iy),' has a tsoil(oce)= ', &
	600	wrf_grid%ltksoil(ix,3,iy),' sst: ',wrf_grid%sst(ix,iy),' tsk: ', &
	601	wrf_grid%tsk(ix,iy)
	602	CALL wrf_error_fatal(TRIM(wrf_err_message))
	603	END IF
	604	END IF
	605	END DO
	606	END DO
	607
	608	!! Soil temperature, humidity, water-content...
	609	! This should go on the namelist... later!
	610	ORCHIDEElevels = (/ 0.098, 0.1955, 0.5865, 1.3685, 2.9326, 6.0606, 12.3167, &
	611	24.8289, 49.8534, 99.9022, 200. /)
	612
	613	datasetname = soilcat
	614	tablename = 'SOILPARM.TBL'
	615
	616	PRINT *,' Lluis reading dataset "' // TRIM(datasetname) // '" from "' // &
	617	TRIM(tablename) // '"'
	618
	619	CALL table_characteristics(datasetname,tablename,begdataset,Nvars,Nrs,Ncs)
	620
	621	! L. Fita, LMD. Feburary 2014
	622	! Again (as in the landuses case), the files content more categeories that the really ones used by the model...
	623	Nsoillayers=SIZE(wrf_grid%soil_layers(ip,:,jp))
	624	Nwrfsoiltypes=SIZE(wrf_grid%soilctop(ip,:,jp))
	625	IF (Nrs > Nwrfsoiltypes) THEN
	626	PRINT *,TRIM(wrnmsg)
	627	PRINT *," '" // TRIM(datasetname)//"' from '"//TRIM(tablename)//"' has ", &
	628	Nrs,' soil layers, meanwhile WRF data has ',Nwrfsoiltypes
	629	PRINT *,' using only the first ', Nwrfsoiltypes, ' of the dataset !!!!'
	630	Nrs=Nsoillayers
	631	ELSE IF (Nrs < Nwrfsoiltypes) THEN
	632	PRINT *,TRIM(errmsg)
	633	PRINT *," '" // TRIM(datasetname)//"' from '"//TRIM(tablename)//"' has ", &
	634	Nrs,' soil layers, meanwhile WRF data has ',Nwrfsoiltypes
	635	WRITE(wrf_err_message,*)' wrong the dataset !!!!'
	636	CALL wrf_error_fatal(TRIM(wrf_err_message))
	637	END IF
	638
	639	IF (ALLOCATED(soilvals)) DEALLOCATE(soilvals)
	640	ALLOCATE(soilvals(Nvars,Nrs,Ncs))
	641
	642	IF (ALLOCATED(soilvalns)) DEALLOCATE(soilvalns)
	643	ALLOCATE(soilvalns(Nrs))
	644
	645	CALL read_table(datasetname, tablename, begdataset, Nvars, Nrs, Ncs, &
	646	soilvals, soilvalns)
	647
	648	DO ix=1,ddimx
	649	DO iy=1,ddimy
	650	!L. Fita, LMD. January 2014 Soil layers are in cm!
	651
	652
	653	! NOTE: wrf_dx, wrf_dy should be matrix !!!
	654	CALL compute_TOTsoil_mass_water_values(datasetname, Nvars, Nrs, Ncs, &
	655	soilvals, soilvalns, ix, iy, Nrs, wrf_grid%soilctop(ix,:,iy), &
	656	wrf_grid%soilcbot(ix,:,iy), Nsoillayers, wrf_grid%sm(ix,:,iy), &
	657	wrf_grid%soil_layers(ix,Nsoillayers:1:-1,iy)/100., wrf_dx, wrf_dy, &
	658	soilTmass, wrf_grid%lwsol(ix,iy))
	659
	660	wrf_grid%lwsol(ix,iy) = wrf_grid%lwsol(ix,iy)/(wrf_dxwrf_dy1000.)
	661
	662	CALL interpolate_layers(SIZE(wrf_grid%soil_layers(ix,:,iy)), &
	663	wrf_grid%soil_layers(ix,SIZE(wrf_grid%soil_layers(ix,:,iy)):1:-1,iy), &
	664	wrf_grid%st(ix,:,iy),11, ORCHIDEElevels,wrf_grid%lotter(ix,:,iy))
	665	END DO
	666	END DO
	667
	668	wrf_grid%lotlic=wrf_grid%lotter
	669	wrf_grid%lotoce=wrf_grid%lotter
	670	wrf_grid%lotsic=wrf_grid%lotter
	671
	672	DEALLOCATE(soilvals)
	673	DEALLOCATE(soilvalns)
	674
	675	! LMDZ does not distinguish among species, use it as a whole
	676	wrf_grid%lclwcon=SUM(wrf_MOIST,4)
	677	DO ix=wims,wime
	678	DO iz=wkms,wkme
	679	DO iy=wjms,wjme
	680	IF (wrf_grid%lclwcon(ix,iz,iy) /= 0.) THEN
	681	wrf_grid%lratqs=(wrf_MOIST(ix,iz,iy,wrf_qrid)+ &
	682	wrf_MOIST(ix,iz,iy,wrf_qsid)+wrf_MOIST(ix,iz,iy,wrf_qiid)+ &
	683	wrf_MOIST(ix,iz,iy,wrf_qhid)+wrf_MOIST(ix,iz,iy,wrf_qgid))/ &
	684	wrf_grid%lclwcon(ix,iz,iy)
	685	END IF
	686	END DO
	687	END DO
	688	END DO
	689
	690	! wrf_grid%znt = wrf_grid%z0
	691	! wrf_grid%albedo = wrf_grid%albbck
	692	wrf_grid%lbils=0.
	693	wrf_grid%lfder=0.
	694	wrf_grid%lrnebcon=0.
	695	DO ix=1,ddimx
	696	DO iy=1,ddimy
	697	ksoils = (/ wrf_grid%lter(ix,iy), wrf_grid%llic(ix,iy), &
	698	wrf_grid%loce(ix,iy), wrf_grid%lsic(ix,iy) /)
	699	ksoilmax = MAXVAL(ksoils)
	700	DO iz = 1,4
	701	IF (ksoilmax == ksoils(iz)) wrf_grid%lnat(ix,iy) = iz
	702	END DO
	703	END DO
	704	END DO
	705	wrf_grid%lrads = wrf_grid%swdown
	706
	707	wrf_grid%lcontrol(1) = wtime_step
	708	defilen = 'config.def'
	709	defvarn = 'nbapp_rad'
	710	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	711	wrf_grid%lcontrol(6) = defint
	712	wrf_grid%lcontrol(2) = NINT (86400./wtime_step/ FLOAT(defint) )
	713	!L. Fita, LMD. February 2014: current value of atmospheric CO2 should be
	714	!! co2_ppm0 + delta_co2_ppm, from 'carbon_cycle_mod.F90'
	715	defvarn = 'co2_ppm'
	716	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	717	wrf_grid%lcontrol(3) = defint
	718	!!c co2_ppm0 : initial value of atmospheric CO2
	719	wrf_grid%lcontrol(16) = defint
	720	defvarn = 'solaire'
	721	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	722	wrf_grid%lcontrol(4) = defreal
	723	wrf_grid%lcontrol(5) = wrf_l_con
	724
	725	! Default value cycle_diurne = .T.
	726	wrf_grid%lcontrol(7) = 1.
	727	defvarn = 'soil_model'
	728	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	729	IF(defbool) wrf_grid%lcontrol(8) = 1.
	730	defilen = 'physiq.def'
	731	defvarn = 'new_oliq'
	732	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	733	IF(defbool) wrf_grid%lcontrol(9) = 1.
	734	defvarn = 'ok_orodr'
	735	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	736	IF(defbool) wrf_grid%lcontrol(10) = 1.
	737	defvarn = 'ok_orolf'
	738	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
	739	IF(defbool) wrf_grid%lcontrol(11) = 1.
	740
	741	wrf_grid%lcontrol(13) = dayref + nday
	742	wrf_grid%lcontrol(14) = annee_ref
	743	wrf_grid%lcontrol(15) = wrf_grid%itimestep
	744
	745	! Let's split rugosity by Ksoils
	746	!! Static roughness length
	747	wrf_grid%lrug=wrf_grid%znt
	748
	749	!! Dynamic roughness length (It should work if we use WRF schemes...)
	750	! wrf_grid%lrug=wrf_grid%znt
	751	! First split by Ksoil types.
	752	! NOTE: in somehow we should read again the table 'VEGPARM.TBL' and recompute the
	753	! roughness length in case that appear a new soil type. (but what left below a
	754	! melted glacier?)
	755	!
	756	LMDZvarmethod='prod'
	757	wrf_grid%lrugksoil(:,1,:) = wrf_grid%znt
	758
	759	CALL lmdz_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,wrf_grid%snow, &
	760	LMDZvarmethod,wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
	761	wrf_grid%lsnowksoil(:,1,:),wrf_grid%lsnowksoil(:,2,:), &
	762	wrf_grid%lsnowksoil(:,3,:),wrf_grid%lsnowksoil(:,4,:))
	763	! All the iced soil types have the same roughtness length 0.001
	764	wrf_grid%lrugksoil(:,2,:)=0.001
	765	wrf_grid%lrugksoil(:,3,:)=0.0001
	766	wrf_grid%lrugksoil(:,4,:)=0.001
	767
	768	! Recomputig again the total roughness length
	769	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
	770	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
	771	wrf_grid%lrugksoil(:,1,:),wrf_grid%lrugksoil(:,2,:), &
	772	wrf_grid%lrugksoil(:,3,:),wrf_grid%lrugksoil(:,4,:),wrf_grid%lrug)
	773
	774	! Air surface moisture for each soil type. First order?
	775	DO iz=1,4
	776	wrf_grid%lqksoil(:,iz,:)=wrf_grid%q2
	777	END DO
	778
	779	! Let's split albedo by Ksoils
	780	wrf_grid%lalbksoil=0.
	781
	782	DO iz=1,4
	783	!! Static roughness length (changing on time in wrflwbdy_d01)
	784	wrf_grid%lalbksoil(:,iz,:)=wrf_grid%albedo
	785	!! Dynamic roughness length (It should work if we use WRF schemes...)
	786	! wrf_grid%lalbksoil(:,iz,:)=wrf_grid%albedo
	787	END DO
	788	! First split by Ksoil types.
	789	! NOTE: in somehow we should read again the table 'VEGPARM.TBL' and recompute the
	790	! albedo in case that appear a new soil type. (but what left below a melted
	791	! glacier?)
	792	!
	793	! All the iced soil types have the same albedo
	794	! For snow/ice there is zmin=0.55, zmax=0.70, tacking the mean value: 0.625
	795	wrf_grid%lalbksoil(:,2,:)=0.625
	796	wrf_grid%lalbksoil(:,3,:)=0.08
	797	wrf_grid%lalbksoil(:,4,:)=0.625
	798
	799	! Recomputig again the total roughness length
	800	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
	801	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
	802	wrf_grid%lalbksoil(:,1,:),wrf_grid%lalbksoil(:,2,:), &
	803	wrf_grid%lalbksoil(:,3,:),wrf_grid%lalbksoil(:,4,:),wrf_grid%albedo)
	804	! Using the same values for the LW albedo
	805	wrf_grid%llwalbksoil = wrf_grid%lalbksoil
	806
	807	! Orography values just beacuse there is no other way right now (they should come
	808	! from geogrid.exe)
	809	wrf_grid%lzmea = wrf_grid%ht
	810	wrf_grid%lzstd = wrf_grid%topostdv
	811	wrf_grid%lzsig = wrf_grid%slope
	812	wrf_grid%lzgam = SQRT((wrf_grid%toposlpx2. + wrf_grid%toposlpy2.))/ &
	813	(1. + wrf_grid%slope)
	814	wrf_grid%lzthe = ACOS(wrf_grid%toposlpx/SQRT(wrf_grid%toposlpx**2. + &
	815	wrf_grid%toposlpy**2.))
	816	wrf_grid%lzpic = wrf_grid%ht*1.10
	817	wrf_grid%lzval = wrf_grid%ht*0.90
	818	wrf_grid%lzrugsrel = wrf_grid%lrugKsoil(:,1,:)
	819
	820	! R = 8.3144621 ([J mol/K] http://en.wikipedia.org/wiki/Gas_constant)
	821	! mu = 28.97 (Molecular Mass of Air [kg/kmol] from http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html)
	822	! R / mu =? r_d
	823
	824	! lmdz_Area = dLon * dLat * cos(lat) * Earth_rad
	825
	826	END IF xtime0
	827
	828	! Grids do not start/end at the same index
	829	wrf_area = 0.
	830	wrf_cu = 0.
	831	wrf_cv = 0.
	832	DO iy=1,ddimy
	833	DO ix=1,ddimx
	834	wrf_area(ix,iy) = wrf_dx*wrf_dy/wrf_mapft(ix,iy)
	835	wrf_cu(ix,iy) = wrf_dx/wrf_mapfu(ix,iy)
	836	wrf_cv(ix,iy) = wrf_dy/wrf_mapfv(ix,iy)
	837	END DO
	838	END DO
	839
	840	! cu = dLon * Earth_rad
	841	! cv = dLat * Earth_rad * cos(Lat)
	842
	843	! lmdz_cu(1) = cu(1)
	844	! lmdz_cv(1) = cv(1)
	845	! DO j=2,dimy
	846	! DO i=1,dimx
	847	! lmdz_cu((j-2)dimx+1+i) = cu((j-1)(dimx+1)+i)
	848	! lmdz_cv((j-2)dimx+1+i) = cv((j-1)(dimx+1)+i)
	849	! ENDDO
	850	! ENDDO
	851	! lmdz_cu(dim_ngridmx) = cu(dimxy+1)
	852	! lmdz_cv(dim_ngridmx) = cv(dimxy-dimx)
	853
	854	AreaValues = RESHAPE(wrf_area, (/ ddimx*ddimy/), ORDER=(/1/))
	855	CUValues = RESHAPE(wrf_cu, (/ ddimx*ddimy/), ORDER=(/1/))
	856	CVValues = RESHAPE(wrf_cv, (/ ddimx*ddimy/), ORDER=(/1/))
	857
	858	CALL lmdz_vars_init(ddimxy)
	859
	860	CALL load_lmdz(wrf_grid,ddimx,ddimy,dimz,wbdyw,ddimxy,4,11,rlon,rlat,zmasq, &
	861	pctsrf,ftsol,ftsoil_rst,qsurf_rst,qsol_rst,falb1,falb2,evap_rst,snow_rst, &
	862	radsol,solsw,sollw,fder_rst,rain_fall,snow_fall,agesno_rst,zmea,zstd,zgam, &
	863	zthe,zpic,zval,rugoro,rugos_rst,restart_runoff,zmax0,f0,wake_s,wake_cstar, &
	864	wake_pe,wake_fip,sst,albedo,t_ancien,q_ancien,u_ancien,v_ancien,clwcon, &
	865	rnebcon,ratqs,ema_work1,ema_work2,wake_deltat,wake_deltaq,fm_therm, &
	866	entr_therm,detr_therm,phy_glo)
	867
	868	CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst, &
	869	rugos_rst, agesno_rst, ftsoil_rst)
	870	CALL fonte_neige_init(restart_runoff)
	871
	872	DO ix=1,ddimx
	873	DO iy=1,ddimy
	874	ixy=ddimx*(iy-1)+ix
	875	IF (wrf_grid%lter(ix,iy) /= pctsrf(ixy,1)) PRINT *,' ' // &
	876	TRIM(subname) // ' ter DOES NOT COINCIDE! ', ix,', ',iy,' lter: ', &
	877	wrf_grid%lter(ix,iy),' pctsrf: ',pctsrf(ixy,1)
	878	IF (wrf_grid%llic(ix,iy) /= pctsrf(ixy,2)) PRINT *,' ' // &
	879	TRIM(subname) // ' lic DOES NOT COINCIDE! ', ix,', ',iy,' llic: ', &
	880	wrf_grid%llic(ix,iy),' pctsrf: ',pctsrf(ixy,2)
	881	IF (wrf_grid%loce(ix,iy) /= pctsrf(ixy,3)) PRINT *,' ' // &
	882	TRIM(subname) // ' oce DOES NOT COINCIDE! ', ix,', ',iy,' loce: ', &
	883	wrf_grid%loce(ix,iy),' pctsrf: ',pctsrf(ixy,3)
	884	IF (wrf_grid%lsic(ix,iy) /= pctsrf(ixy,4)) PRINT *,' ' // &
	885	TRIM(subname) // ' sic DOES NOT COINCIDE! ', ix,', ',iy,' lsic: ', &
	886	wrf_grid%lsic(ix,iy),' pctsrf: ',pctsrf(ixy,4)
	887	END DO
	888	END DO
	889
	890	CALL NOread_limit_init(ddimxy, klon_glo)
	891
	892	sst_glo = sst(1:ddimxy)
	893	! In ocean_forced_mod is used as tsurf_limit = sst(knindex(1:knon))
	894	tsurf_limit = ftsol(:,3)
	895
	896	! In surf_land_bucket_mod is used as z0_new=rugos(knindex(1:knon),1) [z0_new == z0_limit]
	897	z0_limit = rugos_rst(:,1)
	898	alb_limit = albedo
	899
	900	rug_glo = rugos_rst(:,is_ter)
	901	pct_glo(klon_glo,:) = pctsrf(klon_glo,:)
	902
	903	! It is not clear how to deal with 'limit.nc' times...
	904	! temps=wjulday
	905
	906	CALL iniphysiq(ddimx*ddimy, dimz, 86400., wjulday, wtime_step, &
	907	& LatValues, LonValues, wrf_area, CUValues, CVValues, 1./reradius, g, r_d, cp, &
	908	& iflag_phys)
	909
	910	! NO NEED to do this, using wrfinput_d[nn], wrfbdy_d[nn] and wrflowinput_d[nn]
	911	! Checking for 'startphy.nc' presence
	912	!! INQUIRE(FILE='startphy.nc', EXIST=existsfile)
	913	!! IF (.NOT.existsfile) THEN
	914	!! WRITE(wrf_err_message,*)" file 'startphy.nc' with initial conditions " // &
	915	!! & ' for LMDZphys does not exist!!!!'
	916	!! CALL wrf_error_fatal(TRIM(wrf_err_message))
	917	!! END IF
	918
	919	CALL init_ovars_lmdz_NOmodule(ddimxy,dimz,nbsrf)
	920
	921	PRINT *,' Lluis: END of firststep !!!!!'
	922	END IF firststep
	923
	924	ip=icheck_p
	925	jp=jcheck_p
	926	! WRF alarms
	927	! I can not print them !!!!
	928	!!
	929	! wrf_grid%alarms => wrf_alarms
	930	! PRINT *,' Lluis WRF alarms: time-step', wrf_grid%itimestep,' _______'
	931	! DO ix=1,3*INT(MAX_WRF_ALARMS/3),3
	932	! PRINT *,(ix+iy,': ',wrf_alarms(ix+iy), iy=0,2)
	933	! END DO
	934	! PRINT ,(iy,': ',wrf_alarms(iy), iy=3INT(MAX_WRF_ALARMS/3)+1, MAX_WRF_ALARMS)
	935
	936	! What to do with this?
	937	!!! IF ( wrf_restart_alarm ) lmdz_lafin = .TRUE.
	938
	939	INQUIRE(FILE='traceur.def', EXIST=existsfile)
	940	IF (existsfile) THEN
	941	DO funit=10,100
	942	INQUIRE(UNIT=funit, OPENED=is_used)
	943	IF (.NOT. is_used) EXIT
	944	END DO
	945	OPEN(UNIT=10, FILE='traceur.def', STATUS='OLD')
	946	READ(10,*)lmdzmixingratios
	947	CLOSE(UNIT=10)
	948	ELSE
	949
	950	WRITE(message, *)errmsg
	951	CALL wrf_debug(75, message)
	952	WRITE(wrf_err_message,*)' required file traceur.def not found !!'
	953	CALL wrf_error_fatal(TRIM(wrf_err_message))
	954	END IF
	955
	956	! Regarding if it is a lowbdy input time
	957	!!
	958	PRINT *,' ',wrf_grid%itimestep,' wrf_islowbdyin: ',wrf_islowbdyin
	959	! IF ((wrf_islowbdyin) .AND. .NOT.wrf_restart_alarm) THEN
	960	IF ((wrf_islowbdyin)) THEN
	961	PRINT *,' LOWBDY time-step!!!!!!! '
	962	PRINT *,' before: wsst ', wrf_grid%sst(ip,jp),' lsst: ',sst(lp)
	963	CALL get_lowbdy(ddimx, ddimy, dimz, wbdyw, ddimxy, 4, 11, klon_glo, &
	964	wrf_grid, pctsrf, ftsol, rugos_rst, albedo, sst, tsurf_limit, z0_limit, &
	965	alb_limit, rug_glo, pct_glo)
	966	PRINT *,' after: wsst ', wrf_grid%sst(ip,jp),' lsst: ',sst(lp)
	967	! STOP
	968	END IF
	969
	970	IF (ALLOCATED(MixingRatioValues)) DEALLOCATE(MixingRatioValues)
	971	ALLOCATE(MixingRatioValues(ddimxy,dimz,lmdzmixingratios), STAT=ierr)
	972	IF (ierr /= 0) THEN
	973	PRINT *,errmsg
	974	PRINT *,' ' // subname // ': Problem allocating WRF moisture array!', ierr
	975	END IF
	976
	977	IF (ALLOCATED(MixingRatiotendValues)) DEALLOCATE(MixingRatiotendValues)
	978	ALLOCATE(MixingRatiotendValues(ddimxy,dimz,lmdzmixingratios))
	979
	980	! Transforming from WRF theta temperature to LMDZ temperature (in K)
	981	!!
	982	wrf_P = wrf_baseP + wrf_perP
	983	wrf_geopot = wrf_basegeopot + wrf_pergeopot
	984	wrf_sfcgeopot = wrf_grid%ht*g
	985	! de-staggering winds
	986	!!
	987
	988	wrfUdestagg(wims:wime-1,:,wjms:wjme-1) = (wrf_U(wims:wime-1,:,wjms:wjme-1) + &
	989	wrf_U(wims+1:wime,:,wjms:wjme-1))/2.
	990	wrfVdestagg(wims:wime-1,:,wjms:wjme-1) = (wrf_V(wims:wime-1,:,wjms:wjme-1) + &
	991	wrf_V(wims:wime-1,:,wjms+1:wjme))/2.
	992	wrfUtenddestagg(wims:wime-1,:,wjms:wjme-1) = &
	993	(wrf_Utend(wims:wime-1,:,wjms:wjme-1) + wrf_Utend(wims+1:wime,:,wjms:wjme-1))/2.
	994	wrfVtenddestagg(wims:wime-1,:,wjms:wjme-1) = &
	995	(wrf_Vtend(wims:wime-1,:,wjms:wjme-1) + wrf_Vtend(wims:wime-1,:,wjms+1:wjme))/2.
	996
	997	!! NO de-stagger in z
	998	!! wrfWValues = (wrf_w(:,wkms:wkme-1,:)+wrf_w(:,wkms+1:wkme,:))/2.
	999	!! wrfPERMASSValues = (wrf_permass(:,wkms:wkme-1,:)+wrf_permass(:,wkms+1:wkme,:))/2.
	1000	!! DO k=wkms,wkme
	1001	!! wrfMASSValues(:,k,:)=-(wrfPERMASSVAlues(:,k,:)+wrf_basemass)*wrf_dfulleta(k)/g
	1002	!! END DO
	1003	wrfWValues=wrf_w
	1004	DO iz=wkms,wkme
	1005	wrfMASSValues(:,iz,:) = -(wrf_permass+wrf_basemass)*wrf_dfulleta(iz)/g
	1006	END DO
	1007	wrfZFluxMassValues=wrfWValues*wrfMASSValues
	1008
	1009	!! CALL theta_to_temp(THETAVALUES=wrf_T, PRESSVALUES=wrf_P, TEMPVALUES=wrf_Temp &
	1010	!! ! Dimension arguments
	1011	!! & ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
	1012	!! & ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
	1013	!! & ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
	1014	!! & ,I_START=wi_start,I_END=wi_end &
	1015	!! & ,J_START=wj_start,J_END=wj_end &
	1016	!! & ,KTS=wkts, KTE=wkte &
	1017	!! & ,NUM_TILES=wnum_tiles &
	1018	!! & )
	1019
	1020	CALL theta_to_temp(wrf_T, wrf_P, wrf_Temp &
	1021	! Dimension arguments
	1022	& ,wids,wide, wjds,wjde, wkds,wkde &
	1023	& ,wims,wime, wjms,wjme, wkms,wkme &
	1024	& ,wips,wipe, wjps,wjpe, wkps,wkpe &
	1025	& ,wi_start,wi_end &
	1026	& ,wj_start,wj_end &
	1027	& ,wkts, wkte &
	1028	& ,wnum_tiles &
	1029	& )
	1030
	1031	!! CALL theta_to_temp(THETAVALUES=wrf_Ttend, PRESSVALUES=wrf_P, &
	1032	!! & TEMPVALUES=wrf_Temptend &
	1033	!! ! Dimension arguments
	1034	!! & ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
	1035	!! & ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
	1036	!! & ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
	1037	!! & ,I_START=wi_start,I_END=wi_end &
	1038	!! & ,J_START=wj_start,J_END=wj_end &
	1039	!! & ,KTS=wkts, KTE=wkte &
	1040	!! & ,NUM_TILES=wnum_tiles &
	1041	!! & )
	1042
	1043	CALL theta_to_temp(wrf_Ttend, wrf_P, wrf_Temptend &
	1044	! Dimension arguments
	1045	& ,wids,wide, wjds,wjde, wkds,wkde &
	1046	& ,wims,wime, wjms,wjme, wkms,wkme &
	1047	& ,wips,wipe, wjps,wjpe, wkps,wkpe &
	1048	& ,wi_start,wi_end &
	1049	& ,wj_start,wj_end &
	1050	& ,wkts, wkte &
	1051	& ,wnum_tiles &
	1052	& )
	1053
	1054	DO ix=1,ddimx
	1055	DO iy=1,ddimy
	1056	sfcGeopotValues(ddimx(iy-1)+ix) = wrf_grid%ht(ix,iy)g
	1057	END DO
	1058	END DO
	1059
	1060	CALL mat_vect(wrf_Temp, dimx, dimy, dimz, wbdyw, TValues)
	1061	CALL mat_vect(wrf_Temptend, dimx, dimy, dimz, wbdyw, TtendValues)
	1062	CALL mat_vect(wrfUdestagg, dimx, dimy, dimz, wbdyw, UValues)
	1063	CALL mat_vect(wrfUtenddestagg, dimx, dimy, dimz, wbdyw, UtendValues)
	1064	CALL mat_vect(wrfVdestagg, dimx, dimy, dimz, wbdyw, VValues)
	1065	CALL mat_vect(wrfVtenddestagg, dimx, dimy, dimz, wbdyw, VtendValues)
	1066	! CALL mat_vect(wrf_U, dimx, dimy, dimz, wbdyw, UValues)
	1067	! CALL mat_vect(wrf_Utend, dimx, dimy, dimz, wbdyw, UtendValues)
	1068	! CALL mat_vect(wrf_V, dimx, dimy, dimz, wbdyw, VValues)
	1069	! CALL mat_vect(wrf_Vtend, dimx, dimy, dimz, wbdyw, VtendValues)
	1070	CALL mat_vect(wrf_P, dimx, dimy, dimz, wbdyw, HalfPressValues)
	1071	CALL mat_vect_zstagg(wrf_fullpres, dimx, dimy, dimz+1, wbdyw, FullPressValues)
	1072	CALL mat_vect_zstagg(wrf_geopot, dimx, dimy, dimz+1, wbdyw, GeopotValues)
	1073	CALL mat_vect(wrfZFluxMassValues, dimx, dimy, dimz, wbdyw, WFluxMassValues)
	1074
	1075	CALL WRFmoist_LMDZmoist(wrf_MOIST, dimx, dimy, dimz, wnum3dm, wbdyw, wrf_qvid, &
	1076	& wrf_qcid, wrf_qrid, lmdzmixingratios, MixingRatioValues)
	1077
	1078	CALL WRFmoist_LMDZmoist(wrf_MOISTtend, dimx, dimy, dimz, wnum3dm, wbdyw, wrf_qvid,&
	1079	& wrf_qcid, wrf_qrid, lmdzmixingratios, MixingRatiotendValues)
	1080
	1081	preff=100000.
	1082	CALL eta_to_pressure(wrf_halfeta, preff, wrf_Ptop, dimz, lmdzoutP)
	1083
	1084	! We need to provide values to these variables in order to make it works
	1085	presnivs=lmdzoutP
	1086
	1087	! Lluis Testing water species
	1088	MixingRatioValues(:,:,2)=0.
	1089
	1090	IF (wrf_dbg >= 200) THEN
	1091
	1092	!! Values at a given LMDZ matrix point
	1093	lmdzp=lp
	1094	jlmdz=INT(lmdzp/ddimx) + 1
	1095	ilmdz=lmdzp-(jlmdz-1)*ddimx
	1096	!! Values at a given WRF 2D point
	1097	! ilmdz=ip
	1098	! jlmdz=jp
	1099	! lmdzp=ddimx*(jlmdz-1)+ilmdz
	1100
	1101	PRINT *,' WRF+LMDZ lmdzp: ',lmdzp,' : ',ilmdz,jlmdz
	1102	PRINT *,'k wrf_T wrf_Temp TValues_______________'
	1103	DO iz=1,dimz
	1104	PRINT *,iz, wrf_T(ilmdz,iz,jlmdz), wrf_Temp(ilmdz,iz,jlmdz), TValues(lmdzp,iz)
	1105	END DO
	1106	PRINT *,'k wrf_U Uvalues wrf_V VValues _______________'
	1107	DO iz=1,dimz
	1108	PRINT *,iz, wrf_U(ilmdz,iz,jlmdz), UValues(lmdzp,iz), wrf_V(ilmdz,iz,jlmdz), &
	1109	VValues(lmdzp,iz)
	1110	END DO
	1111	PRINT *,'k wrf_fullpres FUllPressValues wrf_sfcgeopot hgt sfcGeopotValues______'
	1112	DO iz=1,dimz+1
	1113	PRINT *,iz,wrf_fullpres(ilmdz,iz,jlmdz), FullPressValues(lmdzp,iz), &
	1114	wrf_sfcgeopot(ilmdz,jlmdz), wrf_grid%ht(ilmdz,jlmdz), sfcGeopotValues(lmdzp)
	1115	END DO
	1116	PRINT *,'k wrf_halfpress HalfPressValues wrf_Geopot GeopotValues lmdz outP ' //&
	1117	'_______________'
	1118	DO iz=1,dimz
	1119	PRINT *,iz, wrf_P(ilmdz,iz,jlmdz), HalfPressValues(lmdzp,iz), &
	1120	wrf_geopot(ilmdz,iz,jlmdz), GeopotValues(lmdzp,iz), lmdzoutP(iz)
	1121	END DO
	1122	PRINT *,'k wrf_MOIST MixingRatioValues _______________'
	1123	DO iz=1,dimz
	1124	PRINT *,iz,wrf_MOIST(ilmdz,iz,jlmdz,wrf_qvid), MixingRatioValues(lmdzp,iz,1),&
	1125	wrf_MOIST(ix,iz,jlmdz,wrf_qcid), MixingRatioValues(lmdzp,iz,2)
	1126	END DO
	1127	PRINT *,'k wrf_mass wrf_W wrf_dfulleta wrfMASSValues wrfZFluxMassValues ' // &
	1128	'WFluxMassValues ____________'
	1129	DO iz=1,dimz
	1130	PRINT *,iz,wrf_permass(ilmdz,jlmdz)+wrf_basemass(ilmdz,jlmdz), &
	1131	wrf_W(ilmdz,iz,jlmdz), wrf_dfulleta(iz), wrfMASSValues(ilmdz,iz,jlmdz), &
	1132	wrfZFluxMassValues(ilmdz,iz,jlmdz), WFluxMassValues(lmdzp,iz)
	1133	END DO
	1134	END IF
	1135
	1136	ip=icheck_p
	1137	jp=jcheck_p
	1138
	1139	!! L. Fita, July 2013.
	1140	!! Misteriously 'physiq' suborutine from LMDZ can not be called using keys
	1141	!!
	1142	!! CALL physiq(NLON=dimxy, NLEV=dimz, &
	1143	!! & DEBUT=lmdz_debut, LAFIN=lmdz_lafin, &
	1144	!! & JD_CUR=wjulday, JH_CUR=wgmt, PDTPHYS=wtime_step, &
	1145	!! & PAPRS=FullPressValues, PPLAY=HalfPressValues, &
	1146	!! & PPHI=GeopotValues, PPHIS=wrf_sfcgeopot, &
	1147	!! & PRESNIVS=lmdzoutP, CLESPHY0=phykeys, &
	1148	!! & U=UValues, V=VValues, T=TValues, QX=MixingRatioValues, &
	1149	!! & FLXMASS_W=WFluxMassValues, &
	1150	!! & D_U=UtendValues, D_V=VtendValues, D_T=TtendValues, &
	1151	!! & D_QX=MixingRatiotendValues, D_PS=wrf_psfctend, &
	1152	!! & DUDYN=lmdz_dudyn, PVTETA=lmdz_PVtheta)
	1153
	1154	!! PRINT *,' Lluis: going to physiq'
	1155
	1156	CALL physiq(ddimxy, dimz, &
	1157	& lmdz_debut, lmdz_lafin, &
	1158	& wjulday, wgmt, wtime_step, &
	1159	& FullPressValues, HalfPressValues, &
	1160	& GeopotValues, sfcGeopotValues, &
	1161	& lmdzoutP, phykeys, &
	1162	& UValues, VValues, TValues, MixingRatioValues, &
	1163	& WFluxMassValues, &
	1164	& UtendValues, VtendValues, TtendValues, &
	1165	& MixingRatiotendValues, psfctendValues, &
	1166	& lmdz_dudyn, lmdz_PVtheta)
	1167
	1168	ip=icheck_p
	1169	jp=jcheck_p
	1170	! Lluis Testing water species
	1171	!! MixingRatiotendValues(:,:,1)=1000.*MixingRatiotendValues(:,:,1)
	1172	!! MixingRatiotendValues(:,:,1)=MixingRatiotendValues(:,:,1)+MixingRatiotendValues(:,:,2)
	1173	MixingRatiotendValues(:,:,2)=0.
	1174
	1175	!! PRINT *,' Lluis: outgoing from physiq'
	1176
	1177	!! Tvalues, [U/V/T]Values, are not given back just inputs
	1178	CALL vect_mat(TtendValues, dimx, dimy, dimz, wbdyw, wrf_Temptend)
	1179	CALL vect_mat(UtendValues, dimx, dimy, dimz, wbdyw, wrfUtenddestagg)
	1180	CALL vect_mat(VtendValues, dimx, dimy, dimz, wbdyw, wrfVtenddestagg)
	1181	! CALL vect_mat(UtendValues, dimx, dimy, dimz, wbdyw, wrf_Utend)
	1182	! CALL vect_mat(VtendValues, dimx, dimy, dimz, wbdyw, wrf_Vtend)
	1183	!! CALL vect_mat(TValues, dimx, dimy, dimz, wbdyw, wrf_Temp)
	1184	!! CALL vect_mat(UValues, dimx, dimy, dimz, wbdyw, wrf_U)
	1185	!! CALL vect_mat(VValues, dimx, dimy, dimz, wbdyw, wrf_V)
	1186	!! CALL vect_mat(HalfPressValues, dimx, dimy, dimz, wbdyw, wrf_P)
	1187	!! CALL vect_mat(FullPressValues, dimx, dimy, dimz, wbdyw, wrf_fullpres)
	1188	!! CALL vect_mat(GeopotValues, dimx, dimy, dimz, wbdyw, wrf_geopot)
	1189
	1190	DO iy=1,ddimy
	1191	DO ix=1,ddimx
	1192	wrf_psfctend(ix,iy) = PsfctendValues(ddimx*(iy-1)+ix)
	1193	END DO
	1194	END DO
	1195
	1196	!! No changes are suffered not need to give it back
	1197	!! CALL temp_to_theta(TEMPVALUES=wrf_Temp, PRESSVALUES=wrf_P, THETAVALUES=wrf_T &
	1198	!! ! Dimension arguments
	1199	!! & ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
	1200	!! & ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
	1201	!! & ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
	1202	!! & ,I_START=wi_start,I_END=wi_end &
	1203	!! & ,J_START=wj_start,J_END=wj_end &
	1204	!! & ,KTS=wkts, KTE=wkte &
	1205	!! & ,NUM_TILES=wnum_tiles &
	1206	!! & )
	1207
	1208	CALL temp_to_theta(TEMPVALUES=wrf_Temptend, PRESSVALUES=wrf_P, &
	1209	& THETAVALUES=wrf_Ttend &
	1210	! Dimension arguments
	1211	& ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
	1212	& ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
	1213	& ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
	1214	& ,I_START=wi_start,I_END=wi_end &
	1215	& ,J_START=wj_start,J_END=wj_end &
	1216	& ,KTS=wkts, KTE=wkte &
	1217	& ,NUM_TILES=wnum_tiles &
	1218	& )
	1219
	1220	wrf_Ttend = wrf_Ttend + 300.
	1221
	1222	! Re-staggering wind tendencies
	1223	!!
	1224
	1225	wrf_Utend(wims,:,wjms:wjme-1) = wrfUtenddestagg(wims,:,wjms:wjme-1)
	1226	wrf_Vtend(wims:wime-1,:,wjms) = wrfVtenddestagg(wims:wime-1,:,wjms)
	1227	wrf_Utend(wims,:,wjme-1) = wrfUtenddestagg(wims,:,wjme-1)
	1228	wrf_Vtend(wime-1,:,wjms) = wrfVtenddestagg(wime-1,:,wjms)
	1229	wrf_Utend(wims+1:wime-1,:,wjms:wjme-1) = &
	1230	& (wrfUtenddestagg(wims:wime-2,:,wjms:wjme-1) + &
	1231	& wrfUtenddestagg(wims+1:wime-1,:,wjms:wjme-1))/2.
	1232	wrf_Vtend(wims:wime-1,:,wjms+1:wjme-1) = &
	1233	& (wrfVtenddestagg(wims:wime-1,:,wjms:wjme-2) + &
	1234	& wrfVtenddestagg(wims:wime-1,:,wjms+1:wjme-1))/2.
	1235	wrf_Utend(wime,:,wjms:wjme-1) = wrfUtenddestagg(wime-1,:,wjms:wjme-1)
	1236	wrf_Vtend(wims:wime-1,:,wjme) = wrfVtenddestagg(wims:wime-1,:,wjme-1)
	1237	wrf_Utend(wime,:,wjme) = wrfUtenddestagg(wime-1,:,wjme-1)
	1238	wrf_Vtend(wime,:,wjme) = wrfVtenddestagg(wime-1,:,wjme-1)
	1239
	1240	!! wrf_Ttend = wrf_Temptend
	1241
	1242	!! wrf_perP = wrf_P - wrf_baseP
	1243	!! wrf_pergeopot = wrf_geopot - wrf_basegeopot
	1244
	1245	!! CALL LMDZmoist_WRFmoist(MixingRatioValues, dimx, dimy, dimz, wnum3dm, wrf_qvid, &
	1246	!! & wrf_qcid, wrf_qrid, lmdzmixingratios, wrf_MOIST)
	1247
	1248	CALL LMDZmoist_WRFmoist(MixingRatiotendValues, dimx, dimy, dimz, wnum3dm, wbdyw, &
	1249	& wrf_qvid, wrf_qcid, wrf_qrid, lmdzmixingratios, wrf_MOISTtend)
	1250
	1251	IF (wrf_dbg >= 200) THEN
	1252	PRINT *,' WRF+LMDZ lmdzp: ',lmdzp
	1253	PRINT *,'k T tendencies _______________'
	1254	DO iz=1,dimz
	1255	PRINT *,iz, TtendValues(lmdzp,iz), wrf_Temptend(ilmdz,iz,jlmdz), &
	1256	wrf_Ttend(ilmdz,iz,jlmdz)
	1257	END DO
	1258	PRINT *,'k U tendencies V tendencies_______________'
	1259	DO iz=1,dimz
	1260	PRINT *,iz, UtendValues(lmdzp,iz), wrf_Utend(ilmdz,iz,jlmdz), &
	1261	VtendValues(lmdzp,iz), wrf_Vtend(ilmdz,iz,jlmdz)
	1262	END DO
	1263	PRINT *,'k MixingRatiotendValues wrf_MOISTtend _______________'
	1264	DO iz=1,dimz
	1265	PRINT *,iz, MixingRatiotendValues(lmdzp,iz,1), &
	1266	wrf_MOISTtend(ilmdz,iz,jlmdz,wrf_qvid), MixingRatiotendValues(lmdzp,iz,2), &
	1267	wrf_MOISTtend(ilmdz,iz,jlmdz,wrf_qcid)
	1268	END DO
	1269	PRINT *,'wrf_sfcgeopot psfc_tendencies sfcGeopotValues_______________'
	1270	PRINT *,wrf_sfcgeopot(ilmdz,jlmdz), wrf_psfctend(ilmdz,jlmdz), &
	1271	sfcGeopotValues(lmdzp)
	1272	END IF
	1273
	1274	PRINT *,' Lluis psfctendValues ______________ '
	1275	DO ix=1,4*INT(ddimxy/4),4
	1276	PRINT *,(ix+iy,':',psfctendValues(ix+iy), iy=0,3)
	1277	END DO
	1278	PRINT ,(iy,':',psfctendValues(iy), iy=4INT(ddimxy/4)+1,ddimxy)
	1279
	1280	! All the tendencies have to be given as decoupled from the dry air mass column! (see dyn_em/module_em.F:rk_update_scalar)
	1281	!!
	1282	DO iz=1,dimz
	1283	wrf_Ttend(:,iz,:) = wrf_Ttend(:,iz,:)*wrf_mut
	1284	wrf_Utend(:,iz,:) = wrf_Utend(:,iz,:)*wrf_muu
	1285	wrf_Vtend(:,iz,:) = wrf_Vtend(:,iz,:)*wrf_muv
	1286	DO iq=1,wnum3dm
	1287	wrf_MOISTtend(:,iz,:,iq) = wrf_MOISTtend(:,iz,:,iq)*wrf_mut
	1288	END DO
	1289	END DO
	1290
	1291	CALL get_lmdz(ddimx,ddimy,dimz,wbdyw,ddimxy,4,11,rlon,rlat,zmasq, &
	1292	pctsrf,ftsol,ftsoil,zxqsurf,qsurf,falb1,falb2,fevap,zxsnow,radsol,solsw, &
	1293	sollw,fder,rain_fall,snow_fall,agesno,zmea,zstd,zgam,zthe,zpic,zval,rugoro, &
	1294	frugs,zmax0,f0,wake_s,wake_cstar,wake_pe,wake_fip,sst,albedo, &
	1295	t_ancien,q_ancien,u_ancien,v_ancien,clwcon,rnebcon,ratqs,ema_work1,ema_work2, &
	1296	wake_deltat,wake_deltaq,fm_therm,entr_therm,detr_therm,phy_glo,wrf_grid)
	1297	! CALL fonte_neige_final(runoff_lic)
	1298
	1299	! Some re-start variables
	1300	PRINT *,' Lluis lrunofflic: ',LBOUND(wrf_grid%lrunofflic), ':', &
	1301	UBOUND(wrf_grid%lrunofflic),' ruonff_lic: ',UBOUND(runoff_lic),' ddimx: ', &
	1302	ddimx,' ddimy: ',ddimy
	1303
	1304	! PRINT ,' iy ix ddimx(iy-1)+ix runoff_lic wrf_grid%lrunofflic_______'
	1305	! DO iy=1,ddimy
	1306	! DO ix=1,ddimx
	1307	! PRINT ,iy,ix,ddimx(iy-1)+ix
	1308	! wrf_grid%lrunofflic(ix,iy) = runoff_lic(ddimx*(iy-1)+ix)
	1309	! PRINT ,iy,ix,ddimx(iy-1)+ix,runoff_lic(ddimx*(iy-1)+ix),wrf_grid%lrunofflic(ix,iy)
	1310	! END DO
	1311	! END DO
	1312
	1313	! STOP
	1314	!! CALL get_lmdz(ddimx,ddimy,dimz,wbdyw,ddimxy,4,11,zmasq,pctsrf,ftsol, &
	1315	!! ftsoil,fqsurf,qsol_d,falb1,falb2,fevap,zxsnow,radsol,solsw,sollw,fder, &
	1316	!! rain_fall,snow_fall,agesno,zmea,zstd,zgam,zthe,zpic,zval,rugoro, &
	1317	!! frugs,run_off_lic,zmax0,f0,wake_s,wake_cstar,wake_pe,wake_fip,sst, &
	1318	!! albedo,t_ancien,q_ancien,u_ancien,v_ancien,clwcon,rnebcon,ratqs, &
	1319	!! ema_work1,ema_work2,wake_deltat,wake_deltaq,fm_therm,entr_therm, &
	1320	!! detr_therm,phy_glo,wrf_grid)
	1321
	1322	PRINT *,' Lluis after get '//TRIM(Spt)//' ltksoil: ',wrf_grid%ltksoil(ip,:,jp), &
	1323	' ftsol: ',ftsol(lp,:)
	1324
	1325	PRINT *,' Lluis after get '//TRIM(Spt)//' fqsurf: ',zxqsurf(lp),' qsurf_rst: ', &
	1326	qsurf_rst(lp,:),' smois :',wrf_grid%smois(ip,:,jp),' lqksoil: ', &
	1327	wrf_grid%lqksoil(ip,:,jp)
	1328
	1329	! Recomputig again the total roughness length
	1330	LMDZvarmethod = 'prod'
	1331
	1332	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
	1333	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
	1334	wrf_grid%lrugksoil(:,1,:),wrf_grid%lrugksoil(:,2,:),wrf_grid%lrugksoil(:,3,:), &
	1335	wrf_grid%lrugksoil(:,4,:),wrf_grid%lrug)
	1336
	1337	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
	1338	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
	1339	wrf_grid%lalbksoil(:,1,:),wrf_grid%lalbksoil(:,2,:),wrf_grid%lalbksoil(:,3,:), &
	1340	wrf_grid%lalbksoil(:,4,:),wrf_grid%albedo)
	1341
	1342	! L. Fita, LMD January 2014.
	1343	! It should work nicely in this way... but memory got crazy...
	1344	CALL get_lmdz_out(ddimx, ddimy, dimz, wbdyw, ddimxy, nbsrf, nsoilmx, &
	1345	& lmdzmixingratios, 1, 2, wrf_l_pbl, wrf_l_con, wrf_l_thermals, wrf_l_wake, &
	1346	& 0, wrf_grid%itimestep, wtime_step, FullPressValues, HalfPressValues, &
	1347	& GeopotValues, sfcGeopotValues, lmdzoutP, MixingRatioValues, &
	1348	& UtendValues, VtendValues, TtendValues, MixingRatiotendValues, &
	1349	& psfctendValues, wrf_grid)
	1350
	1351	ip=icheck_p
	1352	jp=jcheck_p
	1353
	1354	CALL put_lmdz_in_WRFout(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1355	& wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, &
	1356	& lmdzmixingratios, 1, 2, wrf_l_pbl, wrf_l_con, wrf_l_thermals, wrf_l_wake, &
	1357	& wrf_grid%itimestep, wtime_step, FullPressValues, HalfPressValues, &
	1358	& GeopotValues, sfcGeopotValues, lmdzoutP, MixingRatioValues, &
	1359	& UtendValues, VtendValues, TtendValues, MixingRatiotendValues, &
	1360	& psfctendValues, ORCHIDEElevels, ftsoil, WRF_nsoillayers, &
	1361	& wrf_grid%soil_layers/100., wrf_grid%zs, wrf_grid%dzs, wrf_grid%nest_pos, &
	1362	& wrf_grid%q2, wrf_grid%t2, wrf_grid%th2, wrf_grid%psfc, wrf_grid%u10, &
	1363	& wrf_grid%v10, wrf_grid%resm, wrf_grid%zetatop, wrf_grid%tslb,wrf_grid%smois,&
	1364	& wrf_grid%sh2o, wrf_grid%smcrel, wrf_grid%sfcrunoff, wrf_grid%udrunoff, &
	1365	& wrf_grid%grdflx, wrf_grid%acgrdflx, wrf_grid%snow, wrf_grid%snowh, &
	1366	& wrf_grid%canwat, wrf_grid%sst, wrf_grid%sstsk, wrf_grid%lai, &
	1367	& wrf_grid%h_diabatic, wrf_grid%sina, wrf_grid%tsk, wrf_grid%tiso, &
	1368	& wrf_grid%max_msftx, wrf_grid%max_msfty, wrf_grid%rainc, wrf_grid%raincv, &
	1369	& wrf_grid%rainsh, wrf_grid%rainshv, wrf_grid%rainnc, wrf_grid%rainncv, &
	1370	& wrf_grid%snownc, wrf_grid%snowncv, wrf_grid%graupelnc, wrf_grid%graupelncv, &
	1371	& wrf_grid%hailnc, wrf_grid%hailncv, wrf_grid%stepave_count, wrf_grid%cldfra, &
	1372	& wrf_grid%swdown, wrf_grid%glw, wrf_grid%swnorm, wrf_grid%olr, &
	1373	& wrf_grid%emiss, wrf_grid%tmn, wrf_grid%ust, wrf_grid%pblh, wrf_grid%hfx, &
	1374	& wrf_grid%qfx, wrf_grid%lh, wrf_grid%achfx, wrf_grid%aclhf, &
	1375	& wrf_grid%snowc, wrf_grid%sr, wrf_grid%potevp, wrf_grid%snopcx, &
	1376	& wrf_grid%soiltb)
	1377	! & wrf_grid%save_topo_from_real, &
	1378	! & wrf_grid%avg_fuel_frac, wrf_grid%uah, wrf_grid%vah, wrf_grid%seed1, &
	1379	! & wrf_grid%seed2)
	1380
	1381	ip=icheck_p
	1382	jp=jcheck_p
	1383
	1384
	1385	! So it is done in a more larger to type and less elegant way...
	1386	!
	1387	! CALL get_lmdz_out2D_i(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1388	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
	1389	! & wrf_l_wake, wrf_grid%itimestep, wtime_step, FullPressValues2, &
	1390	! & HalfPressValues2, GeopotValues2, &
	1391	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1392	! & UtendValues2, VtendValues2, TtendValues2, &
	1393	! & MixingRatiotendValues2, psfctendValues2, wrf_grid%lages_ter, &
	1394	! & wrf_grid%lages_lic, wrf_grid%lages_oce, wrf_grid%lages_sic, wrf_grid%laire, &
	1395	! & wrf_grid%laireter, wrf_grid%lalb1, wrf_grid%lalb2, wrf_grid%lalbe_ter, &
	1396	! & wrf_grid%lalbe_lic, wrf_grid%lalbe_oce, wrf_grid%lalbe_sic, wrf_grid%lale, &
	1397	! & wrf_grid%lale_bl, wrf_grid%lale_wk, wrf_grid%lalp, wrf_grid%lalp_bl, &
	1398	! & wrf_grid%lalp_wk, wrf_grid%lbils, wrf_grid%lbils_diss, wrf_grid%lbils_ec, &
	1399	! & wrf_grid%lbils_enthalp, wrf_grid%lbils_kinetic, wrf_grid%lbils_latent, &
	1400	! & wrf_grid%lbils_tke, wrf_grid%lcape, wrf_grid%lcape_max, wrf_grid%lcdrh, &
	1401	! & wrf_grid%lcdrm, wrf_grid%lcldh, wrf_grid%lcldl, wrf_grid%lcldm, &
	1402	! & wrf_grid%lcldq, wrf_grid%lcldt, wrf_grid%lcontfracatm, wrf_grid%lcontfracor,&
	1403	! & wrf_grid%ldthmin, wrf_grid%ldtsvdft, wrf_grid%ldtsvdfo, wrf_grid%ldtsvdfg, &
	1404	! & wrf_grid%ldtsvdfi)
	1405
	1406	! CALL get_lmdz_out2D_ii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1407	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
	1408	! & wrf_l_wake, wrf_grid%itimestep, wtime_step, FullPressValues2, &
	1409	! & HalfPressValues2, GeopotValues2, &
	1410	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1411	! & UtendValues2, VtendValues2, TtendValues2, &
	1412	! & MixingRatiotendValues2, psfctendValues2, wrf_grid%levap, wrf_grid%levap_ter,&
	1413	! & wrf_grid%levap_lic, &
	1414	! & wrf_grid%levap_oce, wrf_grid%levap_sic, wrf_grid%levappot_ter, &
	1415	! & wrf_grid%levappot_lic, wrf_grid%levappot_oce, wrf_grid%levappot_sic, &
	1416	! & wrf_grid%lfbase, wrf_grid%lffonte, wrf_grid%lflat, wrf_grid%lflw_ter, &
	1417	! & wrf_grid%lflw_lic, wrf_grid%lflw_oce, wrf_grid%lflw_sic, &
	1418	! & wrf_grid%lfqcalving, wrf_grid%lfqfonte, wrf_grid%lfract_ter, &
	1419	! & wrf_grid%lfract_lic, wrf_grid%lfract_oce, wrf_grid%lfract_sic, &
	1420	! & wrf_grid%lfsnow, wrf_grid%lfsw_ter, wrf_grid%lfsw_lic, wrf_grid%lfsw_oce, &
	1421	! & wrf_grid%lfsw_sic, wrf_grid%lftime_th, wrf_grid%liwp, wrf_grid%llat_ter, &
	1422	! & wrf_grid%llat_lic, wrf_grid%llat_oce, wrf_grid%llat_sic, wrf_grid%llmaxth, &
	1423	! & wrf_grid%llwp, wrf_grid%lmsnow)
	1424
	1425	! CALL get_lmdz_out2D_iii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1426	! & wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, &
	1427	! & wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
	1428	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1429	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1430	! & UtendValues2, VtendValues2, TtendValues2, &
	1431	! & MixingRatiotendValues2, psfctendValues2,wrf_grid%lndayrain,wrf_grid%lnettop,&
	1432	! & wrf_grid%lpbase, wrf_grid%lphis, wrf_grid%lplcl, wrf_grid%lplfc, &
	1433	! & wrf_grid%lpluc, wrf_grid%lplul, wrf_grid%lplulst, wrf_grid%lplulth, &
	1434	! & wrf_grid%lpourc_ter, wrf_grid%lpourc_lic, wrf_grid%lpourc_oce, &
	1435	! & wrf_grid%lpourc_sic, wrf_grid%lprecip, wrf_grid%lprw, wrf_grid%lpsol, &
	1436	! & wrf_grid%lptop, wrf_grid%lqsat2m, wrf_grid%lqsol, wrf_grid%lqsurf, &
	1437	! & wrf_grid%lradsol, wrf_grid%lrh2m, wrf_grid%lrh2m_max, wrf_grid%lrh2m_min, &
	1438	! & wrf_grid%lrugs, wrf_grid%lrugs_ter, wrf_grid%lrugs_lic, wrf_grid%lrugs_oce, &
	1439	! & wrf_grid%lrugs_sic, wrf_grid%lsens, wrf_grid%lsicf, wrf_grid%ls_lcl, &
	1440	! & wrf_grid%lslp, wrf_grid%lsnow, wrf_grid%lsnowl, wrf_grid%lsoll, &
	1441	! & wrf_grid%lsoll0, wrf_grid%lsollwdown, wrf_grid%lsols, wrf_grid%lsols0, &
	1442	! & wrf_grid%ls_pblh, wrf_grid%ls_pblt, wrf_grid%ls_therm)
	1443
	1444	! CALL get_lmdz_out2D_iv(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1445	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
	1446	! & wrf_l_wake, wrf_grid%itimestep, &
	1447	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1448	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1449	! & UtendValues2, VtendValues2, TtendValues2, &
	1450	! & MixingRatiotendValues2, psfctendValues2, wrf_grid%lt2m, &
	1451	! & wrf_grid%lt2m_max, wrf_grid%lt2m_min, wrf_grid%lt2m_ter, wrf_grid%lt2m_lic, &
	1452	! & wrf_grid%lt2m_oce, wrf_grid%lt2m_sic, wrf_grid%ltaux, wrf_grid%ltauy, &
	1453	! & wrf_grid%ltaux_ter, wrf_grid%ltaux_lic, wrf_grid%ltaux_oce, &
	1454	! & wrf_grid%ltaux_sic, wrf_grid%ltauy_ter, wrf_grid%ltauy_lic, &
	1455	! & wrf_grid%ltauy_oce, wrf_grid%ltauy_sic, wrf_grid%lt_oce_sic, wrf_grid%ltopl,&
	1456	! & wrf_grid%ltopl0, wrf_grid%ltops, wrf_grid%ltops0, wrf_grid%ltpot, &
	1457	! & wrf_grid%ltpote, wrf_grid%ltsol, wrf_grid%ltsol_ter, wrf_grid%ltsol_lic, &
	1458	! & wrf_grid%ltsol_oce, wrf_grid%ltsol_sic)
	1459
	1460	! CALL get_lmdz_out2D_v(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1461	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
	1462	! & wrf_l_wake, wrf_grid%itimestep, &
	1463	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1464	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1465	! & UtendValues2, VtendValues2, TtendValues2, &
	1466	! & MixingRatiotendValues2, psfctendValues2,wrf_grid%lu10m, wrf_grid%lu10m_ter, &
	1467	! & wrf_grid%lu10m_lic, wrf_grid%lu10m_oce, wrf_grid%lu10m_sic, wrf_grid%lue, &
	1468	! & wrf_grid%luq, wrf_grid%lustar, wrf_grid%lustar_ter, wrf_grid%lustar_lic, &
	1469	! & wrf_grid%lustar_oce, wrf_grid%lustar_sic, wrf_grid%lv10m, &
	1470	! & wrf_grid%lv10m_ter, wrf_grid%lv10m_lic, wrf_grid%lv10m_oce, &
	1471	! & wrf_grid%lv10m_sic, wrf_grid%lve, wrf_grid%lvq, wrf_grid%lwake_h, &
	1472	! & wrf_grid%lwake_s, wrf_grid%lwbeff, wrf_grid%lwbilo_ter, wrf_grid%lwbilo_lic,&
	1473	! & wrf_grid%lwbilo_oce, wrf_grid%lwbilo_sic, wrf_grid%lwbils_ter, &
	1474	! & wrf_grid%lwbils_lic, wrf_grid%lwbils_oce, wrf_grid%lwbils_sic, &
	1475	! & wrf_grid%lweakinv, wrf_grid%lwind10m, wrf_grid%lwind10max, wrf_grid%lzmax_th)
	1476
	1477	! DO iz=wkms, wkme-2
	1478	! CALL get_lmdz_out3D_z_i(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, wbdyw, &
	1479	! & ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_con, &
	1480	! & wrf_l_thermals, wrf_l_wake, &
	1481	! & wrf_grid%itimestep, &
	1482	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1483	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1484	! & UtendValues2, VtendValues2, TtendValues2, &
	1485	! & MixingRatiotendValues2, psfctendValues2, &
	1486	! & wrf_grid%la_th(wims:wime,iz,wjms:wjme), &
	1487	! & wrf_grid%lbeta_prec(wims:wime,iz,wjms:wjme), &
	1488	! & wrf_grid%ldmc(wims:wime,iz,wjms:wjme), &
	1489	! & wrf_grid%ldnwd(wims:wime,iz,wjms:wjme), &
	1490	! & wrf_grid%ldnwd0(wims:wime,iz,wjms:wjme), &
	1491	! & wrf_grid%ldqajs(wims:wime,iz,wjms:wjme), &
	1492	! & wrf_grid%ldqcon(wims:wime,iz,wjms:wjme), &
	1493	! & wrf_grid%ldqdyn(wims:wime,iz,wjms:wjme), &
	1494	! & wrf_grid%ldqeva(wims:wime,iz,wjms:wjme), &
	1495	! & wrf_grid%ldqlsc(wims:wime,iz,wjms:wjme), &
	1496	! & wrf_grid%ldqlscst(wims:wime,iz,wjms:wjme), &
	1497	! & wrf_grid%ldqlscth(wims:wime,iz,wjms:wjme), &
	1498	! & wrf_grid%ldqphy(wims:wime,iz,wjms:wjme), &
	1499	! & wrf_grid%ldqthe(wims:wime,iz,wjms:wjme), &
	1500	! & wrf_grid%ldqvdf(wims:wime,iz,wjms:wjme), &
	1501	! & wrf_grid%ldqwak(wims:wime,iz,wjms:wjme))
	1502
	1503	! CALL get_lmdz_out3D_z_ii(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, &
	1504	! & wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, &
	1505	! & wrf_l_thermals, wrf_l_wake, &
	1506	! & wrf_grid%itimestep, &
	1507	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1508	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1509	! & UtendValues2, VtendValues2, TtendValues2, &
	1510	! & MixingRatiotendValues2, psfctendValues2, &
	1511	! & wrf_grid%ldtajs(wims:wime,iz,wjms:wjme), &
	1512	! & wrf_grid%ldtcon(wims:wime,iz,wjms:wjme), &
	1513	! & wrf_grid%ldtdis(wims:wime,iz,wjms:wjme), &
	1514	! & wrf_grid%ldtdyn(wims:wime,iz,wjms:wjme), &
	1515	! & wrf_grid%ldtec(wims:wime,iz,wjms:wjme), &
	1516	! & wrf_grid%ldteva(wims:wime,iz,wjms:wjme), &
	1517	! & wrf_grid%ld_th(wims:wime,iz,wjms:wjme), &
	1518	! & wrf_grid%lcldemi(wims:wime,iz,wjms:wjme), &
	1519	! & wrf_grid%lcldtau(wims:wime,iz,wjms:wjme), &
	1520	! & wrf_grid%lclwcon(wims:wime,iz,wjms:wjme), &
	1521	! & wrf_grid%ldtlif(wims:wime,iz,wjms:wjme), &
	1522	! & wrf_grid%ldtlsc(wims:wime,iz,wjms:wjme), &
	1523	! & wrf_grid%ldtlschr(wims:wime,iz,wjms:wjme), &
	1524	! & wrf_grid%ldtlscst(wims:wime,iz,wjms:wjme), &
	1525	! & wrf_grid%ldtlscth(wims:wime,iz,wjms:wjme), &
	1526	! & wrf_grid%ldtlw0(wims:wime,iz,wjms:wjme), &
	1527	! & wrf_grid%ldtlwr(wims:wime,iz,wjms:wjme), &
	1528	! & wrf_grid%ldtoro(wims:wime,iz,wjms:wjme), &
	1529	! & wrf_grid%ldtphy(wims:wime,iz,wjms:wjme), &
	1530	! & wrf_grid%ldtsw0(wims:wime,iz,wjms:wjme), &
	1531	! & wrf_grid%ldtswr(wims:wime,iz,wjms:wjme), &
	1532	! & wrf_grid%ldtthe(wims:wime,iz,wjms:wjme), &
	1533	! & wrf_grid%ldtvdf(wims:wime,iz,wjms:wjme), &
	1534	! & wrf_grid%ldtwak(wims:wime,iz,wjms:wjme))
	1535
	1536	! CALL get_lmdz_out3D_z_iii(wims, wime, wjms, wjme, iz, ddimx, ddimy, &
	1537	! & dimz, wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, &
	1538	! & wrf_l_con, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
	1539	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1540	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1541	! & UtendValues2, VtendValues2, TtendValues2, &
	1542	! & MixingRatiotendValues2, psfctendValues2, &
	1543	! & wrf_grid%lducon(wims:wime,iz,wjms:wjme), &
	1544	! & wrf_grid%ldudyn(wims:wime,iz,wjms:wjme), &
	1545	! & wrf_grid%ldulif(wims:wime,iz,wjms:wjme), &
	1546	! & wrf_grid%lduoro(wims:wime,iz,wjms:wjme), &
	1547	! & wrf_grid%lduvdf(wims:wime,iz,wjms:wjme), &
	1548	! & wrf_grid%ldvcon(wims:wime,iz,wjms:wjme), &
	1549	! & wrf_grid%ldvdyn(wims:wime,iz,wjms:wjme), &
	1550	! & wrf_grid%ldvlif(wims:wime,iz,wjms:wjme), &
	1551	! & wrf_grid%ldvoro(wims:wime,iz,wjms:wjme), &
	1552	! & wrf_grid%ldvvdf(wims:wime,iz,wjms:wjme), &
	1553	! & wrf_grid%le_th(wims:wime,iz,wjms:wjme), &
	1554	! & wrf_grid%levu(wims:wime,iz,wjms:wjme), &
	1555	! & wrf_grid%lfl(wims:wime,iz,wjms:wjme), &
	1556	! & wrf_grid%lgeop(wims:wime,iz,wjms:wjme), &
	1557	! & wrf_grid%lh2o(wims:wime,iz,wjms:wjme), &
	1558	! & wrf_grid%liwcon(wims:wime,iz,wjms:wjme), &
	1559	! & wrf_grid%lkz(wims:wime,iz,wjms:wjme), &
	1560	! & wrf_grid%lkz_max(wims:wime,iz,wjms:wjme), &
	1561	! & wrf_grid%llambda_th(wims:wime,iz,wjms:wjme), &
	1562	! & wrf_grid%llwcon(wims:wime,iz,wjms:wjme), &
	1563	! & wrf_grid%lma(wims:wime,iz,wjms:wjme), &
	1564	! & wrf_grid%lmc(wims:wime,iz,wjms:wjme), &
	1565	! & wrf_grid%lmcd(wims:wime,iz,wjms:wjme), &
	1566	! & wrf_grid%lmass(wims:wime,iz,wjms:wjme))
	1567
	1568	! CALL get_lmdz_out3D_z_iv(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, &
	1569	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
	1570	! & wrf_l_wake, wrf_grid%itimestep, &
	1571	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1572	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1573	! & UtendValues2, VtendValues2, TtendValues2, &
	1574	! & MixingRatiotendValues2, psfctendValues2, &
	1575	! & wrf_grid%loliq(wims:wime,iz,wjms:wjme), &
	1576	! & wrf_grid%lovap(wims:wime,iz,wjms:wjme), &
	1577	! & wrf_grid%lovapinit(wims:wime,iz,wjms:wjme), &
	1578	! & wrf_grid%lpaprs(wims:wime,iz,wjms:wjme), &
	1579	! & wrf_grid%lpr_con_i(wims:wime,iz,wjms:wjme), &
	1580	! & wrf_grid%lpr_con_l(wims:wime,iz,wjms:wjme), &
	1581	! & wrf_grid%lpres(wims:wime,iz,wjms:wjme), &
	1582	! & wrf_grid%lpr_lsc_i(wims:wime,iz,wjms:wjme), &
	1583	! & wrf_grid%lpr_lsc_l(wims:wime,iz,wjms:wjme), &
	1584	! & wrf_grid%lptconv(wims:wime,iz,wjms:wjme), &
	1585	! & wrf_grid%lq_th(wims:wime,iz,wjms:wjme), &
	1586	! & wrf_grid%lratqs(wims:wime,iz,wjms:wjme), &
	1587	! & wrf_grid%lre(wims:wime,iz,wjms:wjme), &
	1588	! & wrf_grid%lref_ice(wims:wime,iz,wjms:wjme), &
	1589	! & wrf_grid%lref_liq(wims:wime,iz,wjms:wjme), &
	1590	! & wrf_grid%lrhum(wims:wime,iz,wjms:wjme), &
	1591	! & wrf_grid%lrld(wims:wime,iz,wjms:wjme), &
	1592	! & wrf_grid%lrldcs(wims:wime,iz,wjms:wjme), &
	1593	! & wrf_grid%lrlu(wims:wime,iz,wjms:wjme), &
	1594	! & wrf_grid%lrlucs(wims:wime,iz,wjms:wjme), &
	1595	! & wrf_grid%lrneb(wims:wime,iz,wjms:wjme), &
	1596	! & wrf_grid%lrnebcon(wims:wime,iz,wjms:wjme), &
	1597	! & wrf_grid%lrnebls(wims:wime,iz,wjms:wjme), &
	1598	! & wrf_grid%lrsd(wims:wime,iz,wjms:wjme), &
	1599	! & wrf_grid%lrsdcs(wims:wime,iz,wjms:wjme), &
	1600	! & wrf_grid%lrsu(wims:wime,iz,wjms:wjme), &
	1601	! & wrf_grid%lrsucs(wims:wime,iz,wjms:wjme))
	1602
	1603	! CALL get_lmdz_out3D_z_v(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, &
	1604	! & wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, 1, &
	1605	! & wrf_l_con, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
	1606	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1607	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1608	! & UtendValues2, VtendValues2, TtendValues2, &
	1609	! & MixingRatiotendValues2, psfctendValues2, &
	1610	! & wrf_grid%lsens_ter(wims:wime,iz,wjms:wjme), &
	1611	! & wrf_grid%lsens_lic(wims:wime,iz,wjms:wjme), &
	1612	! & wrf_grid%lsens_oce(wims:wime,iz,wjms:wjme), &
	1613	! & wrf_grid%lsens_sic(wims:wime,iz,wjms:wjme), &
	1614	! & wrf_grid%ltemp(wims:wime,iz,wjms:wjme), &
	1615	! & wrf_grid%ltheta(wims:wime,iz,wjms:wjme), &
	1616	! & wrf_grid%ltke(wims:wime,iz,wjms:wjme), &
	1617	! & wrf_grid%ltke_max(wims:wime,iz,wjms:wjme), &
	1618	! & wrf_grid%ltke_ter(wims:wime,iz,wjms:wjme), &
	1619	! & wrf_grid%ltke_lic(wims:wime,iz,wjms:wjme), &
	1620	! & wrf_grid%ltke_oce(wims:wime,iz,wjms:wjme), &
	1621	! & wrf_grid%ltke_sic(wims:wime,iz,wjms:wjme), &
	1622	! & wrf_grid%ltnhus(wims:wime,iz,wjms:wjme), &
	1623	! & wrf_grid%ltnhusc(wims:wime,iz,wjms:wjme), &
	1624	! & wrf_grid%ltnhusscpbl(wims:wime,iz,wjms:wjme), &
	1625	! & wrf_grid%ltnt(wims:wime,iz,wjms:wjme), &
	1626	! & wrf_grid%ltntc(wims:wime,iz,wjms:wjme), &
	1627	! & wrf_grid%ltntr(wims:wime,iz,wjms:wjme), &
	1628	! & wrf_grid%ltntscpbl(wims:wime,iz,wjms:wjme), &
	1629	! & wrf_grid%lupwd(wims:wime,iz,wjms:wjme), &
	1630	! & wrf_grid%lvitu(wims:wime,iz,wjms:wjme), &
	1631	! & wrf_grid%lvitv(wims:wime,iz,wjms:wjme), &
	1632	! & wrf_grid%lvitw(wims:wime,iz,wjms:wjme), &
	1633	! & wrf_grid%lvprecip(wims:wime,iz,wjms:wjme), &
	1634	! & wrf_grid%lwake_omg(wims:wime,iz,wjms:wjme), &
	1635	! & wrf_grid%lwdtraina(wims:wime,iz,wjms:wjme), &
	1636	! & wrf_grid%lwdtrainm(wims:wime,iz,wjms:wjme), &
	1637	! & wrf_grid%lzfull(wims:wime,iz,wjms:wjme), &
	1638	! & wrf_grid%lzhalf(wims:wime,iz,wjms:wjme))
	1639
	1640	! END DO
	1641
	1642	! CALL get_lmdz_out3D_ii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1643	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
	1644	! & wrf_l_wake, wrf_grid%itimestep, &
	1645	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1646	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1647	! & UtendValues2, VtendValues2, TtendValues2, &
	1648	! & MixingRatiotendValues2, psfctendValues2, &
	1649	! & wrf_grid%ldtajs, &
	1650	! & wrf_grid%ldtcon, wrf_grid%ldtdis, wrf_grid%ldtdyn, wrf_grid%ldtec, &
	1651	! & wrf_grid%ldteva, wrf_grid%ld_th, wrf_grid%lcldemi, wrf_grid%lcldtau, &
	1652	! & wrf_grid%lclwcon, wrf_grid%ldtlif, wrf_grid%ldtlsc, wrf_grid%ldtlschr, &
	1653	! & wrf_grid%ldtlscst, wrf_grid%ldtlscth, wrf_grid%ldtlw0, wrf_grid%ldtlwr, &
	1654	! & wrf_grid%ldtoro, wrf_grid%ldtphy, wrf_grid%ldtsw0, wrf_grid%ldtswr, &
	1655	! & wrf_grid%ldtthe, wrf_grid%ldtvdf, wrf_grid%ldtwak)
	1656
	1657	! CALL get_lmdz_out3D_iii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1658	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
	1659	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1660	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1661	! & UtendValues2, VtendValues2, TtendValues2, &
	1662	! & MixingRatiotendValues2, psfctendValues2, &
	1663	! & wrf_grid%lducon, &
	1664	! & wrf_grid%ldudyn, wrf_grid%ldulif, wrf_grid%lduoro, wrf_grid%lduvdf, &
	1665	! & wrf_grid%ldvcon, wrf_grid%ldvdyn, wrf_grid%ldvlif, wrf_grid%ldvoro, &
	1666	! & wrf_grid%ldvvdf, wrf_grid%le_th, wrf_grid%levu, wrf_grid%lfl,wrf_grid%lgeop,&
	1667	! & wrf_grid%lh2o, wrf_grid%liwcon, wrf_grid%lkz, wrf_grid%lkz_max, &
	1668	! & wrf_grid%llambda_th, wrf_grid%llwcon, wrf_grid%lma, wrf_grid%lmass)
	1669
	1670	! CALL get_lmdz_out3D_iv(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1671	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
	1672	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1673	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1674	! & UtendValues2, VtendValues2, TtendValues2, &
	1675	! & MixingRatiotendValues2, psfctendValues2, &
	1676	! & wrf_grid%loliq, wrf_grid%lovap, &
	1677	! & wrf_grid%lovapinit, wrf_grid%lpaprs, wrf_grid%lpr_con_i, wrf_grid%lpr_con_l,&
	1678	! & wrf_grid%lpres, wrf_grid%lpr_lsc_i, wrf_grid%lpr_lsc_l, wrf_grid%lptconv, &
	1679	! & wrf_grid%lq_th, wrf_grid%lratqs, wrf_grid%lre, wrf_grid%lref_ice, &
	1680	! & wrf_grid%lref_liq, wrf_grid%lrhum, wrf_grid%lrld, wrf_grid%lrldcs, &
	1681	! & wrf_grid%lrlu, wrf_grid%lrlucs, wrf_grid%lrneb, wrf_grid%lrnebcon, &
	1682	! & wrf_grid%lrnebls, wrf_grid%lrsd, wrf_grid%lrsdcs, wrf_grid%lrsu, &
	1683	! & wrf_grid%lrsucs)
	1684
	1685	! CALL get_lmdz_out3D_v(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
	1686	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
	1687	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
	1688	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
	1689	! & UtendValues2, VtendValues2, TtendValues2, &
	1690	! & MixingRatiotendValues2, psfctendValues2, &
	1691	! & wrf_grid%lsens_ter, wrf_grid%lsens_lic, wrf_grid%lsens_oce, &
	1692	! & wrf_grid%lsens_sic, wrf_grid%ltemp, wrf_grid%ltheta, wrf_grid%ltke, &
	1693	! & wrf_grid%ltke_max, wrf_grid%ltke_ter, wrf_grid%ltke_lic, wrf_grid%ltke_oce, &
	1694	! & wrf_grid%ltke_sic, wrf_grid%ltnhus, wrf_grid%ltnhusc, wrf_grid%ltnhusscpbl, &
	1695	! & wrf_grid%ltnt, wrf_grid%ltntc, wrf_grid%ltntr, wrf_grid%ltntscpbl, &
	1696	! & wrf_grid%lupwd, wrf_grid%lvitu, wrf_grid%lvitv, wrf_grid%lvitw, &
	1697	! & wrf_grid%lwake_deltaq, wrf_grid%lwake_deltat, wrf_grid%lwake_omg, &
	1698	! & wrf_grid%lwdtraina, wrf_grid%lwdtrainm, wrf_grid%lzfull, wrf_grid%lzhalf)
	1699
	1700	DEALLOCATE(MixingRatioValues)
	1701	DEALLOCATE(MixingRatiotendValues)
	1702
	1703	WRITE(message, *)' exiting: module_lmdz_phys.F'
	1704	CALL wrf_debug(75, message)
	1705
	1706	END SUBROUTINE call_lmdz_phys
	1707
	1708	END MODULE module_lmdz_phys

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

Download in other formats:

Original Format