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

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	REAL, 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 (0, 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