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

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

No Ale, Alp, ... values because pphi is at the intermediate level!!!
Removing raw prints of checkig for Ale, Alp

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