source: lmdz_wrf/trunk/WRFV3/phys/module_lmdz_phys.F @ 1393

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

Removiong psfcTendValues printing

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