Context Navigation

module_lmdz_phys.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 98.7 KB

Line
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
93	INTEGER, INTENT(IN) :: wjulday
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
237	! wgmt: WRF gmt hour
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	PRINT *,' '//TRIM(Spt)//': k p ph ___________'
378	DO iz=1,dimz
379	PRINT *,iz,wrf_perP(ip,iz,jp),wrf_perGeopot(ip,iz,jp)
380	END DO
381	PRINT *,'dimx: ',dimx,' dimy: ',dimy,' dimz: ',dimz,' ddimx: ',ddimx,' ddimy: ',ddimy
382
383	phykeys=0.
384	lmdz_dudyn=0.
385	lmdz_PVtheta=0.
386
387	! Checking of first/last time-step for LMDZ
388	!!
389	lmdz_debut = .FALSE.
390	lmdz_lafin = .FALSE.
391
392	! PRINT *,'Lluis: i j ij lon lat lonvals latvals_______'
393	DO iy=1,ddimy
394	DO ix=1,ddimx
395	IF (wrf_Lon(ix,iy) < 0.) THEN
396	LonValues(ddimx(iy-1)+ix) = (360. + wrf_Lon(ix,iy))Pi/180.
397	ELSE
398	LonValues(ddimx(iy-1)+ix) = wrf_Lon(ix,iy)Pi/180.
399	END IF
400	LatValues(ddimx(iy-1)+ix) = wrf_Lat(ix,iy)Pi/180.
401	PsfctendValues(ddimx*(iy-1)+ix) = wrf_psfctend(ix,iy)
402	! PRINT ,ix,iy,(ddimx)(iy-1)+ix,wrf_Lon(ix,iy),LonValues(ddimx*(iy-1)+ix), &
403	! wrf_Lat(ix,iy),LatValues(ddimx*(iy-1)+ix)
404	END DO
405	END DO
406
407	! Reshape can not be used? Different initial/ending indices !
408	! LonValues = RESHAPE(wrf_Lon, (/ dimx*dimy/), ORDER=(/1/))
409	! LatValues = RESHAPE(wrf_Lat, (/ dimx*dimy/), ORDER=(/1/))
410
411	firststep: IF ( wrf_xtime == 0 .OR. wrf_isrestart) THEN
412	lmdz_debut = .TRUE.
413
414	DO funit=10,100
415	INQUIRE(UNIT=funit, OPENED=is_used)
416	IF (.NOT. is_used) EXIT
417	END DO
418	lmdz_clesphy0 = 0.
419
420	CALL conf_gcm(funit, .TRUE., lmdz_clesphy0)
421	! L. Fita, wee need to provide the value
422	annee_ref = anneeref
423	day_ref = dayref
424
425	! Uploading LMDZ general configuration
426	! used instead call conf_gcm, since variables are not found in memory (or at least with a value)
427	defilen = 'gcm.def'
428	defvarn = 'iflag_phys'
429	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
430	iflag_phys = defint
431
432	! Is not working it has an space...
433	! defvarn = 'type_trac'
434	! CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
435	! type_trac = TRIM(defstr)
436
437	! used instead call conf_gcm, since variables are not found in memory (or at least with a value)
438	defilen = 'config.def'
439	defvarn = 'type_ocean'
440	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
441	type_ocean= TRIM(defstr)
442
443	! used instead call conf_gcm, since variables are not found in memory (or at least with a value)
444	! Looking for consistency for some namelist variables and *.def LMDZ values
445	PRINT ,' Checking consistency between namelist and LMDZ .def values'
446	defilen = 'physiq.def'
447	defvarn = 'iflag_pbl'
448	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
449	lpbl = defint
450	defvarn = 'iflag_con'
451	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
452	lcon = defint
453	defvarn = 'iflag_thermals'
454	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
455	lth = defint
456	defvarn = 'iflag_wake'
457	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
458	lwake = defint
459
460	IF (wrf_l_pbl /= lpbl) THEN
461	PRINT *,TRIM(errmsg)
462	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ pbl value: ', &
463	wrf_l_pbl,' differs from iflag_pbl= ',lpbl ,' in physiq.def !!!!'
464	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
465	CALL wrf_error_fatal(TRIM(wrf_err_message))
466	END IF
467	IF (wrf_l_con /= lcon) THEN
468	PRINT *,TRIM(errmsg)
469	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ convection value: ', &
470	wrf_l_con,' differs from iflag_con= ',lcon ,' in physiq.def !!!!'
471	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
472	CALL wrf_error_fatal(TRIM(wrf_err_message))
473	END IF
474	IF (wrf_l_thermals /= lth) THEN
475	PRINT *,TRIM(errmsg)
476	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ thermals value: ', &
477	wrf_l_thermals,' differs from iflag_thermals= ',lth ,' in physiq.def !!!!'
478	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
479	CALL wrf_error_fatal(TRIM(wrf_err_message))
480	END IF
481	IF (wrf_l_wake /= lwake) THEN
482	PRINT *,TRIM(errmsg)
483	PRINT *,' ' // TRIM(subname) // ' namelist LMDZ wakes value: ', &
484	wrf_l_wake,' differs from iflag_wake= ',lwake ,' in physiq.def !!!!'
485	WRITE(wrf_err_message,*)' wrong namelist.input value !!'
486	CALL wrf_error_fatal(TRIM(wrf_err_message))
487	END IF
488
489	! Initialization of the variables
490	!!
491
492	INQUIRE(FILE='traceur.def', EXIST=existsfile)
493
494	! L. Fita, LMD. August, this would be much more nice if it is in config_flags structure....
495	IF (existsfile) THEN
496	DO funit=10,100
497	INQUIRE(UNIT=funit, OPENED=is_used)
498	IF (.NOT. is_used) EXIT
499	END DO
500	OPEN(UNIT=funit, FILE='traceur.def', STATUS='OLD')
501	READ(funit,*)lmdzmixingratios
502	CLOSE(UNIT=funit)
503	ELSE
504	WRITE(message, *)errmsg
505	CALL wrf_debug(75, message)
506	WRITE(wrf_err_message,*)' required file traceur.def not found !!'
507	CALL wrf_error_fatal(TRIM(wrf_err_message))
508	END IF
509
510	CALL infotrac_init
511
512	! Initializing LMDZ matrices...
513
514	CALL Init_Phys_lmdz(1,ddimxy,dimz,1,(/ddimxy/))
515	CALL InitComgeomphy()
516
517	CALL phys_state_var_init(0)
518	!! CALL neige_initialize()
519	!! CALL limit_initialize()
520
521	!! CALL pbl_initialize()
522
523	PRINT *,'Lluis after initializing all LMDZ!'
524
525	rlon=LonValues
526	rlat=LatValues
527
528	xtime0: IF (wrf_xtime == 0.) THEN
529
530	nlanduse=UBOUND(wrf_grid%landusef, 2)
531	CALL WRFlanduse_LMDZKsoil(landcat,wims,wime,wjms,wjme,ddimx,ddimy,nlanduse, &
532	wrf_grid%landusef,wrf_grid%landmask,wrf_grid%lter,wrf_grid%llic, &
533	wrf_grid%loce,wrf_grid%lsic)
534
535	PRINT *,' Lluis ltksoil. '// Spt //' ter: ',wrf_grid%lter(ip,jp),' lic: ',&
536	wrf_grid%llic(ip,jp),' loce: ', wrf_grid%loce(ip,jp),' sic: ', &
537	wrf_grid%lsic(ip,jp)
538
539	! LANDMASK in LMDZ as fractions
540	wrf_grid%lmsk = wrf_grid%lter + wrf_grid%llic
541	! Pseudo WRF-initialization
542
543	! LANDMASK in LMDZ as 1/0
544	! wrf_grid%lmsk=0.
545	! DO ix=1,ddimx
546	! DO iy=1,ddimy
547	! IF (wrf_grid%lter(ix,iy) >= 0.5) THEN
548	! wrf_grid%lmsk(ix,iy)=1.
549	! wrf_grid%lter(ix,iy)=1.
550	! wrf_grid%llic(ix,iy)=0.
551	! wrf_grid%loce(ix,iy)=0.
552	! wrf_grid%lsic(ix,iy)=0.
553	! ELSE
554	! wrf_grid%lmsk(ix,iy)=0.
555	! wrf_grid%lter(ix,iy)=0.
556	! wrf_grid%llic(ix,iy)=0.
557	! wrf_grid%loce(ix,iy)=1.
558	! wrf_grid%lsic(ix,iy)=0.
559	! END IF
560	! END DO
561	! END DO
562	! End of Pseudo WRF-initialization
563	!!!!!
564
565	! Value is given at the first time-step, not any more, since it will be continuously
566	! used by LMDZ
567	PRINT *,' Lluis before '//TRIM(Spt)//' SMOIS: ',wrf_grid%smois(ip,:,jp), &
568	' dzs: ',wrf_grid%dzs,' lwsol: ',wrf_grid%lwsol(ip,jp)
569	! l. Fita, LMD. January 2014
570	! This should be in this way, but, because we can not use any scheme DZS=0., so
571	! direct met_em values (Registry modified) are used. soil layers are in [cm]!
572	! DO iz=1,SIZE(wrf_grid%soil_layers(ip,:,jp))
573	! wrf_grid%lwsol = wrf_grid%lwsol + wrf_grid%sm(:,iz,:)* &
574	! wrf_grid%soil_layers(:,iz,:)/100.
575	! END DO
576	PRINT *,' Lluis after: '//TRIM(Spt)//' SMOIS: ',wrf_grid%smois(ip,:,jp), &
577	' lwsol: ',wrf_grid%lwsol(ip,jp)
578	DO iz=1,4
579	wrf_grid%ltksoil(:,iz,:)=wrf_grid%st(:,1,:)
580	END DO
581	!
582	! Only using SST on that ocean grid points
583	DO ix=1,ddimx
584	DO iy=1,ddimy
585	IF (wrf_grid%loce(ix,iy) > 0.) THEN
586	! Full SST grid points must have a SST value from the wrfinput (WRF 1/0 landmask)
587	! but that LMDZ points with fractions of land/mask could not. Using the closest one
588	! if it is an isolated grid point, using TSK instead
589	IF (wrf_grid%sst(ix,iy) < 200.15) THEN
590	wrf_grid%ltksoil(ix,3,iy) = -9.
591	IF ( (ix > 1).AND.(ix < ddimx) .AND. (iy > 1).AND.(iy < ddimy) ) THEN
592	DO ixx=-1,1
593	DO iyy=-1,1
594	! PRINT *,ixx,iyy,wrf_grid%sst(ix+ixx,iy+iyy)
595	IF (wrf_grid%sst(ix+ixx,iy+iyy) > 200.15) THEN
596	wrf_grid%ltksoil(ix,3,iy) = wrf_grid%sst(ix+ixx,iy+iyy)
597	EXIT
598	END IF
599	END DO
600	IF ( wrf_grid%ltksoil(ix,3,iy) > 200.15) EXIT
601	END DO
602	IF ( wrf_grid%ltksoil(ix,3,iy) == -9.) wrf_grid%ltksoil(ix,3,iy) = &
603	wrf_grid%tsk(ix,iy)
604	ELSE
605	wrf_grid%ltksoil(ix,3,iy) = wrf_grid%tsk(ix,iy)
606	END IF
607	ELSE
608	wrf_grid%ltksoil(ix,3,iy) = wrf_grid%sst(ix,iy)
609	END IF
610	IF ( (wrf_grid%ltksoil(ix,3,iy) < 200.15) .OR. &
611	(wrf_grid%ltksoil(ix,3,iy) > 370.) ) THEN
612	PRINT *,TRIM(errmsg)
613	WRITE(wrf_err_message,*) ' ' // TRIM(subname) // &
614	' SST inconsistency at point: ', ix, ', ', iy,' with loce: ', &
615	wrf_grid%loce(ix,iy),' has a tsoil(oce)= ', &
616	wrf_grid%ltksoil(ix,3,iy),' sst: ',wrf_grid%sst(ix,iy),' tsk: ', &
617	wrf_grid%tsk(ix,iy)
618	CALL wrf_error_fatal(TRIM(wrf_err_message))
619	END IF
620	END IF
621	END DO
622	END DO
623
624	!! Soil temperature, humidity, water-content...
625	! This should go on the namelist... later!
626	ORCHIDEElevels = (/ 0.098, 0.1955, 0.5865, 1.3685, 2.9326, 6.0606, 12.3167, &
627	24.8289, 49.8534, 99.9022, 200. /)
628	PRINT *,' Lluis interpolating:',SIZE(wrf_grid%soil_layers(ip,:,jp)),' : ', &
629	' soil: ', &
630	wrf_grid%soil_layers(ip,SIZE(wrf_grid%soil_layers(ip,:,jp)):1:-1,jp), &
631	' st: ',wrf_grid%st(ip,:,jp), 11, ' \| ', ORCHIDEElevels,' soiltop: ', &
632	wrf_grid%soilctop(ip,:,jp),' soilbot: ',wrf_grid%soilcbot(ip,:,jp)
633
634	datasetname = soilcat
635	tablename = 'SOILPARM.TBL'
636
637	PRINT *,' Lluis reading dataset "' // TRIM(datasetname) // '" from "' // &
638	TRIM(tablename) // '"'
639
640	CALL table_characteristics(datasetname,tablename,begdataset,Nvars,Nrs,Ncs)
641
642	! L. Fita, LMD. Feburary 2014
643	! Again (as in the landuses case), the files content more categeories that the really ones used by the model...
644	Nsoillayers=SIZE(wrf_grid%soil_layers(ip,:,jp))
645	Nwrfsoiltypes=SIZE(wrf_grid%soilctop(ip,:,jp))
646	IF (Nrs > Nwrfsoiltypes) THEN
647	PRINT *,TRIM(wrnmsg)
648	PRINT *," '" // TRIM(datasetname)//"' from '"//TRIM(tablename)//"' has ", &
649	Nrs,' soil layers, meanwhile WRF data has ',Nwrfsoiltypes
650	PRINT *,' using only the first ', Nwrfsoiltypes, ' of the dataset !!!!'
651	Nrs=Nsoillayers
652	ELSE IF (Nrs < Nwrfsoiltypes) THEN
653	PRINT *,TRIM(errmsg)
654	PRINT *," '" // TRIM(datasetname)//"' from '"//TRIM(tablename)//"' has ", &
655	Nrs,' soil layers, meanwhile WRF data has ',Nwrfsoiltypes
656	WRITE(wrf_err_message,*)' wrong the dataset !!!!'
657	CALL wrf_error_fatal(TRIM(wrf_err_message))
658	END IF
659
660	IF (ALLOCATED(soilvals)) DEALLOCATE(soilvals)
661	ALLOCATE(soilvals(Nvars,Nrs,Ncs))
662
663	IF (ALLOCATED(soilvalns)) DEALLOCATE(soilvalns)
664	ALLOCATE(soilvalns(Nrs))
665
666	CALL read_table(datasetname, tablename, begdataset, Nvars, Nrs, Ncs, &
667	soilvals, soilvalns)
668
669	DO ix=1,ddimx
670	DO iy=1,ddimy
671	!L. Fita, LMD. January 2014 Soil layers are in cm!
672
673	IF ( (ix == ip) .AND. (iy == jp)) PRINT *,' Lluis '//TRIM(Spt)// &
674	' soiltop: ', wrf_grid%soilctop(ix,:,iy),' soilbot: ', &
675	wrf_grid%soilcbot(ix,:,iy),' layers: ', &
676	wrf_grid%soil_layers(ix,Nsoillayers:1:-1,iy)/100.,' sm: ', &
677	wrf_grid%sm(ix,:,iy)
678
679	! NOTE: wrf_dx, wrf_dy should be matrix !!!
680	CALL compute_TOTsoil_mass_water_values(datasetname, Nvars, Nrs, Ncs, &
681	soilvals, soilvalns, ix, iy, Nrs, wrf_grid%soilctop(ix,:,iy), &
682	wrf_grid%soilcbot(ix,:,iy), Nsoillayers, wrf_grid%sm(ix,:,iy), &
683	wrf_grid%soil_layers(ix,Nsoillayers:1:-1,iy)/100., wrf_dx, wrf_dy, &
684	soilTmass, wrf_grid%lwsol(ix,iy))
685
686	IF ( (ix == ip) .AND. (iy == jp)) PRINT *,'TOTAL soil mass at '// &
687	TRIM(Spt)//': ', soilTmass,' [kg] water content: ', &
688	wrf_grid%lwsol(ix,iy), ' [kg] equiv. water depth: ', &
689	wrf_grid%lwsol(ix,iy)/(wrf_dxwrf_dy1000.)
690
691	wrf_grid%lwsol(ix,iy) = wrf_grid%lwsol(ix,iy)/(wrf_dxwrf_dy1000.)
692
693	CALL interpolate_layers(SIZE(wrf_grid%soil_layers(ix,:,iy)), &
694	wrf_grid%soil_layers(ix,SIZE(wrf_grid%soil_layers(ix,:,iy)):1:-1,iy), &
695	wrf_grid%st(ix,:,iy),11, ORCHIDEElevels,wrf_grid%lotter(ix,:,iy))
696	IF ( (ix == ip) .AND. (iy == jp)) PRINT *,' Lluis orchidee temperature'//&
697	' at '// TRIM(Spt)//': ',wrf_grid%lotter(ix,:,iy)
698	! IF ( (ix == ip) .AND. (iy == jp)) STOP
699	END DO
700	END DO
701
702	wrf_grid%lotlic=wrf_grid%lotter
703	wrf_grid%lotoce=wrf_grid%lotter
704	wrf_grid%lotsic=wrf_grid%lotter
705
706	DEALLOCATE(soilvals)
707	DEALLOCATE(soilvalns)
708
709	! LMDZ does not distinguish among species, use it as a whole
710	wrf_grid%lclwcon=SUM(wrf_MOIST,4)
711	DO ix=wims,wime
712	DO iz=wkms,wkme
713	DO iy=wjms,wjme
714	IF (wrf_grid%lclwcon(ix,iz,iy) /= 0.) THEN
715	wrf_grid%lratqs=(wrf_MOIST(ix,iz,iy,wrf_qrid)+ &
716	wrf_MOIST(ix,iz,iy,wrf_qsid)+wrf_MOIST(ix,iz,iy,wrf_qiid)+ &
717	wrf_MOIST(ix,iz,iy,wrf_qhid)+wrf_MOIST(ix,iz,iy,wrf_qgid))/ &
718	wrf_grid%lclwcon(ix,iz,iy)
719	END IF
720	END DO
721	END DO
722	END DO
723
724	! wrf_grid%znt = wrf_grid%z0
725	! wrf_grid%albedo = wrf_grid%albbck
726	wrf_grid%lbils=0.
727	wrf_grid%lfder=0.
728	wrf_grid%lrnebcon=0.
729	DO ix=1,ddimx
730	DO iy=1,ddimy
731	ksoils = (/ wrf_grid%lter(ix,iy), wrf_grid%llic(ix,iy), &
732	wrf_grid%loce(ix,iy), wrf_grid%lsic(ix,iy) /)
733	ksoilmax = MAXVAL(ksoils)
734	DO iz = 1,4
735	IF (ksoilmax == ksoils(iz)) wrf_grid%lnat(ix,iy) = iz
736	END DO
737	END DO
738	END DO
739	wrf_grid%lrads = wrf_grid%swdown
740
741	wrf_grid%lcontrol(1) = wtime_step
742	defilen = 'config.def'
743	defvarn = 'nbapp_rad'
744	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
745	wrf_grid%lcontrol(6) = defint
746	wrf_grid%lcontrol(2) = NINT (86400./wtime_step/ FLOAT(defint) )
747	!L. Fita, LMD. February 2014: current value of atmospheric CO2 should be
748	!! co2_ppm0 + delta_co2_ppm, from 'carbon_cycle_mod.F90'
749	defvarn = 'co2_ppm'
750	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
751	wrf_grid%lcontrol(3) = defint
752	!!c co2_ppm0 : initial value of atmospheric CO2
753	wrf_grid%lcontrol(16) = defint
754	defvarn = 'solaire'
755	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
756	wrf_grid%lcontrol(4) = defreal
757	wrf_grid%lcontrol(5) = wrf_l_con
758
759	! Default value cycle_diurne = .T.
760	wrf_grid%lcontrol(7) = 1.
761	defvarn = 'soil_model'
762	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
763	IF(defbool) wrf_grid%lcontrol(8) = 1.
764	defilen = 'physiq.def'
765	defvarn = 'new_oliq'
766	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
767	IF(defbool) wrf_grid%lcontrol(9) = 1.
768	defvarn = 'ok_orodr'
769	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
770	IF(defbool) wrf_grid%lcontrol(10) = 1.
771	defvarn = 'ok_orolf'
772	CALL def_get_scalar_value(defilen, defvarn, defstr, defint, defreal, defbool)
773	IF(defbool) wrf_grid%lcontrol(11) = 1.
774
775	wrf_grid%lcontrol(13) = dayref + nday
776	wrf_grid%lcontrol(14) = annee_ref
777	wrf_grid%lcontrol(15) = wrf_grid%itimestep
778
779	! Let's split rugosity by Ksoils
780	!! Static roughness length
781	wrf_grid%lrug=wrf_grid%znt
782
783	!! Dynamic roughness length (It should work if we use WRF schemes...)
784	! wrf_grid%lrug=wrf_grid%znt
785	! First split by Ksoil types.
786	! NOTE: in somehow we should read again the table 'VEGPARM.TBL' and recompute the
787	! roughness length in case that appear a new soil type. (but what left below a
788	! melted glacier?)
789	!
790	LMDZvarmethod='prod'
791	wrf_grid%lrugksoil(:,1,:) = wrf_grid%znt
792
793	CALL lmdz_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,wrf_grid%snow, &
794	LMDZvarmethod,wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
795	wrf_grid%lsnowksoil(:,1,:),wrf_grid%lsnowksoil(:,2,:), &
796	wrf_grid%lsnowksoil(:,3,:),wrf_grid%lsnowksoil(:,4,:))
797	! All the iced soil types have the same roughtness length 0.001
798	wrf_grid%lrugksoil(:,2,:)=0.001
799	wrf_grid%lrugksoil(:,3,:)=0.0001
800	wrf_grid%lrugksoil(:,4,:)=0.001
801
802	! Recomputig again the total roughness length
803	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
804	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
805	wrf_grid%lrugksoil(:,1,:),wrf_grid%lrugksoil(:,2,:), &
806	wrf_grid%lrugksoil(:,3,:),wrf_grid%lrugksoil(:,4,:),wrf_grid%lrug)
807
808	! Air surface moisture for each soil type. First order?
809	DO iz=1,4
810	wrf_grid%lqksoil(:,iz,:)=wrf_grid%q2
811	END DO
812
813	! Let's split albedo by Ksoils
814	wrf_grid%lalbksoil=0.
815
816	DO iz=1,4
817	!! Static roughness length (changing on time in wrflwbdy_d01)
818	wrf_grid%lalbksoil(:,iz,:)=wrf_grid%albedo
819	!! Dynamic roughness length (It should work if we use WRF schemes...)
820	! wrf_grid%lalbksoil(:,iz,:)=wrf_grid%albedo
821	END DO
822	! First split by Ksoil types.
823	! NOTE: in somehow we should read again the table 'VEGPARM.TBL' and recompute the
824	! albedo in case that appear a new soil type. (but what left below a melted
825	! glacier?)
826	!
827	! All the iced soil types have the same albedo
828	! For snow/ice there is zmin=0.55, zmax=0.70, tacking the mean value: 0.625
829	wrf_grid%lalbksoil(:,2,:)=0.625
830	wrf_grid%lalbksoil(:,3,:)=0.08
831	wrf_grid%lalbksoil(:,4,:)=0.625
832
833	! Recomputig again the total roughness length
834	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
835	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
836	wrf_grid%lalbksoil(:,1,:),wrf_grid%lalbksoil(:,2,:), &
837	wrf_grid%lalbksoil(:,3,:),wrf_grid%lalbksoil(:,4,:),wrf_grid%albedo)
838	! Using the same values for the LW albedo
839	wrf_grid%llwalbksoil = wrf_grid%lalbksoil
840
841	! Orography values just beacuse there is no other way right now (they should come
842	! from geogrid.exe)
843	wrf_grid%lzmea = wrf_grid%ht
844	wrf_grid%lzstd = wrf_grid%topostdv
845	wrf_grid%lzsig = wrf_grid%slope
846	wrf_grid%lzgam = SQRT((wrf_grid%toposlpx2. + wrf_grid%toposlpy2.))/ &
847	(1. + wrf_grid%slope)
848	wrf_grid%lzthe = ACOS(wrf_grid%toposlpx/SQRT(wrf_grid%toposlpx**2. + &
849	wrf_grid%toposlpy**2.))
850	wrf_grid%lzpic = wrf_grid%ht*1.10
851	wrf_grid%lzval = wrf_grid%ht*0.90
852	wrf_grid%lzrugsrel = wrf_grid%lrugKsoil(:,1,:)
853
854	PRINT *,'Lluis Topography values at',ip, jp, ' zmea: ', &
855	wrf_grid%lzmea(ip,jp),' zstd: ',wrf_grid%lzstd(ip,jp),' zsig: ', &
856	wrf_grid%lzsig(ip,jp),' zgam: ',wrf_grid%lzgam(ip,jp),' zthe: ', &
857	wrf_grid%lzthe(ip,jp),' zpic: ',wrf_grid%lzpic(ip,jp),' zval: ', &
858	wrf_grid%lzval(ip,jp),' rugoro: ',wrf_grid%lzrugsrel(ip,jp)
859
860	! R = 8.3144621 ([J mol/K] http://en.wikipedia.org/wiki/Gas_constant)
861	! mu = 28.97 (Molecular Mass of Air [kg/kmol] from http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html)
862	! R / mu =? r_d
863
864	! lmdz_Area = dLon * dLat * cos(lat) * Earth_rad
865
866	END IF xtime0
867
868	! Grids do not start/end at the same index
869	!! PRINT *,' Lluis: Creating area/cu/cv...'
870	wrf_area = 0.
871	wrf_cu = 0.
872	wrf_cv = 0.
873	DO iy=1,ddimy
874	DO ix=1,ddimx
875	wrf_area(ix,iy) = wrf_dx*wrf_dy/wrf_mapft(ix,iy)
876	wrf_cu(ix,iy) = wrf_dx/wrf_mapfu(ix,iy)
877	wrf_cv(ix,iy) = wrf_dy/wrf_mapfv(ix,iy)
878	END DO
879	END DO
880
881	! cu = dLon * Earth_rad
882	! cv = dLat * Earth_rad * cos(Lat)
883
884	! lmdz_cu(1) = cu(1)
885	! lmdz_cv(1) = cv(1)
886	! DO j=2,dimy
887	! DO i=1,dimx
888	! lmdz_cu((j-2)dimx+1+i) = cu((j-1)(dimx+1)+i)
889	! lmdz_cv((j-2)dimx+1+i) = cv((j-1)(dimx+1)+i)
890	! ENDDO
891	! ENDDO
892	! lmdz_cu(dim_ngridmx) = cu(dimxy+1)
893	! lmdz_cv(dim_ngridmx) = cv(dimxy-dimx)
894
895	AreaValues = RESHAPE(wrf_area, (/ ddimx*ddimy/), ORDER=(/1/))
896	CUValues = RESHAPE(wrf_cu, (/ ddimx*ddimy/), ORDER=(/1/))
897	CVValues = RESHAPE(wrf_cv, (/ ddimx*ddimy/), ORDER=(/1/))
898
899	CALL lmdz_vars_init(ddimxy)
900
901	! PRINT *,' Lluis before loading 5 11 ltksoil: ',wrf_grid%ltksoil(5,:,11), &
902	! ' ftsol: ',ftsol(315,:)
903	! PRINT *,' Lluis before loading 5 11 qsurf: ',qsurf(315,:),' qsurf_rst: ', &
904	! qsurf_rst(315,:),' smois :',wrf_grid%smois(5,:,11),' lqksoil: ', &
905	! wrf_grid%lqksoil(5,:,11)
906
907	CALL load_lmdz(wrf_grid,ddimx,ddimy,dimz,wbdyw,ddimxy,4,11,rlon,rlat,zmasq, &
908	pctsrf,ftsol,ftsoil_rst,qsurf_rst,qsol_rst,falb1,falb2,evap_rst,snow_rst, &
909	radsol,solsw,sollw,fder_rst,rain_fall,snow_fall,agesno_rst,zmea,zstd,zgam, &
910	zthe,zpic,zval,rugoro,rugos_rst,restart_runoff,zmax0,f0,wake_s,wake_cstar, &
911	wake_pe,wake_fip,sst,albedo,t_ancien,q_ancien,u_ancien,v_ancien,clwcon, &
912	rnebcon,ratqs,ema_work1,ema_work2,wake_deltat,wake_deltaq,fm_therm, &
913	entr_therm,detr_therm,phy_glo)
914
915	CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst, &
916	rugos_rst, agesno_rst, ftsoil_rst)
917	CALL fonte_neige_init(restart_runoff)
918
919	PRINT *,' Lluis on module restart_runoff: ',UBOUND(restart_runoff)
920	PRINT *,' ffonte_global: ',ALLOCATED(ffonte_global)
921	PRINT *,' fqfonte_global: ',ALLOCATED(fqfonte_global)
922	PRINT *,' fqcalving_global: ',ALLOCATED(fqcalving_global)
923
924	PRINT *,' Lluis after loading '//TRIM(Spt)//': qsurf_rst: ', qsurf_rst(lp,:),&
925	' smois :',wrf_grid%smois(ip,:,jp),' lqksoil: ', wrf_grid%lqksoil(ip,:,jp)
926
927	PRINT *,' Lluis after loading '//TRIM(Spt)//' ltksoil: ', &
928	wrf_grid%ltksoil(ip,:,jp),' ftsol: ',ftsol(lp,:)
929
930	PRINT *,' Lluis after loading '//TRIM(Spt)//' lwsol: ',wrf_grid%lwsol(ip,jp), &
931	' qsol_rst: ',qsol_rst(lp)
932
933	PRINT *,' Lluis lmdz pt= ',lp,' k qsurf_rst qsol_rst evap_rst snow_rst ' // &
934	'fder_rst agesno_rst rugos_rst'
935	DO iz=1,4
936	PRINT ,' ',iz,qsurf_rst(lp,iz),qsol_rst(lp), evap_rst(lp,iz), &
937	snow_rst(lp,iz),fder_rst(lp), agesno_rst(lp,iz),rugos_rst(lp,iz)
938	PRINT *,' ftsoil: ',ftsoil_rst(lp,:,iz)
939	END DO
940
941	DO ix=1,ddimx
942	DO iy=1,ddimy
943	ixy=ddimx*(iy-1)+ix
944	IF (wrf_grid%lter(ix,iy) /= pctsrf(ixy,1)) PRINT *,' ' // &
945	TRIM(subname) // ' ter DOES NOT COINCIDE! ', ix,', ',iy,' lter: ', &
946	wrf_grid%lter(ix,iy),' pctsrf: ',pctsrf(ixy,1)
947	IF (wrf_grid%llic(ix,iy) /= pctsrf(ixy,2)) PRINT *,' ' // &
948	TRIM(subname) // ' lic DOES NOT COINCIDE! ', ix,', ',iy,' llic: ', &
949	wrf_grid%llic(ix,iy),' pctsrf: ',pctsrf(ixy,2)
950	IF (wrf_grid%loce(ix,iy) /= pctsrf(ixy,3)) PRINT *,' ' // &
951	TRIM(subname) // ' oce DOES NOT COINCIDE! ', ix,', ',iy,' loce: ', &
952	wrf_grid%loce(ix,iy),' pctsrf: ',pctsrf(ixy,3)
953	IF (wrf_grid%lsic(ix,iy) /= pctsrf(ixy,4)) PRINT *,' ' // &
954	TRIM(subname) // ' sic DOES NOT COINCIDE! ', ix,', ',iy,' lsic: ', &
955	wrf_grid%lsic(ix,iy),' pctsrf: ',pctsrf(ixy,4)
956	END DO
957	END DO
958
959	CALL NOread_limit_init(ddimxy, klon_glo)
960
961	PRINT *,' Lluis klon_glo,', klon_glo,' klon: ',klon,'lp= ',lp,' sst: ',sst(lp)
962	sst_glo = sst(1:ddimxy)
963	! In ocean_forced_mod is used as tsurf_limit = sst(knindex(1:knon))
964	tsurf_limit = ftsol(:,3)
965	PRINT *,' Lluis pt tsurf_limit_check_____'
966	DO ix=-10,10
967	PRINT *,lp+ix,tsurf_limit(lp+ix),sst(lp+ix)
968	END DO
969
970	! In surf_land_bucket_mod is used as z0_new=rugos(knindex(1:knon),1) [z0_new == z0_limit]
971	z0_limit = rugos_rst(:,1)
972	alb_limit = albedo
973	PRINT *,' Lluis pre-limit_values ftsol: ',ftsol(lp,:),' pctsrf: ',pctsrf(lp,:)
974	PRINT *,' Lluis limit_values lp= ',lp,' tsurf_limit: ',tsurf_limit(lp), &
975	' z0_limit: ', z0_limit(lp),' alb_limit: ',alb_limit(lp)
976
977	rug_glo = rugos_rst(:,is_ter)
978	pct_glo(klon_glo,:) = pctsrf(klon_glo,:)
979
980	! It is not clear how to deal with 'limit.nc' times...
981	! temps=wjulday
982
983	!! PRINT *,' Lluis: calling iniphysiq'
984	CALL iniphysiq(ddimx*ddimy, dimz, 86400., wjulday, wtime_step, &
985	& LatValues, LonValues, wrf_area, CUValues, CVValues, 1./reradius, g, r_d, cp, &
986	& iflag_phys)
987
988	! NO NEED to do this, using wrfinput_d[nn], wrfbdy_d[nn] and wrflowinput_d[nn]
989	! Checking for 'startphy.nc' presence
990	!! INQUIRE(FILE='startphy.nc', EXIST=existsfile)
991	!! IF (.NOT.existsfile) THEN
992	!! WRITE(wrf_err_message,*)" file 'startphy.nc' with initial conditions " // &
993	!! & ' for LMDZphys does not exist!!!!'
994	!! CALL wrf_error_fatal(TRIM(wrf_err_message))
995	!! END IF
996
997	CALL init_ovars_lmdz_NOmodule(ddimxy,dimz,nbsrf)
998
999	PRINT *,' Lluis: END of firststep !!!!!'
1000	END IF firststep
1001
1002	ip = 12
1003	jp = 7
1004	PRINT *,' ',wrf_grid%itimestep,' Lluis at 12 7: ',wrf_grid%lter(ip,jp), wrf_grid%llic(ip,jp), &
1005	wrf_grid%loce(ip,jp), wrf_grid%lsic(ip,jp), ' rugksoil: ', &
1006	wrf_grid%lrugksoil(ip,1,jp), wrf_grid%lrugksoil(ip,2,jp), &
1007	wrf_grid%lrugksoil(ip,3,jp), wrf_grid%lrugksoil(ip,4,jp),' rug: ', &
1008	wrf_grid%lrug(ip,jp)
1009	ip = 6
1010	jp = 9
1011	PRINT *,' ',wrf_grid%itimestep,' Lluis at 6 9: ',wrf_grid%lter(ip,jp), wrf_grid%llic(ip,jp), &
1012	wrf_grid%loce(ip,jp), wrf_grid%lsic(ip,jp), ' rugksoil: ', &
1013	wrf_grid%lrugksoil(ip,1,jp), wrf_grid%lrugksoil(ip,2,jp), &
1014	wrf_grid%lrugksoil(ip,3,jp), wrf_grid%lrugksoil(ip,4,jp),' rug: ', &
1015	wrf_grid%lrug(ip,jp)
1016
1017	ip=icheck_p
1018	jp=jcheck_p
1019	! WRF alarms
1020	! I can not print them !!!!
1021	!!
1022	! wrf_grid%alarms => wrf_alarms
1023	! PRINT *,' Lluis WRF alarms: time-step', wrf_grid%itimestep,' _______'
1024	! DO ix=1,3*INT(MAX_WRF_ALARMS/3),3
1025	! PRINT *,(ix+iy,': ',wrf_alarms(ix+iy), iy=0,2)
1026	! END DO
1027	! PRINT ,(iy,': ',wrf_alarms(iy), iy=3INT(MAX_WRF_ALARMS/3)+1, MAX_WRF_ALARMS)
1028
1029	! What to do with this?
1030	!!! IF ( wrf_restart_alarm ) lmdz_lafin = .TRUE.
1031
1032	INQUIRE(FILE='traceur.def', EXIST=existsfile)
1033	IF (existsfile) THEN
1034	DO funit=10,100
1035	INQUIRE(UNIT=funit, OPENED=is_used)
1036	IF (.NOT. is_used) EXIT
1037	END DO
1038	OPEN(UNIT=10, FILE='traceur.def', STATUS='OLD')
1039	READ(10,*)lmdzmixingratios
1040	CLOSE(UNIT=10)
1041	ELSE
1042
1043	WRITE(message, *)errmsg
1044	CALL wrf_debug(75, message)
1045	WRITE(wrf_err_message,*)' required file traceur.def not found !!'
1046	CALL wrf_error_fatal(TRIM(wrf_err_message))
1047	END IF
1048
1049	! Regarding if it is a lowbdy input time
1050	!!
1051	PRINT *,' ',wrf_grid%itimestep,' wrf_islowbdyin: ',wrf_islowbdyin
1052	! IF ((wrf_islowbdyin) .AND. .NOT.wrf_restart_alarm) THEN
1053	IF ((wrf_islowbdyin)) THEN
1054	PRINT *,' LOWBDY time-step!!!!!!! '
1055	PRINT *,' before: wsst ', wrf_grid%sst(ip,jp),' lsst: ',sst(lp)
1056	CALL get_lowbdy(ddimx, ddimy, dimz, wbdyw, ddimxy, 4, 11, klon_glo, &
1057	wrf_grid, pctsrf, ftsol, rugos_rst, albedo, sst, tsurf_limit, z0_limit, &
1058	alb_limit, rug_glo, pct_glo)
1059	PRINT *,' after: wsst ', wrf_grid%sst(ip,jp),' lsst: ',sst(lp)
1060	! STOP
1061	END IF
1062
1063	IF (ALLOCATED(MixingRatioValues)) DEALLOCATE(MixingRatioValues)
1064	ALLOCATE(MixingRatioValues(ddimxy,dimz,lmdzmixingratios), STAT=ierr)
1065	IF (ierr /= 0) THEN
1066	PRINT *,errmsg
1067	PRINT *,' ' // subname // ': Problem allocating WRF moisture array!', ierr
1068	END IF
1069
1070	!! PRINT *,' Lluis wrf MOIST:',LBOUND(wrf_MOIST),' x ',UBOUND(wrf_MOIST)
1071	!! PRINT *,'Mix ratio: ',LBOUND(MixingRatioValues),' x ', UBOUND(MixingRatioValues)
1072
1073	IF (ALLOCATED(MixingRatiotendValues)) DEALLOCATE(MixingRatiotendValues)
1074	ALLOCATE(MixingRatiotendValues(ddimxy,dimz,lmdzmixingratios))
1075
1076	! Transforming from WRF theta temperature to LMDZ temperature (in K)
1077	!!
1078	wrf_P = wrf_baseP + wrf_perP
1079	wrf_geopot = wrf_basegeopot + wrf_pergeopot
1080	wrf_sfcgeopot = wrf_grid%ht*g
1081	! de-staggering winds
1082	!!
1083
1084	wrfUdestagg(wims:wime-1,:,wjms:wjme-1) = (wrf_U(wims:wime-1,:,wjms:wjme-1) + &
1085	wrf_U(wims+1:wime,:,wjms:wjme-1))/2.
1086	wrfVdestagg(wims:wime-1,:,wjms:wjme-1) = (wrf_V(wims:wime-1,:,wjms:wjme-1) + &
1087	wrf_V(wims:wime-1,:,wjms+1:wjme))/2.
1088	wrfUtenddestagg(wims:wime-1,:,wjms:wjme-1) = &
1089	(wrf_Utend(wims:wime-1,:,wjms:wjme-1) + wrf_Utend(wims+1:wime,:,wjms:wjme-1))/2.
1090	wrfVtenddestagg(wims:wime-1,:,wjms:wjme-1) = &
1091	(wrf_Vtend(wims:wime-1,:,wjms:wjme-1) + wrf_Vtend(wims:wime-1,:,wjms+1:wjme))/2.
1092
1093	!! NO de-stagger in z
1094	!! wrfWValues = (wrf_w(:,wkms:wkme-1,:)+wrf_w(:,wkms+1:wkme,:))/2.
1095	!! wrfPERMASSValues = (wrf_permass(:,wkms:wkme-1,:)+wrf_permass(:,wkms+1:wkme,:))/2.
1096	!! DO k=wkms,wkme
1097	!! wrfMASSValues(:,k,:)=-(wrfPERMASSVAlues(:,k,:)+wrf_basemass)*wrf_dfulleta(k)/g
1098	!! END DO
1099	wrfWValues=wrf_w
1100	DO iz=wkms,wkme
1101	wrfMASSValues(:,iz,:) = -(wrf_permass+wrf_basemass)*wrf_dfulleta(iz)/g
1102	END DO
1103	wrfZFluxMassValues=wrfWValues*wrfMASSValues
1104
1105	!! CALL theta_to_temp(THETAVALUES=wrf_T, PRESSVALUES=wrf_P, TEMPVALUES=wrf_Temp &
1106	!! ! Dimension arguments
1107	!! & ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
1108	!! & ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
1109	!! & ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
1110	!! & ,I_START=wi_start,I_END=wi_end &
1111	!! & ,J_START=wj_start,J_END=wj_end &
1112	!! & ,KTS=wkts, KTE=wkte &
1113	!! & ,NUM_TILES=wnum_tiles &
1114	!! & )
1115
1116	CALL theta_to_temp(wrf_T, wrf_P, wrf_Temp &
1117	! Dimension arguments
1118	& ,wids,wide, wjds,wjde, wkds,wkde &
1119	& ,wims,wime, wjms,wjme, wkms,wkme &
1120	& ,wips,wipe, wjps,wjpe, wkps,wkpe &
1121	& ,wi_start,wi_end &
1122	& ,wj_start,wj_end &
1123	& ,wkts, wkte &
1124	& ,wnum_tiles &
1125	& )
1126
1127	!! CALL theta_to_temp(THETAVALUES=wrf_Ttend, PRESSVALUES=wrf_P, &
1128	!! & TEMPVALUES=wrf_Temptend &
1129	!! ! Dimension arguments
1130	!! & ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
1131	!! & ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
1132	!! & ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
1133	!! & ,I_START=wi_start,I_END=wi_end &
1134	!! & ,J_START=wj_start,J_END=wj_end &
1135	!! & ,KTS=wkts, KTE=wkte &
1136	!! & ,NUM_TILES=wnum_tiles &
1137	!! & )
1138
1139	CALL theta_to_temp(wrf_Ttend, wrf_P, wrf_Temptend &
1140	! Dimension arguments
1141	& ,wids,wide, wjds,wjde, wkds,wkde &
1142	& ,wims,wime, wjms,wjme, wkms,wkme &
1143	& ,wips,wipe, wjps,wjpe, wkps,wkpe &
1144	& ,wi_start,wi_end &
1145	& ,wj_start,wj_end &
1146	& ,wkts, wkte &
1147	& ,wnum_tiles &
1148	& )
1149
1150	DO ix=1,ddimx
1151	DO iy=1,ddimy
1152	sfcGeopotValues(ddimx(iy-1)+ix) = wrf_grid%ht(ix,iy)g
1153	END DO
1154	END DO
1155
1156	CALL mat_vect(wrf_Temp, dimx, dimy, dimz, wbdyw, TValues)
1157	CALL mat_vect(wrf_Temptend, dimx, dimy, dimz, wbdyw, TtendValues)
1158	CALL mat_vect(wrfUdestagg, dimx, dimy, dimz, wbdyw, UValues)
1159	CALL mat_vect(wrfUtenddestagg, dimx, dimy, dimz, wbdyw, UtendValues)
1160	CALL mat_vect(wrfVdestagg, dimx, dimy, dimz, wbdyw, VValues)
1161	CALL mat_vect(wrfVtenddestagg, dimx, dimy, dimz, wbdyw, VtendValues)
1162	! CALL mat_vect(wrf_U, dimx, dimy, dimz, wbdyw, UValues)
1163	! CALL mat_vect(wrf_Utend, dimx, dimy, dimz, wbdyw, UtendValues)
1164	! CALL mat_vect(wrf_V, dimx, dimy, dimz, wbdyw, VValues)
1165	! CALL mat_vect(wrf_Vtend, dimx, dimy, dimz, wbdyw, VtendValues)
1166	CALL mat_vect(wrf_P, dimx, dimy, dimz, wbdyw, HalfPressValues)
1167	CALL mat_vect_zstagg(wrf_fullpres, dimx, dimy, dimz+1, wbdyw, FullPressValues)
1168	CALL mat_vect_zstagg(wrf_geopot, dimx, dimy, dimz+1, wbdyw, GeopotValues)
1169	CALL mat_vect(wrfZFluxMassValues, dimx, dimy, dimz, wbdyw, WFluxMassValues)
1170
1171	CALL WRFmoist_LMDZmoist(wrf_MOIST, dimx, dimy, dimz, wnum3dm, wbdyw, wrf_qvid, &
1172	& wrf_qcid, wrf_qrid, lmdzmixingratios, MixingRatioValues)
1173
1174	CALL WRFmoist_LMDZmoist(wrf_MOISTtend, dimx, dimy, dimz, wnum3dm, wbdyw, wrf_qvid,&
1175	& wrf_qcid, wrf_qrid, lmdzmixingratios, MixingRatiotendValues)
1176
1177	preff=100000.
1178	CALL eta_to_pressure(wrf_halfeta, preff, wrf_Ptop, dimz, lmdzoutP)
1179
1180	! We need to provide values to these variables in order to make it works
1181	presnivs=lmdzoutP
1182
1183	! Lluis Testing water species
1184	MixingRatioValues(:,:,2)=0.
1185
1186	IF (wrf_dbg >= 200) THEN
1187
1188	!! Values at a given LMDZ matrix point
1189	lmdzp=lp
1190	jlmdz=INT(lmdzp/ddimx) + 1
1191	ilmdz=lmdzp-(jlmdz-1)*ddimx
1192	!! Values at a given WRF 2D point
1193	! ilmdz=ip
1194	! jlmdz=jp
1195	! lmdzp=ddimx*(jlmdz-1)+ilmdz
1196
1197	PRINT *,' WRF+LMDZ lmdzp: ',lmdzp,' : ',ilmdz,jlmdz
1198	PRINT *,'k wrf_T wrf_Temp TValues_______________'
1199	DO iz=1,dimz
1200	PRINT *,iz, wrf_T(ilmdz,iz,jlmdz), wrf_Temp(ilmdz,iz,jlmdz), TValues(lmdzp,iz)
1201	END DO
1202	PRINT *,'k wrf_U Uvalues wrf_V VValues _______________'
1203	DO iz=1,dimz
1204	PRINT *,iz, wrf_U(ilmdz,iz,jlmdz), UValues(lmdzp,iz), wrf_V(ilmdz,iz,jlmdz), &
1205	VValues(lmdzp,iz)
1206	END DO
1207	PRINT *,'k wrf_fullpres FUllPressValues wrf_sfcgeopot hgt sfcGeopotValues______'
1208	DO iz=1,dimz+1
1209	PRINT *,iz,wrf_fullpres(ilmdz,iz,jlmdz), FullPressValues(lmdzp,iz), &
1210	wrf_sfcgeopot(ilmdz,jlmdz), wrf_grid%ht(ilmdz,jlmdz), sfcGeopotValues(lmdzp)
1211	END DO
1212	PRINT *,'k wrf_halfpress HalfPressValues wrf_Geopot GeopotValues lmdz outP ' //&
1213	'_______________'
1214	DO iz=1,dimz
1215	PRINT *,iz, wrf_P(ilmdz,iz,jlmdz), HalfPressValues(lmdzp,iz), &
1216	wrf_geopot(ilmdz,iz,jlmdz), GeopotValues(lmdzp,iz), lmdzoutP(iz)
1217	END DO
1218	PRINT *,'k wrf_MOIST MixingRatioValues _______________'
1219	DO iz=1,dimz
1220	PRINT *,iz,wrf_MOIST(ilmdz,iz,jlmdz,wrf_qvid), MixingRatioValues(lmdzp,iz,1),&
1221	wrf_MOIST(ix,iz,jlmdz,wrf_qcid), MixingRatioValues(lmdzp,iz,2)
1222	END DO
1223	PRINT *,'k wrf_mass wrf_W wrf_dfulleta wrfMASSValues wrfZFluxMassValues ' // &
1224	'WFluxMassValues ____________'
1225	DO iz=1,dimz
1226	PRINT *,iz,wrf_permass(ilmdz,jlmdz)+wrf_basemass(ilmdz,jlmdz), &
1227	wrf_W(ilmdz,iz,jlmdz), wrf_dfulleta(iz), wrfMASSValues(ilmdz,iz,jlmdz), &
1228	wrfZFluxMassValues(ilmdz,iz,jlmdz), WFluxMassValues(lmdzp,iz)
1229	END DO
1230	END IF
1231
1232	PRINT *,' Lluis before physiq '//TRIM(Spt)//' ltksoil: ', &
1233	wrf_grid%ltksoil(ip,:,jp), ' ftsol: ',ftsol(lp,:)
1234
1235	PRINT *,' Lluis beofre physiq '//TRIM(Spt)//' qsurf: ',zxqsurf(lp), &
1236	' qsurf_rst: ',qsurf_rst(lp,:),' smois :',wrf_grid%smois(ip,:,jp), &
1237	' lqksoil: ', wrf_grid%lqksoil(ip,:,jp)
1238
1239	PRINT *,' Lluis before physiq: '//TRIM(Spt)//' FullPressValues HalfPressValues '&
1240	// 'GeopotValues UValues VValues TValues: '
1241	DO iz=1,dimz
1242	PRINT *,FullPressValues(lp,iz), HalfPressValues(lp,iz), GeopotValues(lp,iz), &
1243	UValues(lp,iz), VValues(lp,iz), TValues(lp,iz)
1244	END DO
1245
1246	ip = 23
1247	jp = 18
1248	PRINT *,' ',wrf_grid%itimestep,' Lluis before physiq at ',ip,', ', jp,' sst: ', &
1249	wrf_grid%sst(ip,jp), sst((jp-1)ddimx+ip),' ftsol 3: ',ftsol((jp-1)ddimx+ip,is_oce)
1250
1251	ip=icheck_p
1252	jp=jcheck_p
1253
1254	!! L. Fita, July 2013.
1255	!! Misteriously 'physiq' suborutine from LMDZ can not be called using keys
1256	!!
1257	!! CALL physiq(NLON=dimxy, NLEV=dimz, &
1258	!! & DEBUT=lmdz_debut, LAFIN=lmdz_lafin, &
1259	!! & JD_CUR=wjulday, JH_CUR=wgmt, PDTPHYS=wtime_step, &
1260	!! & PAPRS=FullPressValues, PPLAY=HalfPressValues, &
1261	!! & PPHI=GeopotValues, PPHIS=wrf_sfcgeopot, &
1262	!! & PRESNIVS=lmdzoutP, CLESPHY0=phykeys, &
1263	!! & U=UValues, V=VValues, T=TValues, QX=MixingRatioValues, &
1264	!! & FLXMASS_W=WFluxMassValues, &
1265	!! & D_U=UtendValues, D_V=VtendValues, D_T=TtendValues, &
1266	!! & D_QX=MixingRatiotendValues, D_PS=wrf_psfctend, &
1267	!! & DUDYN=lmdz_dudyn, PVTETA=lmdz_PVtheta)
1268
1269	!! PRINT *,' Lluis: going to physiq'
1270
1271	CALL physiq(ddimxy, dimz, &
1272	& lmdz_debut, lmdz_lafin, &
1273	& wjulday, wgmt, wtime_step, &
1274	& FullPressValues, HalfPressValues, &
1275	& GeopotValues, sfcGeopotValues, &
1276	& lmdzoutP, phykeys, &
1277	& UValues, VValues, TValues, MixingRatioValues, &
1278	& WFluxMassValues, &
1279	& UtendValues, VtendValues, TtendValues, &
1280	& MixingRatiotendValues, psfctendValues, &
1281	& lmdz_dudyn, lmdz_PVtheta)
1282
1283	PRINT *,' Lluis after physiq '//TRIM(Spt)//' ltksoil: ', &
1284	wrf_grid%ltksoil(ip,:,jp), ' ftsol: ',ftsol(315,:)
1285	PRINT *,' Lluis after physiq '//TRIM(Spt)//' qsurf: ',zxqsurf(lp), &
1286	' qsurf_rst: ', qsurf_rst(lp,:),' smois :',wrf_grid%smois(ip,:,jp), &
1287	' lqksoil: ', wrf_grid%lqksoil(ip,:,jp)
1288	PRINT *,' Lluis after physiq: '//TRIM(Spt)//' FullPressValues HalfPressValues ' &
1289	// 'GeopotValuesUValues VValues TValues: '
1290	DO iz=1,dimz
1291	PRINT *,FullPressValues(lp,iz), HalfPressValues(lp,iz), GeopotValues(lp,iz), &
1292	UValues(lp,iz), VValues(lp,iz), TValues(lp,iz)
1293	END DO
1294
1295	ip = 23
1296	jp = 18
1297	PRINT *,' ',wrf_grid%itimestep,' Lluis after physiq at ',ip,', ', jp,' sst: ', &
1298	wrf_grid%sst(ip,jp), sst((jp-1)ddimx+ip),' ftsol 3: ',ftsol((jp-1)ddimx+ip,is_oce)
1299
1300	ip=icheck_p
1301	jp=jcheck_p
1302	! Lluis Testing water species
1303	!! MixingRatiotendValues(:,:,1)=1000.*MixingRatiotendValues(:,:,1)
1304	!! MixingRatiotendValues(:,:,1)=MixingRatiotendValues(:,:,1)+MixingRatiotendValues(:,:,2)
1305	MixingRatiotendValues(:,:,2)=0.
1306
1307	!! PRINT *,' Lluis: outgoing from physiq'
1308
1309	!! Tvalues, [U/V/T]Values, are not given back just inputs
1310	CALL vect_mat(TtendValues, dimx, dimy, dimz, wbdyw, wrf_Temptend)
1311	CALL vect_mat(UtendValues, dimx, dimy, dimz, wbdyw, wrfUtenddestagg)
1312	CALL vect_mat(VtendValues, dimx, dimy, dimz, wbdyw, wrfVtenddestagg)
1313	! CALL vect_mat(UtendValues, dimx, dimy, dimz, wbdyw, wrf_Utend)
1314	! CALL vect_mat(VtendValues, dimx, dimy, dimz, wbdyw, wrf_Vtend)
1315	!! CALL vect_mat(TValues, dimx, dimy, dimz, wbdyw, wrf_Temp)
1316	!! CALL vect_mat(UValues, dimx, dimy, dimz, wbdyw, wrf_U)
1317	!! CALL vect_mat(VValues, dimx, dimy, dimz, wbdyw, wrf_V)
1318	!! CALL vect_mat(HalfPressValues, dimx, dimy, dimz, wbdyw, wrf_P)
1319	!! CALL vect_mat(FullPressValues, dimx, dimy, dimz, wbdyw, wrf_fullpres)
1320	!! CALL vect_mat(GeopotValues, dimx, dimy, dimz, wbdyw, wrf_geopot)
1321
1322	DO iy=1,ddimy
1323	DO ix=1,ddimx
1324	wrf_psfctend(ix,iy) = PsfctendValues(ddimx*(iy-1)+ix)
1325	END DO
1326	END DO
1327
1328	!! No changes are suffered not need to give it back
1329	!! CALL temp_to_theta(TEMPVALUES=wrf_Temp, PRESSVALUES=wrf_P, THETAVALUES=wrf_T &
1330	!! ! Dimension arguments
1331	!! & ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
1332	!! & ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
1333	!! & ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
1334	!! & ,I_START=wi_start,I_END=wi_end &
1335	!! & ,J_START=wj_start,J_END=wj_end &
1336	!! & ,KTS=wkts, KTE=wkte &
1337	!! & ,NUM_TILES=wnum_tiles &
1338	!! & )
1339
1340	CALL temp_to_theta(TEMPVALUES=wrf_Temptend, PRESSVALUES=wrf_P, &
1341	& THETAVALUES=wrf_Ttend &
1342	! Dimension arguments
1343	& ,IDS=wids,IDE=wide, JDS=wjds,JDE=wjde, KDS=wkds,KDE=wkde &
1344	& ,IMS=wims,IME=wime, JMS=wjms,JME=wjme, KMS=wkms,KME=wkme &
1345	& ,IPS=wips,IPE=wipe, JPS=wjps,JPE=wjpe, KPS=wkps,KPE=wkpe &
1346	& ,I_START=wi_start,I_END=wi_end &
1347	& ,J_START=wj_start,J_END=wj_end &
1348	& ,KTS=wkts, KTE=wkte &
1349	& ,NUM_TILES=wnum_tiles &
1350	& )
1351
1352	wrf_Ttend = wrf_Ttend + 300.
1353
1354	! Re-staggering wind tendencies
1355	!!
1356
1357	wrf_Utend(wims,:,wjms:wjme-1) = wrfUtenddestagg(wims,:,wjms:wjme-1)
1358	wrf_Vtend(wims:wime-1,:,wjms) = wrfVtenddestagg(wims:wime-1,:,wjms)
1359	wrf_Utend(wims,:,wjme-1) = wrfUtenddestagg(wims,:,wjme-1)
1360	wrf_Vtend(wime-1,:,wjms) = wrfVtenddestagg(wime-1,:,wjms)
1361	wrf_Utend(wims+1:wime-1,:,wjms:wjme-1) = &
1362	& (wrfUtenddestagg(wims:wime-2,:,wjms:wjme-1) + &
1363	& wrfUtenddestagg(wims+1:wime-1,:,wjms:wjme-1))/2.
1364	wrf_Vtend(wims:wime-1,:,wjms+1:wjme-1) = &
1365	& (wrfVtenddestagg(wims:wime-1,:,wjms:wjme-2) + &
1366	& wrfVtenddestagg(wims:wime-1,:,wjms+1:wjme-1))/2.
1367	wrf_Utend(wime,:,wjms:wjme-1) = wrfUtenddestagg(wime-1,:,wjms:wjme-1)
1368	wrf_Vtend(wims:wime-1,:,wjme) = wrfVtenddestagg(wims:wime-1,:,wjme-1)
1369	wrf_Utend(wime,:,wjme) = wrfUtenddestagg(wime-1,:,wjme-1)
1370	wrf_Vtend(wime,:,wjme) = wrfVtenddestagg(wime-1,:,wjme-1)
1371
1372	PRINT *,' Lluis ' // TRIM(Spt) // ' wrf_Utend: ',wrf_Utend(ip,1,jp), &
1373	' wrfUtenddestagg a: ', wrfUtenddestagg(ip-1,1,jp),' wrfUtenddestagg b: ', &
1374	wrfUtenddestagg(ip,1,jp), ' wrfUtenddestagg c: ',wrfUtenddestagg(ip+1,1,jp), &
1375	' UtendValues: ', UtendValues(lp,1),' lp: ',lp
1376	PRINT *,'dim UtendValues: ',1-wbdyw,':',dimx-1-wbdyw+1,', ',dimz+1,', ',1-wbdyw,':',dimy-1-wbdyw+1
1377	!! wrf_Ttend = wrf_Temptend
1378
1379	!! wrf_perP = wrf_P - wrf_baseP
1380	!! wrf_pergeopot = wrf_geopot - wrf_basegeopot
1381
1382	!! CALL LMDZmoist_WRFmoist(MixingRatioValues, dimx, dimy, dimz, wnum3dm, wrf_qvid, &
1383	!! & wrf_qcid, wrf_qrid, lmdzmixingratios, wrf_MOIST)
1384
1385	CALL LMDZmoist_WRFmoist(MixingRatiotendValues, dimx, dimy, dimz, wnum3dm, wbdyw, &
1386	& wrf_qvid, wrf_qcid, wrf_qrid, lmdzmixingratios, wrf_MOISTtend)
1387
1388	IF (wrf_dbg >= 200) THEN
1389	PRINT *,' WRF+LMDZ lmdzp: ',lmdzp
1390	PRINT *,'k T tendencies _______________'
1391	DO iz=1,dimz
1392	PRINT *,iz, TtendValues(lmdzp,iz), wrf_Temptend(ilmdz,iz,jlmdz), &
1393	wrf_Ttend(ilmdz,iz,jlmdz)
1394	END DO
1395	PRINT *,'k U tendencies V tendencies_______________'
1396	DO iz=1,dimz
1397	PRINT *,iz, UtendValues(lmdzp,iz), wrf_Utend(ilmdz,iz,jlmdz), &
1398	VtendValues(lmdzp,iz), wrf_Vtend(ilmdz,iz,jlmdz)
1399	END DO
1400	PRINT *,'k MixingRatiotendValues wrf_MOISTtend _______________'
1401	DO iz=1,dimz
1402	PRINT *,iz, MixingRatiotendValues(lmdzp,iz,1), &
1403	wrf_MOISTtend(ilmdz,iz,jlmdz,wrf_qvid), MixingRatiotendValues(lmdzp,iz,2), &
1404	wrf_MOISTtend(ilmdz,iz,jlmdz,wrf_qcid)
1405	END DO
1406	PRINT *,'wrf_sfcgeopot psfc_tendencies sfcGeopotValues_______________'
1407	PRINT *,wrf_sfcgeopot(ilmdz,jlmdz), wrf_psfctend(ilmdz,jlmdz), &
1408	sfcGeopotValues(lmdzp)
1409	END IF
1410
1411	PRINT *,' Lluis psfctendValues ______________ '
1412	DO ix=1,4*INT(ddimxy/4),4
1413	PRINT *,(ix+iy,':',psfctendValues(ix+iy), iy=0,3)
1414	END DO
1415	PRINT ,(iy,':',psfctendValues(iy), iy=4INT(ddimxy/4)+1,ddimxy)
1416
1417	! All the tendencies have to be given as decoupled from the dry air mass column! (see dyn_em/module_em.F:rk_update_scalar)
1418	!!
1419	DO iz=1,dimz
1420	wrf_Ttend(:,iz,:) = wrf_Ttend(:,iz,:)*wrf_mut
1421	wrf_Utend(:,iz,:) = wrf_Utend(:,iz,:)*wrf_muu
1422	wrf_Vtend(:,iz,:) = wrf_Vtend(:,iz,:)*wrf_muv
1423	DO iq=1,wnum3dm
1424	wrf_MOISTtend(:,iz,:,iq) = wrf_MOISTtend(:,iz,:,iq)*wrf_mut
1425	END DO
1426	END DO
1427
1428	CALL get_lmdz(ddimx,ddimy,dimz,wbdyw,ddimxy,4,11,rlon,rlat,zmasq, &
1429	pctsrf,ftsol,ftsoil,zxqsurf,qsurf,falb1,falb2,fevap,zxsnow,radsol,solsw, &
1430	sollw,fder,rain_fall,snow_fall,agesno,zmea,zstd,zgam,zthe,zpic,zval,rugoro, &
1431	frugs,zmax0,f0,wake_s,wake_cstar,wake_pe,wake_fip,sst,albedo, &
1432	t_ancien,q_ancien,u_ancien,v_ancien,clwcon,rnebcon,ratqs,ema_work1,ema_work2, &
1433	wake_deltat,wake_deltaq,fm_therm,entr_therm,detr_therm,phy_glo,wrf_grid)
1434	! CALL fonte_neige_final(runoff_lic)
1435
1436	! Some re-start variables
1437	PRINT *,' Lluis lrunofflic: ',LBOUND(wrf_grid%lrunofflic), ':', &
1438	UBOUND(wrf_grid%lrunofflic),' ruonff_lic: ',UBOUND(runoff_lic),' ddimx: ', &
1439	ddimx,' ddimy: ',ddimy
1440
1441	! PRINT ,' iy ix ddimx(iy-1)+ix runoff_lic wrf_grid%lrunofflic_______'
1442	! DO iy=1,ddimy
1443	! DO ix=1,ddimx
1444	! PRINT ,iy,ix,ddimx(iy-1)+ix
1445	! wrf_grid%lrunofflic(ix,iy) = runoff_lic(ddimx*(iy-1)+ix)
1446	! PRINT ,iy,ix,ddimx(iy-1)+ix,runoff_lic(ddimx*(iy-1)+ix),wrf_grid%lrunofflic(ix,iy)
1447	! END DO
1448	! END DO
1449
1450	! STOP
1451	!! CALL get_lmdz(ddimx,ddimy,dimz,wbdyw,ddimxy,4,11,zmasq,pctsrf,ftsol, &
1452	!! ftsoil,fqsurf,qsol_d,falb1,falb2,fevap,zxsnow,radsol,solsw,sollw,fder, &
1453	!! rain_fall,snow_fall,agesno,zmea,zstd,zgam,zthe,zpic,zval,rugoro, &
1454	!! frugs,run_off_lic,zmax0,f0,wake_s,wake_cstar,wake_pe,wake_fip,sst, &
1455	!! albedo,t_ancien,q_ancien,u_ancien,v_ancien,clwcon,rnebcon,ratqs, &
1456	!! ema_work1,ema_work2,wake_deltat,wake_deltaq,fm_therm,entr_therm, &
1457	!! detr_therm,phy_glo,wrf_grid)
1458
1459	PRINT *,' Lluis after get '//TRIM(Spt)//' ltksoil: ',wrf_grid%ltksoil(ip,:,jp), &
1460	' ftsol: ',ftsol(lp,:)
1461
1462	PRINT *,' Lluis after get '//TRIM(Spt)//' fqsurf: ',zxqsurf(lp),' qsurf_rst: ', &
1463	qsurf_rst(lp,:),' smois :',wrf_grid%smois(ip,:,jp),' lqksoil: ', &
1464	wrf_grid%lqksoil(ip,:,jp)
1465
1466	! Recomputig again the total roughness length
1467	LMDZvarmethod = 'prod'
1468
1469	PRINT *,' AFTER physic and before recomputing!!!!!'
1470	ip = 12
1471	jp = 7
1472	PRINT *,' ',wrf_grid%itimestep,' Lluis at 12 7: ',wrf_grid%lter(ip,jp), wrf_grid%llic(ip,jp), &
1473	wrf_grid%loce(ip,jp), wrf_grid%lsic(ip,jp), ' rugksoil: ', &
1474	wrf_grid%lrugksoil(ip,1,jp), wrf_grid%lrugksoil(ip,2,jp), &
1475	wrf_grid%lrugksoil(ip,3,jp), wrf_grid%lrugksoil(ip,4,jp),' rug: ', &
1476	wrf_grid%lrug(ip,jp)
1477	ip = 6
1478	jp = 9
1479	PRINT *,' ',wrf_grid%itimestep,' Lluis at 6 9: ',wrf_grid%lter(ip,jp), wrf_grid%llic(ip,jp), &
1480	wrf_grid%loce(ip,jp), wrf_grid%lsic(ip,jp), ' rugksoil: ', &
1481	wrf_grid%lrugksoil(ip,1,jp), wrf_grid%lrugksoil(ip,2,jp), &
1482	wrf_grid%lrugksoil(ip,3,jp), wrf_grid%lrugksoil(ip,4,jp),' rug: ', &
1483	wrf_grid%lrug(ip,jp)
1484
1485	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
1486	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
1487	wrf_grid%lrugksoil(:,1,:),wrf_grid%lrugksoil(:,2,:),wrf_grid%lrugksoil(:,3,:), &
1488	wrf_grid%lrugksoil(:,4,:),wrf_grid%lrug)
1489
1490	PRINT *,' AFTER physic and recomputing!!!!!'
1491	ip = 12
1492	jp = 7
1493	PRINT *,' ',wrf_grid%itimestep,' Lluis at 12 7: ',wrf_grid%lter(ip,jp), wrf_grid%llic(ip,jp), &
1494	wrf_grid%loce(ip,jp), wrf_grid%lsic(ip,jp), ' rugksoil: ', &
1495	wrf_grid%lrugksoil(ip,1,jp), wrf_grid%lrugksoil(ip,2,jp), &
1496	wrf_grid%lrugksoil(ip,3,jp), wrf_grid%lrugksoil(ip,4,jp),' rug: ', &
1497	wrf_grid%lrug(ip,jp)
1498	ip = 6
1499	jp = 9
1500	PRINT *,' ',wrf_grid%itimestep,' Lluis at 6 9: ',wrf_grid%lter(ip,jp), wrf_grid%llic(ip,jp), &
1501	wrf_grid%loce(ip,jp), wrf_grid%lsic(ip,jp), ' rugksoil: ', &
1502	wrf_grid%lrugksoil(ip,1,jp), wrf_grid%lrugksoil(ip,2,jp), &
1503	wrf_grid%lrugksoil(ip,3,jp), wrf_grid%lrugksoil(ip,4,jp),' rug: ', &
1504	wrf_grid%lrug(ip,jp)
1505
1506	ip = 23
1507	jp = 18
1508	PRINT *,' ',wrf_grid%itimestep,' Lluis at ',ip,', ', jp,' sst : ', &
1509	wrf_grid%sst(ip,jp), sst((jp-1)ddimx+ip),' ftsol 3: ',ftsol((jp-1)ddimx+ip,is_oce)
1510
1511	CALL wrf_varKsoil(wims,wime,wjms,wjme,dimx,dimy,wbdyw,LMDZvarmethod, &
1512	wrf_grid%lter,wrf_grid%llic,wrf_grid%loce,wrf_grid%lsic, &
1513	wrf_grid%lalbksoil(:,1,:),wrf_grid%lalbksoil(:,2,:),wrf_grid%lalbksoil(:,3,:), &
1514	wrf_grid%lalbksoil(:,4,:),wrf_grid%albedo)
1515
1516	! L. Fita, LMD January 2014.
1517	! It should work nicely in this way... but memory got crazy...
1518	CALL get_lmdz_out(ddimx, ddimy, dimz, wbdyw, ddimxy, nbsrf, nsoilmx, &
1519	& lmdzmixingratios, 1, 2, wrf_l_pbl, wrf_l_con, wrf_l_thermals, wrf_l_wake, &
1520	& 0, wrf_grid%itimestep, wtime_step, FullPressValues, HalfPressValues, &
1521	& GeopotValues, sfcGeopotValues, lmdzoutP, MixingRatioValues, &
1522	& UtendValues, VtendValues, TtendValues, MixingRatiotendValues, &
1523	& psfctendValues, wrf_grid)
1524
1525	ip=icheck_p
1526	jp=jcheck_p
1527
1528	CALL put_lmdz_in_WRFout(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1529	& wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, &
1530	& lmdzmixingratios, 1, 2, wrf_l_pbl, wrf_l_con, wrf_l_thermals, wrf_l_wake, &
1531	& wrf_grid%itimestep, wtime_step, FullPressValues, HalfPressValues, &
1532	& GeopotValues, sfcGeopotValues, lmdzoutP, MixingRatioValues, &
1533	& UtendValues, VtendValues, TtendValues, MixingRatiotendValues, &
1534	& psfctendValues, ORCHIDEElevels, ftsoil, WRF_nsoillayers, &
1535	& wrf_grid%soil_layers/100., wrf_grid%zs, wrf_grid%dzs, wrf_grid%nest_pos, &
1536	& wrf_grid%q2, wrf_grid%t2, wrf_grid%th2, wrf_grid%psfc, wrf_grid%u10, &
1537	& wrf_grid%v10, wrf_grid%resm, wrf_grid%zetatop, wrf_grid%tslb,wrf_grid%smois,&
1538	& wrf_grid%sh2o, wrf_grid%smcrel, wrf_grid%sfcrunoff, wrf_grid%udrunoff, &
1539	& wrf_grid%grdflx, wrf_grid%acgrdflx, wrf_grid%snow, wrf_grid%snowh, &
1540	& wrf_grid%canwat, wrf_grid%sst, wrf_grid%sstsk, wrf_grid%lai, &
1541	& wrf_grid%h_diabatic, wrf_grid%sina, wrf_grid%tsk, wrf_grid%tiso, &
1542	& wrf_grid%max_msftx, wrf_grid%max_msfty, wrf_grid%rainc, wrf_grid%raincv, &
1543	& wrf_grid%rainsh, wrf_grid%rainshv, wrf_grid%rainnc, wrf_grid%rainncv, &
1544	& wrf_grid%snownc, wrf_grid%snowncv, wrf_grid%graupelnc, wrf_grid%graupelncv, &
1545	& wrf_grid%hailnc, wrf_grid%hailncv, wrf_grid%stepave_count, wrf_grid%cldfra, &
1546	& wrf_grid%swdown, wrf_grid%glw, wrf_grid%swnorm, wrf_grid%olr, &
1547	& wrf_grid%emiss, wrf_grid%tmn, wrf_grid%ust, wrf_grid%pblh, wrf_grid%hfx, &
1548	& wrf_grid%qfx, wrf_grid%lh, wrf_grid%achfx, wrf_grid%aclhf, &
1549	& wrf_grid%snowc, wrf_grid%sr, wrf_grid%potevp, wrf_grid%snopcx, &
1550	& wrf_grid%soiltb)
1551	! & wrf_grid%save_topo_from_real, &
1552	! & wrf_grid%avg_fuel_frac, wrf_grid%uah, wrf_grid%vah, wrf_grid%seed1, &
1553	! & wrf_grid%seed2)
1554
1555	ip = 23
1556	jp = 18
1557	PRINT *,' ',wrf_grid%itimestep,' Lluis after put at ',ip,', ', jp,' sst: ', &
1558	wrf_grid%sst(ip,jp), sst((jp-1)ddimx+ip),' ftsol 3: ',ftsol((jp-1)ddimx+ip,is_oce)
1559
1560	ip=icheck_p
1561	jp=jcheck_p
1562
1563
1564	! So it is done in a more larger to type and less elegant way...
1565	!
1566	! CALL get_lmdz_out2D_i(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1567	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
1568	! & wrf_l_wake, wrf_grid%itimestep, wtime_step, FullPressValues2, &
1569	! & HalfPressValues2, GeopotValues2, &
1570	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1571	! & UtendValues2, VtendValues2, TtendValues2, &
1572	! & MixingRatiotendValues2, psfctendValues2, wrf_grid%lages_ter, &
1573	! & wrf_grid%lages_lic, wrf_grid%lages_oce, wrf_grid%lages_sic, wrf_grid%laire, &
1574	! & wrf_grid%laireter, wrf_grid%lalb1, wrf_grid%lalb2, wrf_grid%lalbe_ter, &
1575	! & wrf_grid%lalbe_lic, wrf_grid%lalbe_oce, wrf_grid%lalbe_sic, wrf_grid%lale, &
1576	! & wrf_grid%lale_bl, wrf_grid%lale_wk, wrf_grid%lalp, wrf_grid%lalp_bl, &
1577	! & wrf_grid%lalp_wk, wrf_grid%lbils, wrf_grid%lbils_diss, wrf_grid%lbils_ec, &
1578	! & wrf_grid%lbils_enthalp, wrf_grid%lbils_kinetic, wrf_grid%lbils_latent, &
1579	! & wrf_grid%lbils_tke, wrf_grid%lcape, wrf_grid%lcape_max, wrf_grid%lcdrh, &
1580	! & wrf_grid%lcdrm, wrf_grid%lcldh, wrf_grid%lcldl, wrf_grid%lcldm, &
1581	! & wrf_grid%lcldq, wrf_grid%lcldt, wrf_grid%lcontfracatm, wrf_grid%lcontfracor,&
1582	! & wrf_grid%ldthmin, wrf_grid%ldtsvdft, wrf_grid%ldtsvdfo, wrf_grid%ldtsvdfg, &
1583	! & wrf_grid%ldtsvdfi)
1584
1585	! CALL get_lmdz_out2D_ii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1586	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
1587	! & wrf_l_wake, wrf_grid%itimestep, wtime_step, FullPressValues2, &
1588	! & HalfPressValues2, GeopotValues2, &
1589	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1590	! & UtendValues2, VtendValues2, TtendValues2, &
1591	! & MixingRatiotendValues2, psfctendValues2, wrf_grid%levap, wrf_grid%levap_ter,&
1592	! & wrf_grid%levap_lic, &
1593	! & wrf_grid%levap_oce, wrf_grid%levap_sic, wrf_grid%levappot_ter, &
1594	! & wrf_grid%levappot_lic, wrf_grid%levappot_oce, wrf_grid%levappot_sic, &
1595	! & wrf_grid%lfbase, wrf_grid%lffonte, wrf_grid%lflat, wrf_grid%lflw_ter, &
1596	! & wrf_grid%lflw_lic, wrf_grid%lflw_oce, wrf_grid%lflw_sic, &
1597	! & wrf_grid%lfqcalving, wrf_grid%lfqfonte, wrf_grid%lfract_ter, &
1598	! & wrf_grid%lfract_lic, wrf_grid%lfract_oce, wrf_grid%lfract_sic, &
1599	! & wrf_grid%lfsnow, wrf_grid%lfsw_ter, wrf_grid%lfsw_lic, wrf_grid%lfsw_oce, &
1600	! & wrf_grid%lfsw_sic, wrf_grid%lftime_th, wrf_grid%liwp, wrf_grid%llat_ter, &
1601	! & wrf_grid%llat_lic, wrf_grid%llat_oce, wrf_grid%llat_sic, wrf_grid%llmaxth, &
1602	! & wrf_grid%llwp, wrf_grid%lmsnow)
1603
1604	! CALL get_lmdz_out2D_iii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1605	! & wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, &
1606	! & wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
1607	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1608	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1609	! & UtendValues2, VtendValues2, TtendValues2, &
1610	! & MixingRatiotendValues2, psfctendValues2,wrf_grid%lndayrain,wrf_grid%lnettop,&
1611	! & wrf_grid%lpbase, wrf_grid%lphis, wrf_grid%lplcl, wrf_grid%lplfc, &
1612	! & wrf_grid%lpluc, wrf_grid%lplul, wrf_grid%lplulst, wrf_grid%lplulth, &
1613	! & wrf_grid%lpourc_ter, wrf_grid%lpourc_lic, wrf_grid%lpourc_oce, &
1614	! & wrf_grid%lpourc_sic, wrf_grid%lprecip, wrf_grid%lprw, wrf_grid%lpsol, &
1615	! & wrf_grid%lptop, wrf_grid%lqsat2m, wrf_grid%lqsol, wrf_grid%lqsurf, &
1616	! & wrf_grid%lradsol, wrf_grid%lrh2m, wrf_grid%lrh2m_max, wrf_grid%lrh2m_min, &
1617	! & wrf_grid%lrugs, wrf_grid%lrugs_ter, wrf_grid%lrugs_lic, wrf_grid%lrugs_oce, &
1618	! & wrf_grid%lrugs_sic, wrf_grid%lsens, wrf_grid%lsicf, wrf_grid%ls_lcl, &
1619	! & wrf_grid%lslp, wrf_grid%lsnow, wrf_grid%lsnowl, wrf_grid%lsoll, &
1620	! & wrf_grid%lsoll0, wrf_grid%lsollwdown, wrf_grid%lsols, wrf_grid%lsols0, &
1621	! & wrf_grid%ls_pblh, wrf_grid%ls_pblt, wrf_grid%ls_therm)
1622
1623	! CALL get_lmdz_out2D_iv(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1624	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
1625	! & wrf_l_wake, wrf_grid%itimestep, &
1626	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1627	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1628	! & UtendValues2, VtendValues2, TtendValues2, &
1629	! & MixingRatiotendValues2, psfctendValues2, wrf_grid%lt2m, &
1630	! & wrf_grid%lt2m_max, wrf_grid%lt2m_min, wrf_grid%lt2m_ter, wrf_grid%lt2m_lic, &
1631	! & wrf_grid%lt2m_oce, wrf_grid%lt2m_sic, wrf_grid%ltaux, wrf_grid%ltauy, &
1632	! & wrf_grid%ltaux_ter, wrf_grid%ltaux_lic, wrf_grid%ltaux_oce, &
1633	! & wrf_grid%ltaux_sic, wrf_grid%ltauy_ter, wrf_grid%ltauy_lic, &
1634	! & wrf_grid%ltauy_oce, wrf_grid%ltauy_sic, wrf_grid%lt_oce_sic, wrf_grid%ltopl,&
1635	! & wrf_grid%ltopl0, wrf_grid%ltops, wrf_grid%ltops0, wrf_grid%ltpot, &
1636	! & wrf_grid%ltpote, wrf_grid%ltsol, wrf_grid%ltsol_ter, wrf_grid%ltsol_lic, &
1637	! & wrf_grid%ltsol_oce, wrf_grid%ltsol_sic)
1638
1639	! CALL get_lmdz_out2D_v(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1640	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
1641	! & wrf_l_wake, wrf_grid%itimestep, &
1642	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1643	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1644	! & UtendValues2, VtendValues2, TtendValues2, &
1645	! & MixingRatiotendValues2, psfctendValues2,wrf_grid%lu10m, wrf_grid%lu10m_ter, &
1646	! & wrf_grid%lu10m_lic, wrf_grid%lu10m_oce, wrf_grid%lu10m_sic, wrf_grid%lue, &
1647	! & wrf_grid%luq, wrf_grid%lustar, wrf_grid%lustar_ter, wrf_grid%lustar_lic, &
1648	! & wrf_grid%lustar_oce, wrf_grid%lustar_sic, wrf_grid%lv10m, &
1649	! & wrf_grid%lv10m_ter, wrf_grid%lv10m_lic, wrf_grid%lv10m_oce, &
1650	! & wrf_grid%lv10m_sic, wrf_grid%lve, wrf_grid%lvq, wrf_grid%lwake_h, &
1651	! & wrf_grid%lwake_s, wrf_grid%lwbeff, wrf_grid%lwbilo_ter, wrf_grid%lwbilo_lic,&
1652	! & wrf_grid%lwbilo_oce, wrf_grid%lwbilo_sic, wrf_grid%lwbils_ter, &
1653	! & wrf_grid%lwbils_lic, wrf_grid%lwbils_oce, wrf_grid%lwbils_sic, &
1654	! & wrf_grid%lweakinv, wrf_grid%lwind10m, wrf_grid%lwind10max, wrf_grid%lzmax_th)
1655
1656	PRINT *, ' Lluis: wrf_grid%la_th', UBOUND(wrf_grid%la_th),' fraca: ',UBOUND(fraca)
1657	! DO iz=wkms, wkme-2
1658	! CALL get_lmdz_out3D_z_i(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, wbdyw, &
1659	! & ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_con, &
1660	! & wrf_l_thermals, wrf_l_wake, &
1661	! & wrf_grid%itimestep, &
1662	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1663	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1664	! & UtendValues2, VtendValues2, TtendValues2, &
1665	! & MixingRatiotendValues2, psfctendValues2, &
1666	! & wrf_grid%la_th(wims:wime,iz,wjms:wjme), &
1667	! & wrf_grid%lbeta_prec(wims:wime,iz,wjms:wjme), &
1668	! & wrf_grid%ldmc(wims:wime,iz,wjms:wjme), &
1669	! & wrf_grid%ldnwd(wims:wime,iz,wjms:wjme), &
1670	! & wrf_grid%ldnwd0(wims:wime,iz,wjms:wjme), &
1671	! & wrf_grid%ldqajs(wims:wime,iz,wjms:wjme), &
1672	! & wrf_grid%ldqcon(wims:wime,iz,wjms:wjme), &
1673	! & wrf_grid%ldqdyn(wims:wime,iz,wjms:wjme), &
1674	! & wrf_grid%ldqeva(wims:wime,iz,wjms:wjme), &
1675	! & wrf_grid%ldqlsc(wims:wime,iz,wjms:wjme), &
1676	! & wrf_grid%ldqlscst(wims:wime,iz,wjms:wjme), &
1677	! & wrf_grid%ldqlscth(wims:wime,iz,wjms:wjme), &
1678	! & wrf_grid%ldqphy(wims:wime,iz,wjms:wjme), &
1679	! & wrf_grid%ldqthe(wims:wime,iz,wjms:wjme), &
1680	! & wrf_grid%ldqvdf(wims:wime,iz,wjms:wjme), &
1681	! & wrf_grid%ldqwak(wims:wime,iz,wjms:wjme))
1682
1683	! CALL get_lmdz_out3D_z_ii(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, &
1684	! & wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, &
1685	! & wrf_l_thermals, wrf_l_wake, &
1686	! & wrf_grid%itimestep, &
1687	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1688	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1689	! & UtendValues2, VtendValues2, TtendValues2, &
1690	! & MixingRatiotendValues2, psfctendValues2, &
1691	! & wrf_grid%ldtajs(wims:wime,iz,wjms:wjme), &
1692	! & wrf_grid%ldtcon(wims:wime,iz,wjms:wjme), &
1693	! & wrf_grid%ldtdis(wims:wime,iz,wjms:wjme), &
1694	! & wrf_grid%ldtdyn(wims:wime,iz,wjms:wjme), &
1695	! & wrf_grid%ldtec(wims:wime,iz,wjms:wjme), &
1696	! & wrf_grid%ldteva(wims:wime,iz,wjms:wjme), &
1697	! & wrf_grid%ld_th(wims:wime,iz,wjms:wjme), &
1698	! & wrf_grid%lcldemi(wims:wime,iz,wjms:wjme), &
1699	! & wrf_grid%lcldtau(wims:wime,iz,wjms:wjme), &
1700	! & wrf_grid%lclwcon(wims:wime,iz,wjms:wjme), &
1701	! & wrf_grid%ldtlif(wims:wime,iz,wjms:wjme), &
1702	! & wrf_grid%ldtlsc(wims:wime,iz,wjms:wjme), &
1703	! & wrf_grid%ldtlschr(wims:wime,iz,wjms:wjme), &
1704	! & wrf_grid%ldtlscst(wims:wime,iz,wjms:wjme), &
1705	! & wrf_grid%ldtlscth(wims:wime,iz,wjms:wjme), &
1706	! & wrf_grid%ldtlw0(wims:wime,iz,wjms:wjme), &
1707	! & wrf_grid%ldtlwr(wims:wime,iz,wjms:wjme), &
1708	! & wrf_grid%ldtoro(wims:wime,iz,wjms:wjme), &
1709	! & wrf_grid%ldtphy(wims:wime,iz,wjms:wjme), &
1710	! & wrf_grid%ldtsw0(wims:wime,iz,wjms:wjme), &
1711	! & wrf_grid%ldtswr(wims:wime,iz,wjms:wjme), &
1712	! & wrf_grid%ldtthe(wims:wime,iz,wjms:wjme), &
1713	! & wrf_grid%ldtvdf(wims:wime,iz,wjms:wjme), &
1714	! & wrf_grid%ldtwak(wims:wime,iz,wjms:wjme))
1715
1716	! CALL get_lmdz_out3D_z_iii(wims, wime, wjms, wjme, iz, ddimx, ddimy, &
1717	! & dimz, wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, &
1718	! & wrf_l_con, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
1719	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1720	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1721	! & UtendValues2, VtendValues2, TtendValues2, &
1722	! & MixingRatiotendValues2, psfctendValues2, &
1723	! & wrf_grid%lducon(wims:wime,iz,wjms:wjme), &
1724	! & wrf_grid%ldudyn(wims:wime,iz,wjms:wjme), &
1725	! & wrf_grid%ldulif(wims:wime,iz,wjms:wjme), &
1726	! & wrf_grid%lduoro(wims:wime,iz,wjms:wjme), &
1727	! & wrf_grid%lduvdf(wims:wime,iz,wjms:wjme), &
1728	! & wrf_grid%ldvcon(wims:wime,iz,wjms:wjme), &
1729	! & wrf_grid%ldvdyn(wims:wime,iz,wjms:wjme), &
1730	! & wrf_grid%ldvlif(wims:wime,iz,wjms:wjme), &
1731	! & wrf_grid%ldvoro(wims:wime,iz,wjms:wjme), &
1732	! & wrf_grid%ldvvdf(wims:wime,iz,wjms:wjme), &
1733	! & wrf_grid%le_th(wims:wime,iz,wjms:wjme), &
1734	! & wrf_grid%levu(wims:wime,iz,wjms:wjme), &
1735	! & wrf_grid%lfl(wims:wime,iz,wjms:wjme), &
1736	! & wrf_grid%lgeop(wims:wime,iz,wjms:wjme), &
1737	! & wrf_grid%lh2o(wims:wime,iz,wjms:wjme), &
1738	! & wrf_grid%liwcon(wims:wime,iz,wjms:wjme), &
1739	! & wrf_grid%lkz(wims:wime,iz,wjms:wjme), &
1740	! & wrf_grid%lkz_max(wims:wime,iz,wjms:wjme), &
1741	! & wrf_grid%llambda_th(wims:wime,iz,wjms:wjme), &
1742	! & wrf_grid%llwcon(wims:wime,iz,wjms:wjme), &
1743	! & wrf_grid%lma(wims:wime,iz,wjms:wjme), &
1744	! & wrf_grid%lmc(wims:wime,iz,wjms:wjme), &
1745	! & wrf_grid%lmcd(wims:wime,iz,wjms:wjme), &
1746	! & wrf_grid%lmass(wims:wime,iz,wjms:wjme))
1747
1748	! CALL get_lmdz_out3D_z_iv(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, &
1749	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
1750	! & wrf_l_wake, wrf_grid%itimestep, &
1751	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1752	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1753	! & UtendValues2, VtendValues2, TtendValues2, &
1754	! & MixingRatiotendValues2, psfctendValues2, &
1755	! & wrf_grid%loliq(wims:wime,iz,wjms:wjme), &
1756	! & wrf_grid%lovap(wims:wime,iz,wjms:wjme), &
1757	! & wrf_grid%lovapinit(wims:wime,iz,wjms:wjme), &
1758	! & wrf_grid%lpaprs(wims:wime,iz,wjms:wjme), &
1759	! & wrf_grid%lpr_con_i(wims:wime,iz,wjms:wjme), &
1760	! & wrf_grid%lpr_con_l(wims:wime,iz,wjms:wjme), &
1761	! & wrf_grid%lpres(wims:wime,iz,wjms:wjme), &
1762	! & wrf_grid%lpr_lsc_i(wims:wime,iz,wjms:wjme), &
1763	! & wrf_grid%lpr_lsc_l(wims:wime,iz,wjms:wjme), &
1764	! & wrf_grid%lptconv(wims:wime,iz,wjms:wjme), &
1765	! & wrf_grid%lq_th(wims:wime,iz,wjms:wjme), &
1766	! & wrf_grid%lratqs(wims:wime,iz,wjms:wjme), &
1767	! & wrf_grid%lre(wims:wime,iz,wjms:wjme), &
1768	! & wrf_grid%lref_ice(wims:wime,iz,wjms:wjme), &
1769	! & wrf_grid%lref_liq(wims:wime,iz,wjms:wjme), &
1770	! & wrf_grid%lrhum(wims:wime,iz,wjms:wjme), &
1771	! & wrf_grid%lrld(wims:wime,iz,wjms:wjme), &
1772	! & wrf_grid%lrldcs(wims:wime,iz,wjms:wjme), &
1773	! & wrf_grid%lrlu(wims:wime,iz,wjms:wjme), &
1774	! & wrf_grid%lrlucs(wims:wime,iz,wjms:wjme), &
1775	! & wrf_grid%lrneb(wims:wime,iz,wjms:wjme), &
1776	! & wrf_grid%lrnebcon(wims:wime,iz,wjms:wjme), &
1777	! & wrf_grid%lrnebls(wims:wime,iz,wjms:wjme), &
1778	! & wrf_grid%lrsd(wims:wime,iz,wjms:wjme), &
1779	! & wrf_grid%lrsdcs(wims:wime,iz,wjms:wjme), &
1780	! & wrf_grid%lrsu(wims:wime,iz,wjms:wjme), &
1781	! & wrf_grid%lrsucs(wims:wime,iz,wjms:wjme))
1782
1783	! CALL get_lmdz_out3D_z_v(wims, wime, wjms, wjme, iz, ddimx, ddimy, dimz, &
1784	! & wbdyw, ddimxy, ip, jp, lp, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, 1, &
1785	! & wrf_l_con, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
1786	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1787	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1788	! & UtendValues2, VtendValues2, TtendValues2, &
1789	! & MixingRatiotendValues2, psfctendValues2, &
1790	! & wrf_grid%lsens_ter(wims:wime,iz,wjms:wjme), &
1791	! & wrf_grid%lsens_lic(wims:wime,iz,wjms:wjme), &
1792	! & wrf_grid%lsens_oce(wims:wime,iz,wjms:wjme), &
1793	! & wrf_grid%lsens_sic(wims:wime,iz,wjms:wjme), &
1794	! & wrf_grid%ltemp(wims:wime,iz,wjms:wjme), &
1795	! & wrf_grid%ltheta(wims:wime,iz,wjms:wjme), &
1796	! & wrf_grid%ltke(wims:wime,iz,wjms:wjme), &
1797	! & wrf_grid%ltke_max(wims:wime,iz,wjms:wjme), &
1798	! & wrf_grid%ltke_ter(wims:wime,iz,wjms:wjme), &
1799	! & wrf_grid%ltke_lic(wims:wime,iz,wjms:wjme), &
1800	! & wrf_grid%ltke_oce(wims:wime,iz,wjms:wjme), &
1801	! & wrf_grid%ltke_sic(wims:wime,iz,wjms:wjme), &
1802	! & wrf_grid%ltnhus(wims:wime,iz,wjms:wjme), &
1803	! & wrf_grid%ltnhusc(wims:wime,iz,wjms:wjme), &
1804	! & wrf_grid%ltnhusscpbl(wims:wime,iz,wjms:wjme), &
1805	! & wrf_grid%ltnt(wims:wime,iz,wjms:wjme), &
1806	! & wrf_grid%ltntc(wims:wime,iz,wjms:wjme), &
1807	! & wrf_grid%ltntr(wims:wime,iz,wjms:wjme), &
1808	! & wrf_grid%ltntscpbl(wims:wime,iz,wjms:wjme), &
1809	! & wrf_grid%lupwd(wims:wime,iz,wjms:wjme), &
1810	! & wrf_grid%lvitu(wims:wime,iz,wjms:wjme), &
1811	! & wrf_grid%lvitv(wims:wime,iz,wjms:wjme), &
1812	! & wrf_grid%lvitw(wims:wime,iz,wjms:wjme), &
1813	! & wrf_grid%lvprecip(wims:wime,iz,wjms:wjme), &
1814	! & wrf_grid%lwake_omg(wims:wime,iz,wjms:wjme), &
1815	! & wrf_grid%lwdtraina(wims:wime,iz,wjms:wjme), &
1816	! & wrf_grid%lwdtrainm(wims:wime,iz,wjms:wjme), &
1817	! & wrf_grid%lzfull(wims:wime,iz,wjms:wjme), &
1818	! & wrf_grid%lzhalf(wims:wime,iz,wjms:wjme))
1819
1820	! END DO
1821
1822	PRINT *,' Lluis Loading done! iz ldtlwr _______________'
1823	DO iz=1,dimz
1824	PRINT *,iz,wrf_grid%ldtlwr(ip,iz,jp)
1825	END DO
1826
1827	! CALL get_lmdz_out3D_ii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1828	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, &
1829	! & wrf_l_wake, wrf_grid%itimestep, &
1830	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1831	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1832	! & UtendValues2, VtendValues2, TtendValues2, &
1833	! & MixingRatiotendValues2, psfctendValues2, &
1834	! & wrf_grid%ldtajs, &
1835	! & wrf_grid%ldtcon, wrf_grid%ldtdis, wrf_grid%ldtdyn, wrf_grid%ldtec, &
1836	! & wrf_grid%ldteva, wrf_grid%ld_th, wrf_grid%lcldemi, wrf_grid%lcldtau, &
1837	! & wrf_grid%lclwcon, wrf_grid%ldtlif, wrf_grid%ldtlsc, wrf_grid%ldtlschr, &
1838	! & wrf_grid%ldtlscst, wrf_grid%ldtlscth, wrf_grid%ldtlw0, wrf_grid%ldtlwr, &
1839	! & wrf_grid%ldtoro, wrf_grid%ldtphy, wrf_grid%ldtsw0, wrf_grid%ldtswr, &
1840	! & wrf_grid%ldtthe, wrf_grid%ldtvdf, wrf_grid%ldtwak)
1841
1842	! CALL get_lmdz_out3D_iii(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1843	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
1844	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1845	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1846	! & UtendValues2, VtendValues2, TtendValues2, &
1847	! & MixingRatiotendValues2, psfctendValues2, &
1848	! & wrf_grid%lducon, &
1849	! & wrf_grid%ldudyn, wrf_grid%ldulif, wrf_grid%lduoro, wrf_grid%lduvdf, &
1850	! & wrf_grid%ldvcon, wrf_grid%ldvdyn, wrf_grid%ldvlif, wrf_grid%ldvoro, &
1851	! & wrf_grid%ldvvdf, wrf_grid%le_th, wrf_grid%levu, wrf_grid%lfl,wrf_grid%lgeop,&
1852	! & wrf_grid%lh2o, wrf_grid%liwcon, wrf_grid%lkz, wrf_grid%lkz_max, &
1853	! & wrf_grid%llambda_th, wrf_grid%llwcon, wrf_grid%lma, wrf_grid%lmass)
1854
1855	! CALL get_lmdz_out3D_iv(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1856	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
1857	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1858	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1859	! & UtendValues2, VtendValues2, TtendValues2, &
1860	! & MixingRatiotendValues2, psfctendValues2, &
1861	! & wrf_grid%loliq, wrf_grid%lovap, &
1862	! & wrf_grid%lovapinit, wrf_grid%lpaprs, wrf_grid%lpr_con_i, wrf_grid%lpr_con_l,&
1863	! & wrf_grid%lpres, wrf_grid%lpr_lsc_i, wrf_grid%lpr_lsc_l, wrf_grid%lptconv, &
1864	! & wrf_grid%lq_th, wrf_grid%lratqs, wrf_grid%lre, wrf_grid%lref_ice, &
1865	! & wrf_grid%lref_liq, wrf_grid%lrhum, wrf_grid%lrld, wrf_grid%lrldcs, &
1866	! & wrf_grid%lrlu, wrf_grid%lrlucs, wrf_grid%lrneb, wrf_grid%lrnebcon, &
1867	! & wrf_grid%lrnebls, wrf_grid%lrsd, wrf_grid%lrsdcs, wrf_grid%lrsu, &
1868	! & wrf_grid%lrsucs)
1869
1870	! CALL get_lmdz_out3D_v(wims, wime, wjms, wjme, wkms, wkme, ddimx, ddimy, dimz, &
1871	! & wbdyw, ddimxy, nbsrf, nsoilmx, lmdzmixingratios, 1, 2, wrf_l_thermals, wrf_l_wake, wrf_grid%itimestep, &
1872	! & wtime_step, FullPressValues2, HalfPressValues2, GeopotValues2, &
1873	! & sfcGeopotValues, lmdzoutP, MixingRatioValues2, &
1874	! & UtendValues2, VtendValues2, TtendValues2, &
1875	! & MixingRatiotendValues2, psfctendValues2, &
1876	! & wrf_grid%lsens_ter, wrf_grid%lsens_lic, wrf_grid%lsens_oce, &
1877	! & wrf_grid%lsens_sic, wrf_grid%ltemp, wrf_grid%ltheta, wrf_grid%ltke, &
1878	! & wrf_grid%ltke_max, wrf_grid%ltke_ter, wrf_grid%ltke_lic, wrf_grid%ltke_oce, &
1879	! & wrf_grid%ltke_sic, wrf_grid%ltnhus, wrf_grid%ltnhusc, wrf_grid%ltnhusscpbl, &
1880	! & wrf_grid%ltnt, wrf_grid%ltntc, wrf_grid%ltntr, wrf_grid%ltntscpbl, &
1881	! & wrf_grid%lupwd, wrf_grid%lvitu, wrf_grid%lvitv, wrf_grid%lvitw, &
1882	! & wrf_grid%lwake_deltaq, wrf_grid%lwake_deltat, wrf_grid%lwake_omg, &
1883	! & wrf_grid%lwdtraina, wrf_grid%lwdtrainm, wrf_grid%lzfull, wrf_grid%lzhalf)
1884
1885	DEALLOCATE(MixingRatioValues)
1886	DEALLOCATE(MixingRatiotendValues)
1887
1888	WRITE(message, *)' exiting: module_lmdz_phys.F'
1889	CALL wrf_debug(75, message)
1890
1891	END SUBROUTINE call_lmdz_phys
1892
1893	END MODULE module_lmdz_phys

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

Download in other formats:

Original Format