source: lmdz_wrf/WRFV3/lmdz/wrf_lmdz_mod.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 119.6 KB
Line 
1MODULE wrf_lmdz_mod
2!
3! Module for the interface 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!   November 2013
7!
8
9! WRF modules
10!    USE module_domain_type
11
12!!!!!!! Subroutines
13! wrf_varKsoil: Subroutine to provide WRF variable values from the 4 LMDZ soil types
14! lmdz_varKsoil: Subroutine to provide variable values to the 4 LMDZ soil types according to different methods
15! LMDZmoist_WRFmoist: Subroutine to transform from LMDZ moisture array to the WRF moisture array
16! WRFmoist_LMDZmoist: Subroutine to transform from WRF moisture array to the LMDZ moisture array
17! eta_to_pressure: Subroutine to transform from eta levels to pressure levels
18! temp_to_theta: Subroutine to transform from temperatures to WRF potential temperatures
19! theta_to_temp: Subroutine to transform from WRF potential temperatures to temperature
20! temp_to_theta1: Subroutine to transform from temperature to potential temperature
21! theta_to_temp1: Subroutine to transform from potential temperature to temperature
22! def_get_scalar_value: Subroutine to get a scalar value from a .def file
23!   NOTE: For that scalar values with scientific notation (EN, or e-N) it will not work...
24!     next version!
25! value_kind: Subroutine to determine which kind of value is given
26! string_split: Subtoutine to split a string according to a character and returns the 'Nsec' of the split
27! load_lmdz: Subroutine to load LMDZ variables with values from WRF
28! load_lmdz_limit: Subroutine to load LMDZ limit variables with values from WRF
29! get_lmdz: Subroutine to get LMDZ values and fill WRF variables
30! get_lowbdy: Subroutine to load LMDZ limit variables with values from lowbdy WRF
31! vect_mat: Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
32! vect_mat_zstagg: Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
33! mat_vect: Subroutine to transform a 3D matrix to a LMDZ physic 1D vector
34! mat_vect_zstagg: Subroutine to transform a 3D matrix z-staggered to a LMDZ physic 1D vector
35! WRFlanduse_LMDZKsoil: Subroutine to retrieve LMDZ Ksoil classification using WRF landuse
36! interpolate_layers: Subroutine to interpolate values between different layers
37!   (constant values along the depth of a layer)
38! interpolate_lin: Program to interpolate values y=a+b*x with: (from clWRF modifications)
39! table_characteristics: Subroutine to read values from a WRF .TBL
40! read_table: Subroutine to read values from a WRF .TBL
41! compute_soil_mass: Subroutine to compute the total soil mass (kg) from a point with
42!   different percentages of soil types
43! compute_TOTsoil_mass_water: Subroutine to compute total soil mass and water from a
44!   given grid point using a given data-set from SOILPARM.TBL
45! compute_TOTsoil_mass_water_values: Subroutine to compute total soil mass and water
46!   from a given grid point using the values from a given data-set of SOILPARM.TBL
47! get_lmdz_out: Subroutine to get all the LMDZ output variables into the WRF
48!    Registry structure
49! put_lmdz_in_WRFout: Subroutine to get LMDZ output and get into WRF standard output variables
50
51! LMDZ modules
52!    USE indice_sol_mod
53
54!    IMPLICIT NONE
55
56!    INCLUDE "dimsoil.h"
57
58  CONTAINS
59
60  SUBROUTINE wrf_varKsoil(is,ie,js,je,dx,dy,wbdy,kmethod,ter,lic,oce,sic,vter,vlic,  &
61    voce,vsic,var)
62! Subroutine to provide WRF variable values from the 4 LMDZ soil types
63
64    IMPLICIT NONE
65
66    INTEGER, INTENT(IN)                                  :: is,ie,js,je,dx,dy,wbdy
67    CHARACTER(LEN=50), INTENT(IN)                        :: kmethod
68    REAL, DIMENSION(is:ie,js:je), INTENT(IN)             :: ter,lic,oce,sic
69    REAL, DIMENSION(is:ie,js:je), INTENT(IN)             :: vter,vlic,voce,vsic
70    REAL, DIMENSION(is:ie,js:je), INTENT(OUT)            :: var
71
72! Local
73    INTEGER                                              :: i,j,ddx,ddy
74    CHARACTER(LEN=50)                                    :: errmsg, fname
75    REAL*8, DIMENSION(is:ie,js:je)                       :: prop
76
77!!!!!!! Variables
78! d[x/y]: dimension of the matrices
79! kmethod: method to use
80!   'prod': product, variable is distributed by percentage of the soil type vsoil=var*soil
81! ter: percentage of continent
82! lic: percentage of frozen continent
83! oce: percentage of ocean
84! sic: percentage of frozen ocean
85! vter: variable on continent
86! vlic: variable on frozen continent
87! voce: variable on ocean
88! vsic: variable on frozen ocean
89! prop: proportion between soil fractions
90! var: variable to be reinterpreted from LMDZ soil types
91
92    fname="'wrf_varksoil'"
93    errmsg='ERROR -- error -- ERROR -- error'
94
95    prop=-1.d0
96    var=0.
97    ddx=dx-2*wbdy
98    ddy=dy-2*wbdy
99
100    methodkind: IF (TRIM(kmethod) == 'prod') THEN
101      DO i=1,ddx
102        DO j=1,ddy
103          var(i,j)=vter(i,j)*ter(i,j)+vlic(i,j)*lic(i,j)+voce(i,j)*oce(i,j)+vsic(i,j)*sic(i,j)
104        END DO
105      END DO
106
107    ELSE
108      PRINT *,TRIM(errmsg)
109      PRINT *,'   ' // TRIM(fname) // ": method '" // TRIM(kmethod) // "' does not exist !!!!!"
110      STOP
111
112    END IF methodkind
113
114  END SUBROUTINE wrf_varKsoil
115
116  SUBROUTINE lmdz_varKsoil(is,ie,js,je,dx,dy,wbdy,var,kmethod,ter,lic,oce,sic,vter,  &
117    vlic,voce,vsic)
118! Subroutine to provide variable values to the 4 LMDZ soil types according to
119!   different methods
120
121    IMPLICIT NONE
122
123    INTEGER, INTENT(IN)                                  :: is,ie,js,je,dx,dy,wbdy
124    REAL, DIMENSION(is:ie,js:je), INTENT(IN)             :: var
125    CHARACTER(LEN=50), INTENT(IN)                        :: kmethod
126    REAL, DIMENSION(is:ie,js:je), INTENT(IN)             :: ter,lic,oce,sic
127    REAL, DIMENSION(is:ie,js:je), INTENT(OUT)            :: vter,vlic,voce,vsic
128
129! Local
130    INTEGER                                              :: i,j,ddx,ddy
131    CHARACTER(LEN=50)                                    :: errmsg, fname
132    REAL*8, DIMENSION(is:ie,js:je)                       :: prop
133
134!!!!!!! Variables
135! d[x/y]: dimension of the matrices
136! var: variable to be reinterpreted in LMDZ soil types
137! kmethod: method to use
138!   'NOlic': NO lic, variable is not located in ice continent soil type: lic
139!   'NOlnd': NO land, variable is distributed only to sea soil types: oce, sic
140!   'NOneg': NO negative, variable is distributed only to free ice soil types: ter, oce
141!   'NOoce': NO oce, variable is not located in ocean soil type: oce
142!   'NOpos': NO positive, variable is distributed only to iced soil types: lic, sic
143!   'NOsea': NO sea, variable is distributed only to land soil types: ter, lic
144!   'NOsic': NO sic, variable is not located in ice ocean soil type: sic
145!   'NOter': NO ter, variable is not located in land soil type: ter
146!   'prod': product, variable is distributed by percentage of the soil type vsoil=var*soil
147!   'lic:' lic, variable only is located in ice land soil type: lic
148!   'oce:' oce, variable only is located in ocean soil type: oce
149!   'sic:' sic, variable only is located in ice ocean soil type: sic
150!   'ter:' ter, variable only is located in land soil type: ter
151! ter: percentage of continent
152! lic: percentage of frozen continent
153! oce: percentage of ocean
154! sic: percentage of frozen ocean
155! vter: variable on continent
156! vlic: variable on frozen continent
157! voce: variable on ocean
158! vsic: variable on frozen ocean
159! prop: proportion between soil fractions
160
161    fname="'lmdz_varksoil'"
162    errmsg='ERROR -- error -- ERROR -- error'
163
164    ddx=dx-2*wbdy
165    ddy=dy-2*wbdy
166
167    prop=-1.d0
168    vter=0.
169    vlic=0.
170    voce=0.
171    vsic=0.
172
173    methodkind: IF (TRIM(kmethod) == 'NOlic') THEN
174      WHERE (ter+oce+sic /= 0.) prop=1.d0/(ter+oce+sic)
175      DO i=1,ddx
176        DO j=1,ddy
177          IF (prop(i,j) >= 0.) THEN
178            vter(i,j)=var(i,j)*ter(i,j)*prop(i,j)
179            voce(i,j)=var(i,j)*oce(i,j)*prop(i,j)
180            vsic(i,j)=var(i,j)*sic(i,j)*prop(i,j)
181          END IF
182        END DO
183      END DO
184
185    ELSEIF (TRIM(kmethod) == 'NOlnd') THEN
186      WHERE (oce+sic /= 0.) prop=oce*1.d0/(oce+sic)
187      DO i=1,ddx
188        DO j=1,ddy
189          IF (prop(i,j) >= 0.) THEN
190            voce(i,j)=var(i,j)*prop(i,j)
191            vsic(i,j)=var(i,j)*(1.-prop(i,j))
192          END IF
193        END DO
194      END DO
195
196    ELSEIF (TRIM(kmethod) == 'NOneg') THEN
197      WHERE (ter+oce /= 0.) prop=ter*1.d0/(ter+oce)
198      DO i=1,ddx
199        DO j=1,ddy
200          IF (prop(i,j) >= 0.) THEN
201            vter(i,j)=var(i,j)*prop(i,j)
202            voce(i,j)=var(i,j)*(1.-prop(i,j))
203          END IF
204        END DO
205      END DO
206
207    ELSEIF (TRIM(kmethod) == 'NOoce') THEN
208      WHERE (ter+lic+sic /= 0.) prop=1.d0/(ter+lic+sic)
209      DO i=1,ddx
210        DO j=1,ddy
211          IF (prop(i,j) >= 0.) THEN
212            vter(i,j)=var(i,j)*ter(i,j)*prop(i,j)
213            vlic(i,j)=var(i,j)*lic(i,j)*prop(i,j)
214            vsic(i,j)=var(i,j)*sic(i,j)*prop(i,j)
215          END IF
216        END DO
217      END DO
218
219    ELSEIF (TRIM(kmethod) == 'NOpos') THEN
220      WHERE (lic+sic /=0.) prop=lic*1.d0/(lic+sic)
221      DO i=1,ddx
222        DO j=1,ddy
223          IF (prop(i,j) >= 0.) THEN
224            vlic(i,j)=var(i,j)*prop(i,j)
225            vsic(i,j)=var(i,j)*(1.-prop(i,j))
226          END IF
227        END DO
228      END DO
229
230    ELSEIF (TRIM(kmethod) == 'NOsea') THEN
231      WHERE (ter+lic /= 0.) prop=ter*1.d0/(ter+lic)
232      DO i=1,ddx
233        DO j=1,ddy
234          IF (prop(i,j) >= 0.) THEN
235            vter(i,j)=var(i,j)*prop(i,j)
236            vlic(i,j)=var(i,j)*(1.-prop(i,j))
237          END IF
238        END DO
239      END DO
240
241    ELSEIF (TRIM(kmethod) == 'NOsic') THEN
242      WHERE (ter+lic+oce /= 0.) prop=1.d0/(ter+lic+oce)
243      DO i=1,ddx
244        DO j=1,ddy
245          IF (prop(i,j) >= 0.) THEN
246            vter(i,j)=var(i,j)*ter(i,j)*prop(i,j)
247            vlic(i,j)=var(i,j)*lic(i,j)*prop(i,j)
248            voce(i,j)=var(i,j)*oce(i,j)*prop(i,j)
249          END IF
250        END DO
251      END DO
252
253    ELSEIF (TRIM(kmethod) == 'NOter') THEN
254      WHERE (lic+sic+oce /= 0.) prop=1.d0/(lic+sic+oce)
255      DO i=1,ddx
256        DO j=1,ddy
257          IF (prop(i,j) >= 0.) THEN
258            vlic(i,j)=var(i,j)*lic(i,j)*prop(i,j)
259            vsic(i,j)=var(i,j)*sic(i,j)*prop(i,j)
260            voce(i,j)=var(i,j)*oce(i,j)*prop(i,j)
261          END IF
262        END DO
263      END DO
264
265    ELSEIF (TRIM(kmethod) == 'lic') THEN
266      DO i=1,ddx
267        DO j=1,ddy
268          IF (lic(i,j) /= 0.) THEN
269            vlic(i,j)=var(i,j)
270          END IF
271        END DO
272      END DO
273
274    ELSEIF (TRIM(kmethod) == 'oce') THEN
275      DO i=1,ddx
276        DO j=1,ddy
277          IF (oce(i,j) /= 0.) THEN
278            voce(i,j)=var(i,j)
279          END IF
280        END DO
281      END DO
282
283    ELSEIF (TRIM(kmethod) == 'prod') THEN
284      DO i=1,ddx
285        DO j=1,ddy
286          vter(i,j)=var(i,j)*ter(i,j)*1.d0
287          vlic(i,j)=var(i,j)*lic(i,j)*1.d0
288          voce(i,j)=var(i,j)*oce(i,j)*1.d0
289          vsic(i,j)=var(i,j)*sic(i,j)*1.d0
290        END DO
291      END DO
292
293    ELSEIF (TRIM(kmethod) == 'sic') THEN
294      DO i=1,ddx
295        DO j=1,ddy
296          IF (sic(i,j) /= 0.) THEN
297            vsic(i,j)=var(i,j)
298          END IF
299        END DO
300      END DO
301
302    ELSEIF (TRIM(kmethod) == 'ter') THEN
303      DO i=1,ddx
304        DO j=1,ddy
305          IF (ter(i,j) /= 0.) THEN
306            vter(i,j)=var(i,j)
307          END IF
308        END DO
309      END DO
310
311    ELSE
312      PRINT *,TRIM(errmsg)
313      PRINT *,'   ' // TRIM(fname) // ": method '" // TRIM(kmethod) // "' does not exist !!!!!"
314      STOP
315    END IF methodkind
316
317  END SUBROUTINE lmdz_varKsoil
318
319
320  SUBROUTINE LMDZmoist_WRFmoist(lmdzMixRat, dx, dy, dz, Nwmoist, wbdy, wqvid, wqcid, &
321    & wqrid, Nlmdzq, wMOIST)
322! Subroutine to transform from LMDZ moisture array to the WRF moisture array
323
324    IMPLICIT NONE
325
326    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
327    INTEGER, INTENT(IN)                                  :: Nwmoist, wqvid, wqcid,   &
328      wqrid
329    INTEGER, INTENT(IN)                                  :: Nlmdzq
330    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz,Nlmdzq), INTENT(IN) ::                &
331      lmdzMixRat
332    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1,Nwmoist), INTENT(OUT) :: wMOIST
333
334! Local
335    INTEGER                                              :: i,j,k, iq
336    INTEGER                                              :: ddx, ddy
337    INTEGER, ALLOCATABLE, DIMENSION(:)                   :: widlmdz
338    REAL, DIMENSION(1-wbdy:dx-wbdy+1,1-wbdy:dy-wbdy+1)   :: qcr_ratio, wcrMOIST
339
340!!!!!!! Variables
341! widlmdz: vector with the id of the WRF tracers to be also included in the LMDZ moist array
342! d[x,y,z]: Spatial dimension of the WRF moisture array
343! Nwmoist: number of species in wMOIST
344! wqvid: index for the water vapor mixing ratio
345! wqcid: index for the cloud water mixing ratio
346! wqrid: index for the rain water mixing ratio
347! Nlmdzq: number of species in LMDZ
348! lmdxMixRat: LMDZ mixing ratio array
349! qcr_ratio: ratio between cloud/rain liquid water in WRF moisture array
350! wcrMOIST: WRF array combination of cloud and rain liquid water
351! wMOIST: WRF moisture array
352
353!!!!!! L. Fita, July 203
354! At this stage, we only care about water vapor, cloud liquid water and rain liquid water.
355
356    ddx=dx-2*wbdy
357    ddy=dy-2*wbdy
358
359    wMOIST=0.
360    wcrMOIST=0.
361
362    DO k=1,dz
363      DO j=1,ddy
364        DO i=1,ddx
365          wMOIST(i,k,j,wqvid) = lmdzMixRat(ddx*(j-1)+i,k,1)
366          wMOIST(i,k,j,wqcid) = lmdzMixRat(ddx*(j-1)+i,k,2)
367!!          wcrMOIST(i,j) = lmdzMixRat(ddx*(j-1)+1,k,2)
368        END DO
369      END DO
370!! L. Fita, July 2013
371!! We can not know how to split liquid water between cloud and rain. Thus, we use the previous
372!!   time-step proportion
373!!      qcr_ratio=wMOIST(:,k,:,wqcid)/(wMOIST(:,k,:,wqcid) + wMOIST(:,k,:,wqrid))
374!!      wMOIST(:,k,:,wqcid) = wcrMOIST*qcr_ratio
375!!      wMOIST(:,k,:,wqrid) = wcrMOIST*(1.-qcr_ratio)
376    ENDDO
377
378  END SUBROUTINE LMDZmoist_WRFmoist
379
380
381  SUBROUTINE WRFmoist_LMDZmoist(wMOIST, dx, dy, dz, Nwmoist, wbdy, wqvid, wqcid,     &
382    & wqrid, Nlmdzq, lmdzMixRat)
383! Subroutine to transform from WRF moisture array to the LMDZ moisture array
384
385    IMPLICIT NONE
386
387    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
388    INTEGER, INTENT(IN)                                  :: Nwmoist, wqvid, wqcid, wqrid
389    INTEGER, INTENT(IN)                                  :: Nlmdzq
390    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1,Nwmoist), INTENT(IN) :: wMOIST
391    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz,Nlmdzq), INTENT(OUT) :: lmdzMixRat
392
393! Local
394    INTEGER                                              :: i, j, k, iq, dx2, dy2
395    INTEGER                                              :: ddx, ddy
396!    INTEGER, ALLOCATABLE, DIMENSION(:)                   :: widlmdz
397!    REAL, DIMENSION(dx,dy)                               :: wcrMOIST
398    REAL, DIMENSION(1-wbdy:dx-wbdy+1,1-wbdy:dy-wbdy+1)   :: qcr_ratio, wcrMOIST
399
400
401!!!!!!! Variables
402! wMOIST: WRF moisture array
403! d[x,y,z]: Spatial dimension of the WRF moisture array
404! Nwmoist: number of species in wMOIST
405! wqvid: index for the water vapor mixing ratio
406! wqcid: index for the cloud water mixing ratio
407! wqrid: index for the rain water mixing ratio
408! wcrMOIST: WRF array combination of cloud and rain liquid water
409! Nlmdzq: number of species in LMDZ
410! lmdxMixRat: LMDZ mixing ratio array
411! widlmdz: vector with the id of the WRF tracers to be also included in the LMDZ moist array
412
413    ddx=dx-2*wbdy
414    ddy=dy-2*wbdy
415
416    lmdzMixRat=0.
417
418    DO k=1,dz
419! Only taking water vapor and cloud water mixing ratios
420      DO j=1,ddy
421        DO i=1,ddx
422          lmdzMixRat(ddx*(j-1)+i,k,1) = wMOIST(i,k,j,wqvid)
423          lmdzMixRat(ddx*(j-1)+i,k,2) = wMOIST(i,k,j,wqcid)
424        END DO
425      END DO
426    ENDDO
427
428!!    IF (Nlmdzq > 2) THEN
429! Whenever we would like to use more WRF water spicies some sort of mapping should be needed.
430!   Maybe  something like this:
431!      IF (ALLOCATED(widlmdz)) DEALLOCATE(widlmdz)
432!      ALLOCATE(widlmdz(Nlmdzq-2))
433
434!!      DO iq=3,Nlmdzq
435!!        DO k=1,dz
436!!          lmdzMixRat(:,k,iq) = (iq-3.+(dz-k-1.)*1./(dz+1.))*1.
437!!        END DO
438!!      END DO
439!!    END IF
440
441!    DEALLOCATE(widlmdz)
442
443  END SUBROUTINE WRFmoist_LMDZmoist
444
445
446  SUBROUTINE eta_to_pressure(etaP, psfc, ptop, dz, presP)
447! Subroutine to transform from eta levels to pressure levels
448
449    IMPLICIT NONE
450
451    INTEGER, INTENT(IN)                                  :: dz
452    REAL, INTENT(IN)                                     :: psfc, ptop
453    REAL, DIMENSION(dz), INTENT(IN)                      :: etaP
454    REAL, DIMENSION(dz), INTENT(OUT)                     :: presP
455
456!!!!!!Variables
457! etaP: vector with the eta values
458! psfc: pressure at the surface (in Pa)
459! ptop: pressure at the top (in Pa)
460! dz: number of vertical levels
461! presP: equivalent vector of pressures (in Pa)
462
463    presP = etaP*(psfc - ptop) + ptop
464
465  END SUBROUTINE eta_to_pressure
466
467
468  SUBROUTINE temp_to_theta(tempValues, pressValues, thetaValues,                     &
469    &                  ids,ide, jds,jde, kds,kde,                                    &
470    &                  ims,ime, jms,jme, kms,kme,                                    &
471    &                  ips,ipe, jps,jpe, kps,kpe,                                    & ! patch  dims
472    &                  i_start,i_end,j_start,j_end,kts,kte,num_tiles                 &
473   )
474! Subroutine to transform from temperatures to WRF potential temperatures
475
476    USE module_model_constants
477
478    IMPLICIT NONE
479
480!-- thetaValues   potential temperature values from WRF (300. to be added)
481!-- pressValues   pressure
482!-- tempValues    temperature values in K
483!
484!-- ids           start index for i in domain
485!-- ide           end index for i in domain
486!-- jds           start index for j in domain
487!-- jde           end index for j in domain
488!-- kds           start index for k in domain
489!-- kde           end index for k in domain
490!-- ims           start index for i in memory
491!-- ime           end index for i in memory
492!-- jms           start index for j in memory
493!-- jme           end index for j in memory
494!-- ips           start index for i in patch
495!-- ipe           end index for i in patch
496!-- jps           start index for j in patch
497!-- jpe           end index for j in patch
498!-- kms           start index for k in memory
499!-- kme           end index for k in memory
500!-- i_start       start indices for i in tile
501!-- i_end         end indices for i in tile
502!-- j_start       start indices for j in tile
503!-- j_end         end indices for j in tile
504!-- kts           start index for k in tile
505!-- kte           end index for k in tile
506!-- num_tiles     number of tiles
507!
508!======================================================================
509
510    INTEGER,      INTENT(IN)                             ::                          &
511                                      ids,ide, jds,jde, kds,kde,                     &
512                                      ims,ime, jms,jme, kms,kme,                     &
513                                      ips,ipe, jps,jpe, kps,kpe,                     &
514                                                        kts,kte,                     &
515                                                      num_tiles
516
517    INTEGER, DIMENSION(num_tiles), INTENT(IN)            ::                          &
518       &           i_start,i_end,j_start,j_end
519
520
521    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: tempValues,            &
522       &           pressValues
523    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: thetaValues
524
525! Local
526    REAL, PARAMETER                                      :: pres0 = 100000.
527
528    thetaValues = tempValues * ( pres0 / pressValues ) ** (r_d / cp)
529
530    thetaValues = thetaValues - 300.
531
532  END SUBROUTINE temp_to_theta
533
534  SUBROUTINE theta_to_temp(thetaValues, pressValues, tempValues,                     &
535    &                  ids,ide, jds,jde, kds,kde,                                    &
536    &                  ims,ime, jms,jme, kms,kme,                                    &
537    &                  ips,ipe, jps,jpe, kps,kpe,                                    & ! patch  dims
538    &                  i_start,i_end,j_start,j_end,kts,kte,num_tiles                 &
539   )
540! Subroutine to transform from WRF potential temperatures to temperature
541
542    USE module_model_constants
543
544    IMPLICIT NONE
545
546!-- thetaValues   potential temperature values from WRF (300. to be added)
547!-- pressValues   pressure
548!-- tempValues    temperature values in K
549!
550!-- ids           start index for i in domain
551!-- ide           end index for i in domain
552!-- jds           start index for j in domain
553!-- jde           end index for j in domain
554!-- kds           start index for k in domain
555!-- kde           end index for k in domain
556!-- ims           start index for i in memory
557!-- ime           end index for i in memory
558!-- jms           start index for j in memory
559!-- jme           end index for j in memory
560!-- ips           start index for i in patch
561!-- ipe           end index for i in patch
562!-- jps           start index for j in patch
563!-- jpe           end index for j in patch
564!-- kms           start index for k in memory
565!-- kme           end index for k in memory
566!-- i_start       start indices for i in tile
567!-- i_end         end indices for i in tile
568!-- j_start       start indices for j in tile
569!-- j_end         end indices for j in tile
570!-- kts           start index for k in tile
571!-- kte           end index for k in tile
572!-- num_tiles     number of tiles
573!
574!======================================================================
575
576    INTEGER,      INTENT(IN)                             ::                          &
577                                      ids,ide, jds,jde, kds,kde,                     &
578                                      ims,ime, jms,jme, kms,kme,                     &
579                                      ips,ipe, jps,jpe, kps,kpe,                     &
580                                                        kts,kte,                     &
581                                                      num_tiles
582
583    INTEGER, DIMENSION(num_tiles), INTENT(IN)            ::                          &
584       &           i_start,i_end,j_start,j_end
585
586
587    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: thetaValues,           &
588       &           pressvalues
589    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: tempValues
590
591! Local
592    REAL, PARAMETER                                      :: pres0 = 100000.
593    INTEGER                                              :: i,j,k
594
595!!    PRINT *,'  thetaValues: ', LBOUND(thetaValues),' x ', UBOUND(thetaValues)
596
597!!    DO k=kms,kme-2
598!!      DO j=jms,jme
599!!        DO i=ims,ime
600!!          IF (thetaValues(i,k,j) == 0. .AND. i > 0 .AND. i < 32 &
601!!            .AND. j > 0 .AND. j < 32) THEN
602!!            PRINT *,'  Zth: T ',i,k,j,thetaValues(i,k,j)
603!!          END IF
604!!        END DO
605!!      END DO
606!!    END DO
607
608    tempValues = ( thetaValues + 300.)*( pressValues / pres0 ) ** (r_d / cp)
609
610!!    DO k=kms,kme-2
611  END SUBROUTINE theta_to_temp
612
613  SUBROUTINE temp_to_theta1(tempValue, pressValue, thetaValue)
614! Subroutine to transform from temperature to potential temperature
615
616    USE module_model_constants
617
618    IMPLICIT NONE
619
620!!!!!!! Variables
621! thetaValue   potential temperature value
622! pressValue   pressure
623! tempValue    temperature value in K
624!
625    REAL, INTENT(IN)                                     :: tempValue, pressValue
626    REAL, INTENT(OUT)                                    :: thetaValue
627
628! Local
629    REAL, PARAMETER                                      :: pres0 = 100000.
630
631    thetaValue = tempValue * ( pres0 / pressValue ) ** (r_d / cp)
632
633  END SUBROUTINE temp_to_theta1
634
635  SUBROUTINE theta_to_temp1(thetaValue, pressValue, tempValue)
636! Subroutine to transform from potential temperature to temperature
637
638    USE module_model_constants
639
640    IMPLICIT NONE
641
642!!!!!!! Variables
643! thetaValue   potential temperature
644! pressValue   pressure
645! tempValue    temperature value in K
646
647    REAL, INTENT(IN)                                     :: thetaValue, pressvalue
648    REAL, INTENT(OUT)                                    :: tempValue
649
650! Local
651    REAL, PARAMETER                                      :: pres0 = 100000.
652
653    tempValue = thetaValue*( pressValue / pres0 ) ** (r_d / cp)
654
655  END SUBROUTINE theta_to_temp1
656
657  SUBROUTINE def_get_scalar_value(filen, valn, ds, di, dr, db)
658! Subroutine to get a scalar value from a .def file
659!   NOTE: For that scalar values with scientific notation (EN, or e-N) it will not work...
660!     next version!
661
662    IMPLICIT NONE
663
664    CHARACTER(LEN=100), INTENT(IN)                       :: filen, valn
665    CHARACTER(LEN=100), INTENT(OUT)                      :: ds
666    INTEGER, INTENT(OUT)                                 :: di
667    REAL, INTENT(OUT)                                    :: dr
668    LOGICAL, INTENT(OUT)                                 :: db
669
670! Local
671    INTEGER                                              :: funit, ios, Lvaln
672    INTEGER                                              :: iline, lineval
673    CHARACTER(LEN=1)                                     :: kindvalue
674    CHARACTER(LEN=100)                                   :: varname,varname0
675    CHARACTER(LEN=50)                                    :: subname, errmsg, value,  &
676      value0
677    LOGICAL                                              :: fexists, is_used
678
679!!!!!!! Variables
680! filen: Name of the def file to get the value
681! valn: name of the value to get
682! ds: string value (up to 100 characters)
683! di: integer value
684! dr: real value
685! db: logical value
686! value: value get from the file
687! kindvalue: kind of the value ('s': string, 'i': integer, 'r': real, 'b':boolean )
688
689    subname = 'def_get_scalar_value'
690    errmsg = 'ERROR -- error -- ERROR -- error'
691
692    INQUIRE(FILE=TRIM(filen), EXIST=fexists)
693
694    IF (fexists) THEN
695      Lvaln = LEN_TRIM(valn)
696
697      ds = 'None'
698      di = 0
699      dr = 0.
700      db = .FALSE.
701
702      DO funit=10,100
703        INQUIRE(UNIT=funit, OPENED=is_used)
704        IF (.NOT. is_used) EXIT
705      END DO
706      OPEN(funit, FILE=TRIM(filen), STATUS='old', FORM='formatted', IOSTAT=ios)
707      IF ( ios /= 0 ) THEN
708        PRINT *,errmsg
709        PRINT *, '  ' // TRIM(subname) // ": ERROR opening '" // TRIM(filen) // "'"
710        PRINT *,'ios: ',ios
711        STOP
712      END IF
713!      iunit = get_unused_unit()
714!      IF ( iunit <= 0 ) THEN
715!        IF ( wrf_dm_on_monitor() ) THEN
716!          CALL wrf_error_fatal('Error in module_ra_cam_support: could not find a free Fortran unit.')
717!        END IF
718!      END IF
719      lineval = 0
720      DO iline=1, 10000
721        READ(funit,*,END=125)varname0
722        varname = TRIM(varname0)
723        IF (varname(1:Lvaln) == TRIM(valn)) THEN
724          lineval = iline
725          EXIT
726        END IF
727      END DO
728
729 125  CONTINUE
730      IF (lineval == 0) THEN
731        PRINT *,errmsg
732        PRINT *,'  ' // TRIM(subname) // ": In file '" // TRIM(filen) //             &
733          & "' does not exist a variable called '" // TRIM(valn) // '" !!!'
734        STOP
735      ELSE
736        CLOSE(funit)
737        OPEN(funit, FILE=TRIM(filen), STATUS='old', FORM='formatted', IOSTAT=ios)
738        DO iline=1,lineval-1
739          READ(funit,*)varname
740        END DO
741        READ(funit,*)value0
742        CLOSE(UNIT=funit)
743
744        CALL string_split(value0, '=', 2, value)
745        CALL string_split(value, ':', 2, value)
746        CALL value_kind(value, kindvalue)
747
748!        PRINT *,'value: ',TRIM(value), ' kind: ',kindvalue
749
750        IF (kindvalue == 'i') THEN
751          READ(value, '(I50)')di
752        ELSE IF (kindvalue == 'r') THEN
753          READ(value, '(F50.25)')dr
754        ELSE IF (kindvalue == 'e') THEN
755          READ (value, '(E50.25)')dr
756        ELSE IF (kindvalue == 's') THEN
757          ds = TRIM(value)
758        ELSE
759          IF (TRIM(value) == 'y' .OR. TRIM(value) == 'yes' .OR. TRIM(value) == 'Y'   &
760            .OR. TRIM(value) == 'Yes' .OR. TRIM(value) == 'YES') db = .TRUE.
761
762          IF (TRIM(value) == 'n' .OR. TRIM(value) == 'no' .OR. TRIM(value) == 'N'    &
763            .OR. TRIM(value) == 'No' .OR. TRIM(value) == 'NO') db = .FALSE.
764        END IF
765      END IF
766    ELSE
767      PRINT *,errmsg
768      PRINT *,'  ' // TRIM(subname) // ": File '" // TRIM(filen) // "' does not exist !!"
769      STOP
770    END IF
771
772  END SUBROUTINE def_get_scalar_value
773
774  SUBROUTINE value_kind(value, kindval)
775! Subroutine to determine which kind of value is given
776
777    IMPLICIT NONE
778
779    CHARACTER(LEN=50), INTENT(IN)                        :: value
780    CHARACTER(LEN=1), INTENT(OUT)                        :: kindval
781
782! Local
783    INTEGER                                              :: ic
784    INTEGER                                              :: Lval, Nc, NnumVal, NstrVal
785    LOGICAL                                              :: is_point, is_exp
786
787!!!!!!! Variables
788! value: value to determine the cahracteristics
789! kindval: kind of value. 's': string, 'i': integer, 'r': real, 'e': real scientific notation,
790!   'b':boolean
791! Lval: length of characters of the value
792! NnumVal: number of numeric characers
793! NstrVal: number of string characters
794! is_point: boolean variable indicating the presence of a '.'
795! is_exp: boolean variable indicating the presence of a 'E' or 'e' (exponential)
796
797
798    IF (TRIM(value) == 'y' .OR. TRIM(value) == 'n' .OR. TRIM(value) == 'yes' .OR.    &
799      TRIM(value) == 'no' .OR. TRIM(value) == 'Y' .OR. TRIM(value) == 'N' .OR.       &
800      TRIM(value) == 'Yes' .OR. TRIM(value) == 'No' .OR. TRIM(value) == 'YES' .OR.   &
801      TRIM(value) == 'NO') THEN
802!      PRINT *,'    boolean!'
803      kindval = 'b'
804    ELSE
805      Lval = LEN_TRIM(value)
806      NnumVal = 0
807      NstrVal = 0
808      is_point = .FALSE.
809
810      DO ic=1,Lval
811        Nc = ICHAR(value(ic:ic))
812!        PRINT *, 'char: ',value(ic:ic), ' Nc:',Nc
813        IF (Nc >= 48 .AND. Nc <= 57) THEN
814          Nnumval = Nnumval + 1
815        ELSE
816          NstrVal = NstrVal + 1
817          IF (value(ic:ic) == '.') is_point = .TRUE.
818          IF ((value(ic:ic) == 'e') .OR. (value(ic:ic) == 'E')) is_exp = .TRUE.
819        END IF
820      END DO
821
822!      PRINT *,'NnumVal:', NnumVal, 'NstrVal: ', NstrVal, 'point: ', is_point
823      IF (Nnumval == Lval) THEN
824!        PRINT *,'    integer!'
825        kindval = 'i'
826      ELSE IF ( (Nnumval == Lval -1) .AND. (is_point .EQV. .TRUE.) .AND.             &
827        (NstrVal == 1)) THEN
828!        PRINT *,'    real!'
829        kindval = 'r'
830      ELSE IF ( ((Nnumval == Lval -2) .OR.(Nnumval == Lval -3)) .AND.                &
831        (is_point .EQV. .TRUE.) .AND. (is_exp .EQV. .TRUE.) .AND.                    &
832        ((NstrVal == 2) .OR. (NstrVal == 3)) ) THEN
833!        PRINT *,'    real exponential!'
834        kindval = 'e'
835      ELSE
836!        PRINT *,'    string!'
837        kindval = 's'
838      END IF
839    END IF
840
841  END SUBROUTINE value_kind
842
843  SUBROUTINE string_split(string, charsplit, Nsec, sec)
844! Subtoutine to split a string according to a character and returns the 'Nsec' of the split
845
846  IMPLICIT NONE
847
848  CHARACTER(LEN=50), INTENT(IN)                          :: string
849  CHARACTER(LEN=1), INTENT(IN)                           :: charsplit
850  INTEGER, INTENT(IN)                                    :: Nsec
851  CHARACTER(LEN=50), INTENT(OUT)                         :: sec
852
853! Local
854  INTEGER                                                :: ic
855  INTEGER                                                :: poschar, Lstring
856  INTEGER                                                :: isec, inisec
857  CHARACTER(LEN=50)                                      :: str, secchar
858
859!!!!!!! Variables
860! string: string to split
861! charsplit: character to use to split
862! Nsec: number of section to return
863! sec: returned section of the string
864
865!  PRINT *,'string: ',string
866
867  str = TRIM(string)
868  poschar=0
869
870  Lstring = LEN_TRIM(str)
871  isec = 1
872  poschar = INDEX(str, charsplit)
873  inisec = 1
874
875  IF (poschar == 0) THEN
876    sec = str
877
878  ELSE
879    DO ic=poschar,50 
880      secchar = str(inisec:poschar-1)
881      IF (isec == Nsec) sec = secchar
882      inisec = poschar + 1
883      poschar = INDEX(str(inisec:Lstring), charsplit)
884      IF (poschar == 0) poschar = Lstring + 2
885      isec = isec + 1
886    END DO
887
888  END IF
889
890  END SUBROUTINE string_split
891
892  SUBROUTINE get_lmdz(ddx,ddy,ddz,bdy,dlmdz,Nks,NzO,plon,plat,pzmasq,ppctsrf,pftsol, &
893    pftsoil,pqsol,pqsurf,pfalb1,pfalb2,pevap,psnow,pradsol,psolsw,psollw,pfder,      &
894    prain_fall,psnow_fall,pagesno,pzmea,pzstd,pzgam,pzthe,pzpic,pzval,prugoro,       &
895    prugos,pzmax0,pf0,pwake_s,pwake_cstar,pwake_pe,pwake_fip,psst,                   &
896    palbedo,pt_ancien,pq_ancien,pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,         &
897    pema_work1,pema_work2,pwake_deltat,pwake_deltaq,pfm_therm,pentr_therm,           &
898    pdetr_therm,pnat,pwrf_grid)
899! Subroutine to get LMDZ variables with values from WRF
900
901! WRF modules
902    USE module_domain_type
903
904! LMDZ modules
905    USE indice_sol_mod
906
907    IMPLICIT NONE
908
909    TYPE(domain) , TARGET                                :: pwrf_grid
910    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
911      dlmdz, Nks, NzO
912    REAL, DIMENSION(dlmdz), INTENT(IN)                   :: plon, plat, pzmasq,      &
913      pqsol, pradsol, psolsw, psollw, pfder, prain_fall, psnow_fall, pzmea, pzstd,   &
914      pzgam, pzthe, pzpic, pzval, prugoro, pzmax0, pf0, pwake_s,                     &
915      pwake_cstar, pwake_pe, pwake_fip, psst, palbedo, pnat
916! What to do with these variables? ,pnat,pbils,prug
917    REAL, DIMENSION(dlmdz,Nks), INTENT(IN)               :: ppctsrf,pftsol,pqsurf,   &
918      pfalb1,pfalb2,pevap,psnow,prugos,pagesno
919    REAL, DIMENSION(dlmdz,ddz), INTENT(IN)               :: pt_ancien,pq_ancien,     &
920      pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,pema_work1,pema_work2,pwake_deltat,&
921      pwake_deltaq,pentr_therm,pdetr_therm
922    REAL, DIMENSION(dlmdz,NzO,Nks), INTENT(IN)           :: pftsoil
923    REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN)             :: pfm_therm
924
925! Local
926    INTEGER                                              :: ix,iy,iz,ixy,dx,dy
927    CHARACTER(LEN=50)                                    :: LMDZvarmethod
928
929! Filling all 2D variables
930!!
931    PRINT *,'  Lluis before get_lmdz: 5 11 ltksoil: ',pwrf_grid%ltksoil(5,:,11),     &
932      ' ftsol: ', pftsol(315,:)
933
934    DO iy=1,ddy
935      DO ix=1,ddx
936        ixy=ddx*(iy-1)+ix
937! It does not change on time: pzmask
938! It might change due to ice dynamics?
939        pwrf_grid%lter(ix,iy) = ppctsrf(ixy,is_ter)
940        pwrf_grid%llic(ix,iy) = ppctsrf(ixy,is_lic)
941        pwrf_grid%loce(ix,iy) = ppctsrf(ixy,is_oce)
942        pwrf_grid%lsic(ix,iy) = ppctsrf(ixy,is_sic)
943        DO iz=1,Nks
944          pwrf_grid%ltksoil(ix,iz,iy) = pftsol(ixy,iz)
945        END DO
946        DO iz=1,NzO
947          pwrf_grid%lotter(ix,iz,iy) = pftsoil(ixy,iz,1)
948          pwrf_grid%lotlic(ix,iz,iy) = pftsoil(ixy,iz,2)
949          pwrf_grid%lotoce(ix,iz,iy) = pftsoil(ixy,iz,3)
950          pwrf_grid%lotsic(ix,iz,iy) = pftsoil(ixy,iz,4)
951        END DO
952        DO iz=1,Nks
953          pwrf_grid%lqksoil(ix,iz,iy) = pqsurf(ixy,iz)
954        END DO
955        pwrf_grid%lwsol(ix,iy) = pqsol(ixy)
956        DO iz=1,Nks
957          pwrf_grid%lalbksoil(ix,iz,iy) = pfalb1(ixy,iz)
958          pwrf_grid%llwalbksoil(ix,iz,iy) = pfalb2(ixy,iz)
959          pwrf_grid%levapksoil(ix,iz,iy) = pevap(ixy,iz)
960          pwrf_grid%lsnowksoil(ix,iz,iy) = psnow(ixy,iz)
961        END DO
962        pwrf_grid%lrads(ix,iy) = pradsol(ixy)
963        pwrf_grid%lsolsw(ix,iy) = psolsw(ixy)
964        pwrf_grid%lsollw(ix,iy) = psollw(ixy)
965        pwrf_grid%lfder(ix,iy) = pfder(ixy)
966        pwrf_grid%lrain(ix,iy) = prain_fall(ixy)
967        pwrf_grid%lsnow(ix,iy) = psnow_fall(ixy)
968        DO iz=1,Nks
969          pwrf_grid%lrugksoil(ix,iz,iy) = prugos(ixy,iz)
970          pwrf_grid%lagesnoksoil(ix,iz,iy) = pagesno(ixy,iz)
971        END DO
972! These variables do not change on time!
973!!
974!    pzmea, pzstd, pzsig, pzgam, pzthe, pzpic, pzval, prugsrel
975        pwrf_grid%lrugsea(ix,iy) = prugos(ixy,is_oce)
976!        pwrf_grid%lrunofflic(ix,iy)= prun_off_lic(ixy)
977        pwrf_grid%lzmax0(ix,iy) = pzmax0(ixy)
978        pwrf_grid%lf0(ix,iy) = pf0(ixy)
979        pwrf_grid%lwake_s(ix,iy) = pwake_s(ixy)
980        pwrf_grid%lwake_cstar(ix,iy) = pwake_cstar(ixy)
981        pwrf_grid%lwake_pe(ix,iy) = pwake_pe(ixy)
982        pwrf_grid%lwake_fip(ix,iy) = pwake_fip(ixy)
983        pwrf_grid%lnat(ix,iy) = pnat(ixy)
984        pwrf_grid%sst(ix,iy) = psst(ixy)
985!        pwrf_grid%lbils(ix,iy) = pbils(ixy)
986        pwrf_grid%albedo(ix,iy) = palbedo(ixy)
987!        pwrf_grid%lrug(ix,iy) = prug(ixy)
988      END DO
989    END DO
990
991    dx=ddx+2*bdy
992    dy=ddy+2*bdy
993
994    CALL vect_mat(pclwcon, dx, dy, ddz, bdy, pwrf_grid%lclwcon)
995    CALL vect_mat(prnebcon, dx, dy, ddz, bdy, pwrf_grid%lrnebcon)
996    CALL vect_mat(pratqs, dx, dy, ddz, bdy, pwrf_grid%lratqs)
997    CALL vect_mat(pema_work1, dx, dy, ddz, bdy, pwrf_grid%lema_work1)
998    CALL vect_mat(pema_work2, dx, dy, ddz, bdy, pwrf_grid%lema_work2)
999    CALL vect_mat(pwake_deltat, dx, dy, ddz, bdy, pwrf_grid%lwake_deltat)
1000    CALL vect_mat(pwake_deltaq, dx, dy, ddz, bdy, pwrf_grid%lwake_deltaq)
1001    CALL vect_mat(pfm_therm, dx, dy, ddz+1, bdy, pwrf_grid%lfm_therm)
1002    CALL vect_mat(pentr_therm, dx, dy, ddz, bdy, pwrf_grid%lentr_therm)
1003    CALL vect_mat(pdetr_therm, dx, dy, ddz, bdy, pwrf_grid%ldetr_therm)
1004
1005    PRINT *,'  Lluis after get_lmdz: 5 11 ltksoil: ',pwrf_grid%ltksoil(5,:,11),      &
1006      ' ftsol: ', pftsol(315,:)
1007
1008    RETURN
1009
1010  END SUBROUTINE get_lmdz
1011
1012  SUBROUTINE get_lowbdy(ddx, ddy, ddz, bdy, dlmdz, Nks, Nz0, pkglo, pwrf_grid,       &
1013    ppctsrf, pftsol, prugos, palbedo, psst, ptsurf_lim, pz0_lim, palb_lim, prug_glo, &
1014    ppct_glo)
1015! Subroutine to load LMDZ limit variables with values from lowbdy WRF
1016
1017! WRF modules
1018!    USE module_model_constants
1019    USE module_domain_type
1020
1021! LMDZ modules
1022    USE indice_sol_mod
1023
1024    IMPLICIT NONE
1025
1026    TYPE(domain) , TARGET                                :: pwrf_grid
1027    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
1028      dlmdz, Nks, Nz0, pkglo
1029    REAL, DIMENSION(dlmdz), INTENT(INOUT)                :: psst, palbedo,           &
1030      ptsurf_lim, pz0_lim, palb_lim
1031! What to do with these variables? ,pnat,pbils
1032    REAL, DIMENSION(dlmdz,Nks), INTENT(INOUT)            :: ppctsrf, pftsol, prugos
1033    REAL, DIMENSION(pkglo), INTENT(OUT)                  :: prug_glo
1034    REAL, DIMENSION(pkglo,Nks), INTENT(OUT)              :: ppct_glo
1035
1036! Local
1037    INTEGER                                              :: ix,iy,iz,ixy,dx,dy,ixx,iyy
1038    CHARACTER(LEN=50)                                    :: LMDZvarmethod,fname,errmsg
1039
1040    fname='get_lowbdy'
1041    errmsg='ERROR -- error -- ERROR -- error'
1042
1043! Filling all 2D variables
1044!!
1045    DO iy=1,ddy
1046      DO ix=1,ddx
1047        ixy=ddx*(iy-1)+ix
1048        ppctsrf(ixy,is_ter) = pwrf_grid%lter(ix,iy)
1049        ppctsrf(ixy,is_lic) = pwrf_grid%llic(ix,iy)
1050        IF ( pwrf_grid%xice(ix,iy) /= pwrf_grid%lsic(ix,iy) ) THEN
1051! No ocean/frozen ocean fraction before
1052          IF (pwrf_grid%loce(ix,iy) == 0.) pwrf_grid%lrugksoil(ix,3,iy) = 0.0001
1053          IF (pwrf_grid%lsic(ix,iy) == 0.) pwrf_grid%lrugksoil(ix,4,iy) = 0.001
1054          pwrf_grid%loce(ix,iy) = 1. - pwrf_grid%lsic(ix,iy)
1055          pwrf_grid%lsic(ix,iy) = pwrf_grid%lsic(ix,iy)
1056          ppctsrf(ixy,is_oce) = pwrf_grid%loce(ix,iy)
1057          ppctsrf(ixy,is_sic) = pwrf_grid%lsic(ix,iy)
1058! Recomputing total roughness
1059!!
1060          pwrf_grid%lrug(ix,iy)=pwrf_grid%lter(ix,iy)*pwrf_grid%lrugksoil(ix,1,iy) +&
1061            pwrf_grid%llic(ix,iy)*pwrf_grid%lrugksoil(ix,2,iy) +                    &
1062            pwrf_grid%loce(ix,iy)*pwrf_grid%lrugksoil(ix,3,iy) +                    &
1063            pwrf_grid%lsic(ix,iy)*pwrf_grid%lrugksoil(ix,4,iy)
1064        END IF
1065! Full SST grid points must have a SST value from the wrfinput (WRF 1/0 landmask)
1066!   but that LMDZ points with fractions of land/mask could not. Using the closest one
1067!   if it is an isolated grid point, using TSK instead
1068        IF (pwrf_grid%sst(ix,iy) < 200.15) THEN
1069          pwrf_grid%ltksoil(ix,3,iy) = -9.
1070          IF ( (ix > 1).AND.(ix < ddx) .AND. (iy > 1).AND.(iy < ddy) ) THEN
1071            DO ixx=-1,1
1072              DO iyy=-1,1
1073!               PRINT *,ixx,iyy,wrf_grid%sst(ix+ixx,iy+iyy)
1074                IF (pwrf_grid%sst(ix+ixx,iy+iyy) > 200.15) THEN
1075                  pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%sst(ix+ixx,iy+iyy)
1076                  EXIT
1077                END IF
1078              END DO
1079              IF ( pwrf_grid%ltksoil(ix,3,iy) > 200.15) EXIT
1080            END DO
1081            IF ( pwrf_grid%ltksoil(ix,3,iy) == -9.) pwrf_grid%ltksoil(ix,3,iy) =     &
1082              pwrf_grid%tsk(ix,iy)
1083            ELSE
1084              pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%tsk(ix,iy)
1085            END IF
1086        ELSE
1087          pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%sst(ix,iy)
1088        END IF
1089        IF ( (pwrf_grid%ltksoil(ix,3,iy) < 200.15) .OR.                              &
1090          (pwrf_grid%ltksoil(ix,3,iy) > 370.) ) THEN
1091          PRINT *,TRIM(errmsg)
1092!          WRITE(wrf_err_message,*) '  ' // TRIM(fname) //                            &
1093!            ' SST inconsistency at point: ', ix, ', ', iy,' with loce: ',            &
1094!            wrf_grid%loce(ix,iy),' has a tsoil(oce)= ',                              &
1095!            wrf_grid%ltksoil(ix,3,iy),' sst: ',wrf_grid%sst(ix,iy),' tsk: ',         &
1096!            wrf_grid%tsk(ix,iy)
1097!          CALL wrf_error_fatal(TRIM(wrf_err_message))
1098          PRINT *, '  ' // TRIM(fname) //                                            &
1099          ' SST inconsistency at point: ', ix, ', ', iy,' with loce: ',              &
1100            pwrf_grid%loce(ix,iy),' has a tsoil(oce)= ',                             &
1101            pwrf_grid%ltksoil(ix,3,iy),' sst: ', pwrf_grid%sst(ix,iy),' tsk: ',      &
1102            pwrf_grid%tsk(ix,iy)
1103        END IF
1104        psst(ixy) = pwrf_grid%ltksoil(ix,3,iy)
1105        palbedo(ixy) = pwrf_grid%albbck(ix,iy)
1106        DO iz=1,Nks
1107          prugos(ixy,iz) = pwrf_grid%lrugksoil(ix,iz,iy)
1108        END DO
1109
1110! In ocean_forced_mod is used as tsurf_lim = sst(knindex(1:knon)). Let's complicate it
1111!   a little bit...
1112        ptsurf_lim(ixy)= psst(ixy)
1113! In surf_land_bucket_mod is used as z0_new=rugos(knindex(1:knon),1) [z0_new == z0_limit]
1114        pz0_lim(ixy) = prugos(ixy,1)
1115        palb_lim(ixy) = palbedo(ixy)
1116
1117        prug_glo(ixy) = prugos(ixy,is_ter)
1118        ppct_glo(ixy,:) = ppctsrf(ixy,:)
1119
1120      END DO
1121    END DO
1122
1123    RETURN
1124
1125  END SUBROUTINE get_lowbdy
1126
1127  SUBROUTINE load_lmdz(pwrf_grid,ddx,ddy,ddz,bdy,dlmdz,Nks,NzO,prlon,prlat,pzmasq,   &
1128    ppctsrf,pftsol,pftsoil,pqsurf,pqsol,pfalb1,pfalb2,pevap,psnow,pradsol,psolsw,    &
1129    psollw,pfder,prain_fall,psnow_fall,pagesno,pzmea,pzstd,pzgam,pzthe,pzpic,pzval,  &
1130    prugoro,prugos,prun_off_lic,pzmax0,pf0,pwake_s,pwake_cstar,pwake_pe,pwake_fip,   &
1131    psst,palbedo,pt_ancien,pq_ancien,pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,    &
1132    pema_work1,pema_work2,pwake_deltat,pwake_deltaq,pfm_therm,pentr_therm,           &
1133    pdetr_therm,pnat)
1134! Subroutine to load LMDZ variables with values from WRF
1135
1136! WRF modules
1137!    USE module_model_constants
1138    USE module_domain_type
1139
1140! LMDZ modules
1141    USE indice_sol_mod
1142
1143    IMPLICIT NONE
1144
1145    TYPE(domain) , TARGET                                :: pwrf_grid
1146    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
1147      dlmdz, Nks, NzO
1148    REAL, DIMENSION(dlmdz), INTENT(OUT)                  :: prlon, prlat, pzmasq,    &
1149      pqsol, pradsol, psolsw, psollw, pfder, prain_fall, psnow_fall, pzmea, pzstd,   &
1150      pzgam, pzthe, pzpic, pzval, prugoro, prun_off_lic, pzmax0, pf0, pwake_s,       &
1151      pwake_cstar, pwake_pe, pwake_fip, psst, palbedo,pnat
1152    REAL, DIMENSION(dlmdz,Nks), INTENT(OUT)              :: ppctsrf,pftsol,pqsurf,   &
1153      pfalb1,pfalb2,pevap,psnow,prugos,pagesno
1154    REAL, DIMENSION(dlmdz,ddz), INTENT(OUT)              :: pt_ancien,pq_ancien,     &
1155      pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,pema_work1,pema_work2,pwake_deltat,&
1156      pwake_deltaq,pentr_therm,pdetr_therm
1157    REAL, DIMENSION(dlmdz,NzO,Nks), INTENT(OUT)          :: pftsoil
1158    REAL, DIMENSION(dlmdz,ddz+1), INTENT(OUT)            :: pfm_therm
1159
1160
1161! Local
1162    INTEGER                                              :: ix,iy,iz,ixy,dx,dy
1163    CHARACTER(LEN=50)                                    :: LMDZvarmethod,fname
1164
1165    fname='load_lmdz'
1166
1167! Filling all 2D variables
1168!!
1169    PRINT *,'  ' // TRIM(fname) // '; LUBOUND(grid%lmsk): ',LBOUND(pwrf_grid%lmsk), UBOUND(pwrf_grid%lmsk)
1170    PRINT *,'  ' // TRIM(fname) // '; ddx ddy: ',ddy,ddx
1171
1172    DO iy=1,ddy
1173      DO ix=1,ddx
1174        ixy=ddx*(iy-1)+ix
1175        prlon(ixy) = pwrf_grid%xlong(ix,iy)
1176        prlat(ixy) = pwrf_grid%xlat(ix,iy)
1177        pzmasq(ixy) = pwrf_grid%lmsk(ix,iy)
1178        ppctsrf(ixy,is_ter) = pwrf_grid%lter(ix,iy)
1179        ppctsrf(ixy,is_lic) = pwrf_grid%llic(ix,iy)
1180        ppctsrf(ixy,is_oce) = pwrf_grid%loce(ix,iy)
1181        ppctsrf(ixy,is_sic) = pwrf_grid%lsic(ix,iy)
1182        DO iz=1,Nks
1183          pftsol(ixy,iz) = pwrf_grid%ltksoil(ix,iz,iy)
1184        END DO
1185        DO iz=1,NzO
1186          pftsoil(ixy,iz,1) = pwrf_grid%lotter(ix,iz,iy)
1187          pftsoil(ixy,iz,2) = pwrf_grid%lotlic(ix,iz,iy)
1188          pftsoil(ixy,iz,3) = pwrf_grid%lotoce(ix,iz,iy)
1189          pftsoil(ixy,iz,4) = pwrf_grid%lotsic(ix,iz,iy)
1190        END DO
1191        DO iz=1,Nks
1192          pqsurf(ixy,iz) = pwrf_grid%lqksoil(ix,iz,iy)
1193        END DO
1194        pqsol(ixy) = pwrf_grid%lwsol(ix,iy)
1195        DO iz=1,Nks
1196          pfalb1(ixy,iz) = pwrf_grid%lalbksoil(ix,iz,iy)
1197          pfalb2(ixy,iz) = pwrf_grid%llwalbksoil(ix,iz,iy)
1198          pevap(ixy,iz) = pwrf_grid%levapksoil(ix,iz,iy)
1199          psnow(ixy,iz) = pwrf_grid%lsnowksoil(ix,iz,iy)
1200        END DO
1201! No way to get back the independent radiative values (unless using LMDZ 'hist*'
1202!   files). Use it in the LMDZ's added way also from grid (no need to initialize it,
1203!   no radiation on t=0)
1204        pradsol(ixy) = pwrf_grid%lrads(ix,iy)
1205        psolsw(ixy) = pwrf_grid%lsolsw(ix,iy)
1206        psollw(ixy) = pwrf_grid%lsollw(ix,iy)
1207        pfder(ixy) = pwrf_grid%lfder(ix,iy)
1208! LMDZ does not provide a way (in the restart) to distingusih between convective,
1209!    non-convective. Thus Use it as a whole
1210        prain_fall(ixy) = pwrf_grid%lrain(ix,iy)
1211        psnow_fall(ixy) = pwrf_grid%lsnow(ix,iy)
1212        DO iz=1,Nks
1213          prugos(ixy,iz) = pwrf_grid%lrugksoil(ix,iz,iy)
1214          pagesno(ixy,iz) = pwrf_grid%lagesnoksoil(ix,iz,iy)
1215        END DO
1216        pzmea(ixy) = pwrf_grid%lzmea(ix,iy)
1217        pzstd(ixy) = pwrf_grid%lzstd(ix,iy)
1218        pzgam(ixy) = pwrf_grid%lzgam(ix,iy)
1219        pzthe(ixy) = pwrf_grid%lzthe(ix,iy)
1220        pzpic(ixy) = pwrf_grid%lzpic(ix,iy)
1221        pzval(ixy) = pwrf_grid%lzval(ix,iy)
1222        prugoro(ixy) = pwrf_grid%lzrugsrel(ix,iy)
1223!        prugos(ixy,is_oce) = pwrf_grid%lrugsea(ix,iy)
1224        prun_off_lic(ixy) = pwrf_grid%lrunofflic(ix,iy)
1225        pzmax0(ixy) = pwrf_grid%lzmax0(ix,iy)
1226        pf0(ixy) = pwrf_grid%lf0(ix,iy)
1227        pwake_s(ixy) = pwrf_grid%lwake_s(ix,iy)
1228        pwake_cstar(ixy) = pwrf_grid%lwake_cstar(ix,iy)
1229        pwake_pe(ixy) = pwrf_grid%lwake_pe(ix,iy)
1230        pwake_fip(ixy) = pwrf_grid%lwake_fip(ix,iy)
1231!! limit.nc
1232! nat == pctsrf, not needed ?
1233        pnat(ixy) = pwrf_grid%lnat(ix,iy)
1234        psst(ixy) = pwrf_grid%sst(ix,iy)
1235! not used, not defined....
1236!        pbils(ixy) = pwrf_grid%lbils(ix,iy)
1237        palbedo(ixy) = pwrf_grid%albedo(ix,iy)
1238! already done in previous steps
1239!        prugos(ixy) = pwrf_grid%lrug(ix,iy)
1240
1241      END DO
1242    END DO
1243
1244    dx=ddx+2*bdy
1245    dy=ddy+2*bdy
1246
1247    CALL mat_vect(pwrf_grid%t_2, dx, dy, ddz, bdy, pt_ancien)
1248! New variable added in the Registry qv_2
1249    CALL mat_vect(pwrf_grid%qv_2, dx, dy, ddz, bdy, pq_ancien)
1250    CALL mat_vect(pwrf_grid%u_2, dx, dy, ddz, bdy, pu_ancien)
1251    CALL mat_vect(pwrf_grid%v_2, dx, dy, ddz, bdy, pv_ancien)
1252    CALL mat_vect(pwrf_grid%lclwcon, dx, dy, ddz, bdy, pclwcon)
1253    CALL mat_vect(pwrf_grid%lrnebcon, dx, dy, ddz, bdy, prnebcon)
1254    CALL mat_vect(pwrf_grid%lratqs, dx, dy, ddz, bdy, pratqs)
1255    CALL mat_vect(pwrf_grid%lema_work1, dx, dy, ddz, bdy, pema_work1)
1256    CALL mat_vect(pwrf_grid%lema_work2, dx, dy, ddz, bdy, pema_work2)
1257    CALL mat_vect(pwrf_grid%lwake_deltat, dx, dy, ddz, bdy, pwake_deltat)
1258    CALL mat_vect(pwrf_grid%lwake_deltaq, dx, dy, ddz, bdy, pwake_deltaq)
1259    CALL mat_vect(pwrf_grid%lfm_therm, dx, dy, ddz+1, bdy, pfm_therm)
1260    CALL mat_vect(pwrf_grid%lentr_therm, dx, dy, ddz, bdy, pentr_therm)
1261    CALL mat_vect(pwrf_grid%ldetr_therm, dx, dy, ddz, bdy, pdetr_therm)
1262
1263    RETURN
1264
1265  END SUBROUTINE load_lmdz
1266
1267  SUBROUTINE vect_mat(vect, dx, dy, dz, wbdy, mat)
1268! Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
1269
1270    IMPLICIT NONE
1271
1272    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1273    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(IN) :: vect
1274    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1), INTENT(OUT) :: mat
1275
1276! Local
1277    INTEGER                                              :: i,j,k
1278    INTEGER                                              :: ddx, ddy
1279    INTEGER                                              :: lp, il, jl
1280
1281! Variables
1282! d[x/y/z]: dimensions
1283! mat: matrix to transform
1284! vect: vector to produce
1285    mat=0.
1286
1287    ddx=dx-2*wbdy
1288    ddy=dy-2*wbdy
1289!    lp=321
1290!    jl=INT(lp/ddx)+1
1291!    il=lp-ddx*(jl-1)
1292
1293    DO k=1,dz
1294      DO j=1,ddy
1295        DO i=1,ddx
1296          mat(i,k,j)=vect(ddx*(j-1)+i,k)
1297        END DO
1298      END DO
1299!      PRINT *,k,'  Lluis vect: ',vect(lp,k),mat(il,k,jl),' lp: ',lp,' il jl :',il,jl
1300    END DO
1301
1302  END SUBROUTINE vect_mat
1303
1304
1305  SUBROUTINE vect_mat_zstagg(vect, dx, dy, dz, wbdy, mat)
1306! Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
1307
1308    IMPLICIT NONE
1309
1310    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1311    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(IN) :: vect
1312    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz,1-wbdy:dy-wbdy+1), INTENT(OUT) :: mat
1313
1314! Local
1315    INTEGER                                              :: i,j,k
1316    INTEGER                                              :: ddx, ddy
1317
1318! Variables
1319! d[x/y/z]: dimensions
1320! mat: matrix to transform
1321! vect: vector to produce
1322
1323    ddx=dx-2*wbdy
1324    ddy=dy-2*wbdy
1325
1326    mat=0.
1327
1328    DO k=1,dz
1329      DO j=1,ddy
1330        DO i=1,ddx
1331          mat(i,k,j)=vect(ddx*(j-1)+i,k)
1332        END DO
1333      END DO
1334    END DO
1335
1336  END SUBROUTINE vect_mat_zstagg
1337
1338
1339  SUBROUTINE mat_vect(mat, dx, dy, dz, wbdy, vect)
1340! Subroutine to transform a 3D matrix to a LMDZ physic 1D vector
1341
1342    IMPLICIT NONE
1343
1344    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1345    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1), INTENT(IN) :: mat
1346    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(OUT) ::                      &
1347      vect
1348
1349! Local
1350    INTEGER                                              :: i,j,k
1351    INTEGER                                              :: ddx,ddy
1352
1353! Variables
1354! d[x/y/z]: dimensions
1355! mat: matrix to transform
1356! vect: vector to produce
1357
1358! Removing HALO
1359    ddx=dx-2*wbdy
1360    ddy=dy-2*wbdy
1361
1362    vect=0.
1363
1364    DO k=1,dz
1365      DO j=1,ddy
1366        DO i=1,ddx
1367          vect(ddx*(j-1)+i,k)= mat(i,k,j)
1368        END DO
1369      END DO
1370    END DO
1371
1372  END SUBROUTINE mat_vect
1373
1374  SUBROUTINE mat_vect_zstagg(mat, dx, dy, dz, wbdy, vect)
1375! Subroutine to transform a 3D matrix z-staggered to a LMDZ physic 1D vector
1376
1377    IMPLICIT NONE
1378
1379    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1380    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz,1-wbdy:dy-wbdy+1), INTENT(IN) :: mat
1381    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(OUT) ::                      &
1382      vect
1383
1384! Local
1385    INTEGER                                              :: i,j,k
1386    INTEGER                                              :: ddx,ddy
1387
1388! Variables
1389! d[x/y/z]: dimensions
1390! mat: matrix to transform
1391! vect: vector to produce
1392
1393! Removing HALO
1394    ddx=dx-2*wbdy
1395    ddy=dy-2*wbdy
1396
1397    DO k=1,dz
1398      DO j=1,ddy
1399        DO i=1,ddx
1400          vect(ddx*(j-1)+i,k)= mat(i,k,j)
1401        END DO
1402      END DO
1403    END DO
1404
1405  END SUBROUTINE mat_vect_zstagg
1406
1407  SUBROUTINE WRFlanduse_LMDZKsoil(lndcn,is,ie,js,je,dx,dy,ncat,lndu,mask,kt,kl,ko,ks)
1408! Subroutine to retrieve LMDZ Ksoil classification using WRF landuse
1409
1410    IMPLICIT NONE
1411
1412    CHARACTER(LEN=50), INTENT(IN)                        :: lndcn
1413    INTEGER, INTENT(IN)                                  :: is,ie,js,je,dx,dy,ncat
1414    REAL, DIMENSION(is:ie,ncat,js:je), INTENT(IN)        :: lndu
1415    REAL, DIMENSION(is:ie,js:je), INTENT(IN)             :: mask
1416    REAL, DIMENSION(is:ie,js:je), INTENT(OUT)            :: kt,kl,ko,ks
1417
1418! Local
1419    INTEGER                                              :: i,j,k,ij,Nnones,Ncattot
1420    CHARACTER(LEN=50)                                    :: fname, errmsg, wrnmsg
1421    REAL, DIMENSION(is:ie,js:je)                         :: totlndu
1422    INTEGER, ALLOCATABLE, DIMENSION(:)                   :: ter_wk, lic_wk, oce_wk,  &
1423      sic_wk
1424    LOGICAL, ALLOCATABLE, DIMENSION(:,:)                 :: wk
1425    INTEGER                                              :: Nter_wk, Nlic_wk,        &
1426        Noce_wk, Nsic_wk, klks_equal
1427    REAL*8, DIMENSION(is:ie,js:je)                       :: Ksoilval
1428
1429
1430!!!!!!!! Variables
1431! lndcn: land-use data-set (same as in VEGPARM.TBL)
1432! dx,dy: horizontal dimension of the domain
1433! ncat: number of categories of the data-set
1434! lndu: WRF landuse matrix
1435! mask: land-sea mask from WRF
1436! k[t/l/o/s]: LMDZ ter, lic, oce, sic matrices
1437! Ncattot: theoretical number of laduses for a given data-base
1438! [ter/lic/oce/sic]_wk: vectors with the indices of the landuse for LMDZ Ksoil
1439! N[ter/lic/oce/sic]_wk: number of indices of the landuse for LMDZ Ksoil
1440! wk: boolean matrix for each Ksoil indicating 'true' if the WRF landuse is for Ksoil
1441
1442    fname="'WRFlanduse_LMDZKsoil'"
1443    errmsg='ERROR -- error -- ERROR -- error'
1444    wrnmsg='WARNING -- warning -- WARNING -- warning'
1445   
1446! Checking if at all the grid points SUM(landuse) == 1.
1447!!
1448    PRINT *,'  ' // TRIM(fname) // " using '" // TRIM(lndcn) // "' WRF land types"
1449
1450    totlndu = SUM(lndu, 2)
1451
1452    Nnones=0
1453    DO i=1,dx
1454      DO j=1,dy
1455        IF (totlndu(i,j) /= 1. ) Nnones = Nnones+1
1456      END DO
1457    END DO
1458
1459    IF (Nnones /= 0) THEN
1460      PRINT *,TRIM(errmsg)
1461      PRINT *,'    ' // TRIM(fname) // ': ', Nnones,' grid points do not have a ' // &
1462        'TOTlanduse == 1. from WRF !!!!!'
1463      DO i=1,dx
1464        DO j=1,dy
1465          IF (totlndu(i,j) /= 1.) PRINT *,'      ',i,', ',j,' total land use: ',     &
1466            totlndu(i,j),' indiv. values :', lndu(i,:,j)
1467        END DO
1468      END DO
1469    END IF
1470
1471    SELECT CASE (lndcn)
1472      CASE ('OLD')
1473        Ncattot = 13
1474        Nter_wk = 11
1475        Nlic_wk = 1
1476        Noce_wk = 1
1477        Nsic_wk = 1
1478
1479        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1480        ALLOCATE(ter_wk(Nter_wk))
1481        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1482        ALLOCATE(lic_wk(Nlic_wk))
1483        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1484        ALLOCATE(oce_wk(Noce_wk))
1485        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1486        ALLOCATE(sic_wk(Nsic_wk))
1487
1488! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1489!!
1490!  1,'Urban land'                    2,'Agriculture'
1491!  3,'Range-grassland'               4,'Deciduous forest'
1492!  5,'Coniferous forest'             6,'Mixed forest and wet land'
1493!  7,'Water'                         8,'Marsh or wet land'
1494!  9,'Desert'                       10,'Tundra'
1495! 11,'Permanent ice'                12,'Tropical or subtropical forest'
1496! 13,'Savannah'
1497
1498        ter_wk = (/ 1,2,3,4,5,6,8,9,10,12,13 /)
1499        lic_wk = (/ 11 /)
1500        oce_wk = (/ 7 /)
1501        sic_wk = (/ 11 /)
1502
1503      CASE ('USGS')
1504        Ncattot = 33
1505        Nter_wk = 22
1506        Nlic_wk = 1
1507        Noce_wk = 1
1508        Nsic_wk = 1
1509
1510        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1511        ALLOCATE(ter_wk(Nter_wk))
1512        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1513        ALLOCATE(lic_wk(Nlic_wk))
1514        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1515        ALLOCATE(oce_wk(Noce_wk))
1516        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1517        ALLOCATE(sic_wk(Nsic_wk))
1518
1519! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1520!!
1521! 1,'Urban and Built-Up Land'        2,'Dryland Cropland and Pasture'
1522! 3,'Irrigated Cropland and Pasture' 4,'Mixed Dryland/Irrigated Cropland and Pasture'
1523! 5,'Cropland/Grassland Mosaic'      6,'Cropland/Woodland Mosaic'
1524! 7,'Grassland'                      8,'Shrubland'
1525! 9,'Mixed Shrubland/Grassland'     10,'Savanna'
1526!11,'Deciduous Broadleaf Forest'    12,'Deciduous Needleleaf Forest'
1527!13,'Evergreen Broadleaf Forest'    14,'Evergreen Needleleaf Forest'
1528!15,'Mixed Forest'                  16,'Water Bodies'
1529!17,'Herbaceous Wetland'            18,'Wooded Wetland'
1530!19,'Barren or Sparsely Vegetated'  20,'Herbaceous Tundra'
1531!21,'Wooded Tundra'                 22,'Mixed Tundra'
1532!23,'Bare Ground Tundra'            24,'Snow or Ice'
1533!25,'Playa'                         26,'Lava'
1534!27,'White Sand'                    28,'Unassigned'
1535!29,'Unassigned'                    30,'Unassigned'
1536!31,'Low Intensity Residential '    32,'High Intensity Residential'
1537!33,'Industrial or Commercial'     
1538
1539! Correct values
1540!        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,18,19,20,21,22,23,25,26,  &
1541!          27,31,32,33 /)
1542! How ever, USGS by default has 24 so...
1543        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,18,19,20,21,22,23 /)
1544        lic_wk = (/ 24 /)
1545        oce_wk = (/ 16 /)
1546        sic_wk = (/ 24 /)
1547
1548      CASE ('MODIFIED_IGBP_MODIS_NOAH')
1549        Ncattot = 33
1550        Nter_wk = 28
1551        Nlic_wk = 1
1552        Noce_wk = 1
1553        Nsic_wk = 1
1554
1555        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1556        ALLOCATE(ter_wk(Nter_wk))
1557        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1558        ALLOCATE(lic_wk(Nlic_wk))
1559        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1560        ALLOCATE(oce_wk(Noce_wk))
1561        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1562        ALLOCATE(sic_wk(Nsic_wk))
1563
1564! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1565!!
1566!  1,'Evergreen Needleleaf Forest'             2,'Evergreen Broadleaf Forest'
1567!  3,'Deciduous Needleleaf Forest'             4,'Deciduous Broadleaf Forest'
1568!  5,'Mixed Forests'                           6,'Closed Shrublands'
1569!  7,'Open Shrublands'                         8,'Woody Savannas'
1570!  9,'Savannas'                               10,'Grasslands'
1571! 11,'Permanent wetlands'                     12,'Croplands'
1572! 13,'Urban and Built-Up'                     14,'cropland/natural vegetation mosaic'
1573! 15,'Snow and Ice'                           16,'Barren or Sparsely Vegetated'
1574! 17,'Water'                                  18,'Wooded Tundra'
1575! 19,'Mixed Tundra'                           20,'Barren Tundra'
1576! 21,'Unassigned'                             22,'Unassigned'
1577! 23,'Unassigned'                             24,'Unassigned'
1578! 25,'Unassigned'                             26,'Unassigned'
1579! 27,'Unassigned'                             28,'Unassigned'
1580! 29,'Unassigned'                             30,'Unassigned'
1581! 31,'Low Intensity Residential '             32,'High Intensity Residential'
1582! 33,'Industrial or Commercial'
1583
1584        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,19,20,31,32,33 /)
1585        lic_wk = (/ 15 /)
1586        oce_wk = (/ 17 /)
1587        sic_wk = (/ 15 /)
1588
1589      CASE ('SiB')
1590        Ncattot = 16
1591        Nter_wk = 14
1592        Nlic_wk = 1
1593        Noce_wk = 1
1594        Nsic_wk = 1
1595
1596        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1597        ALLOCATE(ter_wk(Nter_wk))
1598        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1599        ALLOCATE(lic_wk(Nlic_wk))
1600        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1601        ALLOCATE(oce_wk(Noce_wk))
1602        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1603        ALLOCATE(sic_wk(Nsic_wk))
1604
1605! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1606!!
1607!  1,'Evergreen Broadleaf Trees'        2,'Broadleaf Deciduous Trees'
1608!  3,'Deciduous and Evergreen Trees'    4,'Evergreen Needleleaf Trees'
1609!  5,'Deciduous Needleleaf Trees'       6,'Ground Cover with Trees and Shrubs'
1610!  7,'Groundcover Only'              8,'Broadleaf Shrubs with Perennial Ground Cover'
1611!  9,'Broadleaf Shrubs with Bare Soil' 10,'Groundcover with Dwarf Trees and Shrubs'
1612! 11,'Bare Soil'                       12,'Agriculture or C3 Grassland'
1613! 13,'Persistent Wetland'              14,'Dry Coastal Complexes'
1614! 15,'Water'                           16,'Ice Cap and Glacier'
1615
1616        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14 /)
1617        lic_wk = (/ 16 /)
1618        oce_wk = (/ 15 /)
1619        sic_wk = (/ 16 /)
1620
1621      CASE ('LW12')
1622        Ncattot = 3
1623        Nter_wk = 1
1624        Nlic_wk = 1
1625        Noce_wk = 1
1626        Nsic_wk = 1
1627
1628        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1629        ALLOCATE(ter_wk(Nter_wk))
1630        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1631        ALLOCATE(lic_wk(Nlic_wk))
1632        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1633        ALLOCATE(oce_wk(Noce_wk))
1634        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1635        ALLOCATE(sic_wk(Nsic_wk))
1636
1637! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1638!!
1639!  1,'Land'                              2,'Water'
1640!  3,'Snow or Ice'
1641
1642        ter_wk = (/ 1 /)
1643        lic_wk = (/ 3 /)
1644        oce_wk = (/ 2 /)
1645        sic_wk = (/ 3 /)
1646
1647      CASE DEFAULT
1648        PRINT *,TRIM(errmsg)
1649        PRINT *,'  ' // TRIM(fname) // ": landuse data-base '" // TRIM(lndcn) //     &
1650          "' not ready !!!"
1651        STOP
1652
1653    END SELECT
1654
1655! Checking the correct number of categories
1656    IF (Ncattot /= ncat) THEN
1657      PRINT *,TRIM(wrnmsg)
1658      PRINT *,'  ' // TRIM(fname) // ": '" // TRIM(lndcn) // "' landuse data-base "  &
1659        // ' has ',Ncattot,' land-uses types ',ncat,' passed !!!!'
1660!      STOP
1661    END IF
1662
1663! Mimic WRF landuse with 4 LMDZ Ksoils
1664!!
1665    IF (ALLOCATED(wk)) DEALLOCATE(wk)
1666!    ALLOCATE(wk(4,ncattot))
1667    ALLOCATE(wk(4,ncat))
1668
1669    wk = .FALSE.
1670
1671    DO i=1,Nter_wk
1672      wk(1,ter_wk(i)) = .TRUE.
1673    END DO
1674    DO i=1,Nlic_wk
1675      wk(2,lic_wk(i)) = .TRUE.
1676    END DO
1677    DO i=1,Noce_wk
1678      wk(3,oce_wk(i)) = .TRUE.
1679    END DO
1680    DO i=1,Nsic_wk
1681      wk(4,sic_wk(i)) = .TRUE.
1682    END DO
1683
1684    klks_equal = 0
1685    DO i=1,dx
1686      DO j=1,dy
1687! Using double precission numbers to avoid truncation errors...       
1688        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(1,:))
1689        kt(i,j) = Ksoilval(1,1)
1690        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(2,:))
1691        kl(i,j) = Ksoilval(1,1)
1692        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(3,:))
1693        ko(i,j) = Ksoilval(1,1)
1694        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(4,:))
1695        ks(i,j) = Ksoilval(1,1)
1696
1697        IF (kl(i,j) == ks(i,j)) klks_equal = klks_equal+1
1698      END DO
1699    END DO
1700
1701! Assigning lic/sic to all that ice/snow according to the percentage of ter/oce
1702!   Is there any other better way? Satellites, do not know what is below the ice ?!
1703
1704    IF (klks_equal == dx*dy) THEN
1705!   Assuming that kl=ks=ice
1706      DO i=1,dx
1707        DO j=1,dy
1708          IF ((kt(i,j)+ko(i,j)) /= 0.) THEN
1709            Ksoilval(1,1) = kl(i,j)*kt(i,j)/(kt(i,j)+ko(i,j))
1710            kl(i,j) = Ksoilval(1,1)
1711            Ksoilval(1,1) = ks(i,j)*ko(i,j)/(kt(i,j)+ko(i,j))
1712            ks(i,j) = Ksoilval(1,1)
1713          ELSE
1714            Ksoilval(1,1) = kl(i,j)/(kl(i,j)+ks(i,j))
1715            Ksoilval(1,2) = ks(i,j)/(kl(i,j)+ks(i,j))
1716            kl(i,j) = Ksoilval(1,1)
1717            ks(i,j) = Ksoilval(1,2)
1718          END IF
1719! Introducing WRF land-sea mask
1720          IF (mask(i,j) == 1.) THEN
1721            kl(i,j)=kl(i,j)+ks(i,j)
1722            ks(i,j)=0.
1723          ELSE
1724            ks(i,j)=kl(i,j)+ks(i,j)
1725            kl(i,j)=0.
1726          END IF
1727        END DO
1728      END DO
1729    ELSE
1730      PRINT *,'  ' // TRIM(fname) // ': land-use data provided different lic & sic !'
1731    END IF
1732
1733! Checking if at all the grid points kt+kl+ko+ks == 1.
1734!!
1735    totlndu = kt+kl+ko+ks
1736
1737    Nnones=0
1738    DO i=1,dx
1739      DO j=1,dy
1740        IF (totlndu(i,j) /= 1. ) Nnones = Nnones+1
1741        IF (kt(i,j) + ko(i,j) == 0.)THEN
1742          PRINT *,TRIM(wrnmsg)
1743          PRINT *,'  ' // TRIM(fname) // ': point without land or water category !!'
1744          PRINT *,'i j WRF_ter ter ______________'
1745          DO k=1,Nter_wk
1746            PRINT *,i,j,lndu(i,ter_wk(k),j), kt(i,j)
1747          END DO
1748          PRINT *,'i j WRF_oce oce ______________'
1749          DO k=1,Noce_wk
1750            PRINT *,i,j,lndu(i,oce_wk(k),j), ko(i,j)
1751          END DO
1752          PRINT *,'  wrf%landusef: ', lndu(i,:,j)
1753!          STOP
1754        END IF
1755      END DO
1756    END DO
1757
1758    ij=1
1759    IF (Nnones /= 0) THEN
1760      PRINT *,TRIM(errmsg)
1761      PRINT *,'  ' // TRIM(fname) // ': ', Nnones,' grid points do not have a ' //   &
1762        'TOTlanduse == 1.!!!!!'
1763      PRINT *,'pt i j ter+lic+oce+sic : ter lic oce sic ______________'
1764      DO i=1,dx
1765        DO j=1,dy
1766          IF (totlndu(i,j) /= 1.) THEN
1767            PRINT *,ij,'; ',i,j,totlndu(i,j),' : ',kt(i,j),kl(i,j),ko(i,j),ks(i,j)
1768            ij=ij+1
1769          END IF
1770        END DO
1771      END DO
1772    END IF
1773
1774    DEALLOCATE(ter_wk)
1775    DEALLOCATE(lic_wk)
1776    DEALLOCATE(oce_wk)
1777    DEALLOCATE(sic_wk)
1778    DEALLOCATE(wk)
1779
1780    RETURN
1781
1782  END SUBROUTINE WRFlanduse_LMDZKsoil
1783
1784  SUBROUTINE interpolate_layers(NvalsA,posA,valsA,NvalsB,posB,valsB)
1785! Subroutine to interpolate values between different layers (constant values along the depth of a layer)
1786
1787    IMPLICIT NONE
1788
1789    INTEGER, INTENT(IN)                                  :: NvalsA,NvalsB
1790    REAL, DIMENSION(NvalsA), INTENT(IN)                  :: posA,valsA
1791    REAL, DIMENSION(NvalsB), INTENT(IN)                  :: posB
1792    REAL, DIMENSION(NvalsB), INTENT(OUT)                 :: valsB
1793
1794! Local
1795    INTEGER                                              :: i,k
1796    REAL                                                 :: xA,yA,xB,yB
1797
1798!!!!!!! Variables
1799! A: original values used to interpolate
1800! B: values to obtain
1801! Nvals[A/B]: number of values
1802! pos[A/B]: depths of each layer
1803! vals[A/B]: values of each layer
1804
1805!!!!!!! Functions/Subroutines
1806! interpolate_lin: subroutine to linearly interpolate
1807
1808    i=1
1809    DO k=1,NvalsB
1810      IF (posB(k) > posA(i)) i = i+1
1811      IF (i-1 < 1) THEN
1812        xA = posA(i)
1813        yA = valsA(i)
1814        xB = posA(i+1)
1815        yB = valsA(i+1)
1816      ELSEIF ( i > NvalsA) THEN
1817        xA = posA(i-2)
1818        yA = valsA(i-2)
1819        xB = posA(i-1)
1820        yB = valsA(i-1)
1821      ELSE
1822        xA = posA(i-1)
1823        yA = valsA(i-1)
1824        xB = posA(i)
1825        yB = valsA(i)
1826      END IF
1827      IF (posB(k) <= xB .AND. i-1 > 1) THEN
1828        CALL interpolate_lin(xA,yA,xB,yB,posB(k),valsB(k))
1829      ELSEIF (posB(k) <= xA .AND. i-1 < 1) THEN
1830        valsB(k) = yA + (yB - yA)*(posB(k) - xA)/(xB - xA)
1831      ELSE
1832        valsB(k) = yB + (yB - yA)*(posB(k) - xB)/(xB - xA)
1833      END IF
1834!      PRINT *,'i ',i,'  Interpolating for: ',posB(k),' with xA: ', xA,' yA: ',yA,' & ',   &
1835!       ' xB: ', xB,' yB: ',yB,' --> ',valsB(k)
1836    END DO
1837
1838  END SUBROUTINE interpolate_layers
1839
1840  SUBROUTINE interpolate_lin(x1,y1,x2,y2,x0,y)
1841! Program to interpolate values y=a+b*x with: (from clWRF modifications)
1842!       a=y1
1843!       b=(y2-y1)/(x2-x1)
1844!       x=(x2-x0)
1845
1846    IMPLICIT NONE
1847
1848    REAL, INTENT (IN)                                    :: x1,x2,x0
1849!    REAL(r8), INTENT (IN)                                :: y1,y2
1850!    REAL(r8), INTENT (OUT)                               :: y
1851!    REAL(r8)                                             :: a,b,x
1852    REAL, INTENT (IN)                                   :: y1,y2
1853    REAL, INTENT (OUT)                                  :: y
1854    REAL                                                :: a,b,x
1855
1856    a=y1
1857    b=(y2-y1)/(x2-x1)
1858
1859    IF (x0 .GE. x1) THEN
1860      x=x0-x1
1861    ELSE
1862      x=x1-x0
1863      b=-b
1864    ENDIF
1865
1866    y=a+b*x
1867   
1868!    PRINT *,'y=a+b*x'
1869!    IF (b .LT. 0.) THEN
1870!      PRINT *,'y=',y1,' ',b,'*x'
1871!    ELSE
1872!      PRINT *,'y=',y1,'+',b,'*x'
1873!    ENDIF
1874!    PRINT *,'y(',x,')=',y
1875
1876  END SUBROUTINE interpolate_lin
1877
1878  SUBROUTINE table_characteristics(datasetn,tablen,datasetline,Nvar,Nr,Nc)
1879! Subroutine to read values from a WRF .TBL
1880
1881    IMPLICIT NONE
1882
1883    CHARACTER(LEN=50), INTENT(IN)                        :: datasetn
1884    CHARACTER(LEN=100), INTENT(IN)                       :: tablen
1885    INTEGER, INTENT(OUT)                                 :: datasetline, Nvar, Nr, Nc
1886
1887! Local
1888    INTEGER                                              :: il, iline, ic
1889    CHARACTER(LEN=50)                                    :: fname, errmsg
1890    CHARACTER(LEN=200)                                   :: lineval
1891    LOGICAL                                              :: existsfile, is_used
1892    INTEGER                                              :: funit, Llineval, icoma
1893
1894!!!!!! Variables
1895! datasetn: name of the data set to use
1896! tablen: name of the table to read
1897! Nvar: number of variations (g.e. seasons)
1898! Nr: table row number
1899! Nc: table column number
1900! datasetline: line in the table where the dataset starts
1901
1902    fname='table_characteristics'
1903    errmsg='ERROR -- error -- ERROR -- error'
1904
1905    INQUIRE(FILE=TRIM(tablen), EXIST=existsfile)
1906
1907    IF (existsfile) THEN
1908      DO funit=10,100
1909        INQUIRE(UNIT=funit, OPENED=is_used)
1910        IF (.NOT. is_used) EXIT
1911      END DO
1912
1913    ELSE
1914      PRINT *,TRIM(errmsg)
1915      PRINT *,TRIM(fname),': file "'// TRIM(tablen) // '" does not exist !!!'
1916      STOP
1917    END IF
1918
1919    OPEN(UNIT=funit, FILE=TRIM(tablen), STATUS='OLD')
1920
1921! Looking for the place where data-set starts
1922!!
1923    datasetline = -1
1924    iline=1
1925    DO il=1,100000
1926      READ(funit,*,END=100)lineval
1927      IF (TRIM(lineval) == TRIM(datasetn)) datasetline = il
1928    END DO
1929
1930100 CONTINUE
1931    IF (datasetline == -1) THEN
1932      PRINT *,TRIM(errmsg)
1933      PRINT *,TRIM(fname),': in file "'// TRIM(tablen) // '" dataset: "'//           &
1934        TRIM(datasetn) // '" does not exist !!!'
1935      STOP
1936    END IF
1937
1938! Looking for number of values
1939!!
1940    REWIND(UNIT=funit)
1941    DO ic=1,datasetline
1942      READ(funit,*)lineval
1943    END DO
1944    READ(funit,*)Nr,Nvar
1945
1946    IF (Nvar > 1) READ(funit,*)lineval
1947
1948    READ(funit,'(200A)')lineval
1949    Llineval=LEN(TRIM(lineval))
1950
1951!    PRINT *,'  lineval: ******' // TRIM(lineval) // '******', Llineval
1952
1953    Nc = 0
1954    ic = 1
1955    DO WHILE (ic <= Llineval)
1956!      PRINT *,lineval(ic:Llineval)
1957      icoma = SCAN(lineval(ic:Llineval),',')
1958      IF (icoma > 1) THEN
1959        Nc = Nc + 1
1960!        PRINT *,'  Nc: ',Nc,' ic ',ic,' icoma: ',icoma
1961      ic = ic + icoma
1962      ELSE
1963        ic = ic + Llineval
1964      END IF
1965    END DO
1966    PRINT *,'  dataset :"' // TRIM(datasetn) // '" at line ',datasetline,            &
1967      ' in table :"'// TRIM(tablen) // '" has N variations: ',Nvar,                  &
1968      ' Number of types: ',Nr, ' Number of values: ',Nc
1969
1970    CLOSE(funit)
1971
1972    RETURN
1973
1974  END SUBROUTINE table_characteristics
1975
1976  SUBROUTINE read_table(datasetn,tablen,datasetline,Nvar,Nr,Nc,vals,valns)
1977! Subroutine to read values from a WRF .TBL
1978
1979    IMPLICIT NONE
1980
1981    CHARACTER(LEN=50), INTENT(IN)                        :: datasetn
1982    CHARACTER(LEN=100), INTENT(IN)                       :: tablen
1983    INTEGER, INTENT(IN)                                  :: datasetline, Nvar, Nr, Nc
1984    REAL, DIMENSION(Nvar,Nr,Nc), INTENT(OUT)             :: vals
1985    CHARACTER(LEN=50), DIMENSION(Nr), INTENT(OUT)        :: valns
1986
1987
1988! Local
1989    INTEGER                                              :: iv, ir, ic
1990    CHARACTER(LEN=50)                                    :: fname, errmsg, line
1991    LOGICAL                                              :: existsfile, is_used
1992    INTEGER                                              :: funit
1993
1994!!!!!!! Variables
1995! datasetn: name of the data set to use
1996! tablen: name of the table to read
1997! Nvar: number of variations (g.e. seasons)
1998! Nr: table row number
1999! Nc: table column number
2000! datasetline: line in the table where the dataset starts
2001! vals: values of the data-set
2002! valns: assigned name of each value in the dataset
2003
2004    fname='read_table'
2005    errmsg='ERROR -- error -- ERROR -- error'
2006
2007    INQUIRE(FILE=TRIM(tablen), EXIST=existsfile)
2008
2009    IF (existsfile) THEN
2010      DO funit=10,100
2011        INQUIRE(UNIT=funit, OPENED=is_used)
2012        IF (.NOT. is_used) EXIT
2013      END DO
2014
2015    ELSE
2016      PRINT *,TRIM(errmsg)
2017        PRINT *,TRIM(fname),': file "'// TRIM(tablen) // '" does not exist !!!'
2018        STOP
2019    END IF
2020
2021    OPEN(UNIT=funit, FILE=TRIM(tablen), STATUS='OLD')
2022
2023    DO ic=1,datasetline+1
2024      READ(funit,*)line
2025    END DO
2026
2027    DO iv=1, Nvar
2028      IF (Nvar > 1) READ(funit,*)line
2029      DO ir=1, Nr
2030        READ(funit,*)(vals(iv,ir,ic), ic=1,Nc),valns(ir)
2031      END DO
2032    END DO
2033
2034!    PRINT *,'  values dataset :"' // TRIM(datasetn) // '" in table :"' //TRIM(tablen)
2035!    DO iv=1,Nvar
2036!      PRINT *,'   variation: ',iv
2037!      DO ir=1,Nr
2038!        PRINT *,'    ',(vals(iv,ir,ic), ic=1,Nc),valns(ir)
2039!      END DO
2040!    END DO
2041
2042    CLOSE(funit)
2043
2044    RETURN
2045
2046  END SUBROUTINE read_table
2047
2048
2049  SUBROUTINE compute_soil_mass(Nvar,Nr,Nc,vals,valns,soilt,Nnost,NOsoilt,smois,ddx,  &
2050    ddy,ddz,soilm,soilq)
2051! Subroutine to compute the total soil mass (kg) from a point with different
2052!   percentages of soil types
2053
2054    IMPLICIT NONE
2055
2056    INTEGER, INTENT(IN)                                  :: Nvar, Nr, Nc, Nnost
2057    REAL, DIMENSION(Nvar,Nr,Nc), INTENT(IN)              :: vals
2058    CHARACTER(LEN=50), DIMENSION(Nr), INTENT(IN)         :: valns
2059    REAL, DIMENSION(Nr), INTENT(IN)                      :: soilt
2060    CHARACTER(LEN=50), DIMENSION(Nnost), INTENT(IN)      :: NOsoilt
2061    REAL, INTENT(IN)                                     :: smois, ddx, ddy, ddz
2062    REAL, INTENT(OUT)                                    :: soilm,soilq
2063
2064! Local
2065    INTEGER                                              :: iv, ir, ic
2066    CHARACTER(LEN=50)                                    :: fname, errmsg
2067    LOGICAL                                              :: is_land
2068
2069!!!!!!! Variables
2070! Nvar: number of variations (g.e. seasons)
2071! Nr: table row number
2072! Nc: table column number
2073! vals: values of the data-set (col2: density [g cm-3])
2074! valns: assigned name of each value in the dataset
2075! soilt: vector with the percentages of soil types
2076! Nnost: number of soil types which are not land
2077! NOsoilt: labels of the non-land soil types
2078! smois: soil moisture [m3 m-3]
2079! ddx,ddy,ddz: volume of the grid cell [m]
2080! soilm: total soil mass [kg]
2081! soilq: total soil water [kg kg-1]
2082
2083    fname='read_table'
2084    errmsg='ERROR -- error -- ERROR -- error'
2085
2086    soilm = 0.
2087    soilq = 0.
2088
2089!    PRINT *,'  soilmoisture: ',smois
2090    DO ir=1,Nr
2091      is_land = .TRUE.
2092      IF (soilt(ir) /= 0.) THEN
2093        DO iv=1, Nnost
2094          IF (TRIM(valns(ir)) == TRIM(NOsoilt(iv))) is_land = .FALSE.
2095        END DO
2096        IF (is_land) THEN
2097!          PRINT *,'  soil moisture range: ', vals(1,ir,3), vals(1,ir,5)
2098          soilm = soilm + soilt(ir)*vals(1,ir,2)*ddx*ddy*ddz*1000.
2099          soilq = soilq + soilt(ir)*vals(1,ir,2)*smois*ddx*ddy*ddz*1000.
2100!          PRINT *,'  Percentage of soil "' // TRIM(valns(ir)) //'": ',soilt(ir)*100.,&
2101!            ' % density: ',vals(1,ir,2),' partial mass: ',                           &
2102!            soilt(ir)*vals(1,ir,2)*ddx*ddy*ddz*1000. , ' partial water: ',           &
2103!            soilt(ir)*vals(1,ir,2)*smois*ddx*ddy*ddz*1000.
2104        END IF
2105      END IF
2106    END DO
2107! Remember 1 kg H20 = 1 l = 1mm / m2
2108!!  PRINT *,'  TOTAL. soil mass: ',soilm,' [kg] soil water: ',soilq,' [kg]',           &
2109!!    ' equiv. water depth in soil ',soilq/(ddx*ddy*1000.),' [m]'
2110
2111  END SUBROUTINE compute_soil_mass
2112
2113  SUBROUTINE compute_TOTsoil_mass_water(datname, tabname, Nsoiltypes, soiltypest,    &
2114    soiltypesb, Nsoillayers, smois, soillayers, sizeX, sizeY, soilTOTmass,           &
2115    soilTOThumidity)
2116! Subroutine to compute total soil mass and water from a given grid point using a
2117!   given data-set from SOILPARM.TBL
2118
2119      IMPLICIT NONE
2120
2121      CHARACTER(LEN=50), INTENT(IN)                      :: datname
2122      CHARACTER(LEN=100), INTENT(IN)                     :: tabname
2123      INTEGER, INTENT(IN)                                :: Nsoiltypes, Nsoillayers
2124      REAL, DIMENSION(Nsoiltypes), INTENT(IN)            :: soiltypest, soiltypesb
2125      REAL, DIMENSION(Nsoillayers), INTENT(IN)           :: smois, soillayers
2126      REAL, INTENT(IN)                                   :: sizeX, sizeY
2127      REAL, INTENT(OUT)                                  :: soilTOTmass,             &
2128         soilTOThumidity
2129
2130! Local
2131      INTEGER                                            :: ix,iy,iz
2132      REAL, ALLOCATABLE, DIMENSION(:,:,:)                :: values
2133      CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:)       :: namevalues,NOsoiltypes
2134      INTEGER                                            :: linedataset,Nvariations, &
2135        Nrows, Ncols
2136      INTEGER                                            :: NNOsoiltypes
2137      REAL, ALLOCATABLE, DIMENSION(:)                    :: soiltypes
2138      REAL                                               :: dz,soilmass,soilhumidity
2139      CHARACTER(LEN=50)                                  :: fname, errmsg
2140
2141!!!!!!! Variables
2142! datname: name of the data set to use
2143! tabname: name of the table to read
2144! values: mmatrix with the name of the values: variation (Nvariations, g.e. seasons),
2145!   table row (Nrows),  table column (Ncols)
2146! linedataset: line inside data-set where values start
2147! namevalues: string assigned to each table row
2148! Nsoiltypes: number of soil types
2149! soiltypes[t/b]: percentage of the soil types at a given grid point on top/bottom
2150! soiltypes: percentage of the soil types at a given grid point (lineray interpolated)
2151! NOsoiltypes: vector with the names of soil types which are not land
2152
2153      fname='"compute_TOTsoil_mass_water"'
2154      errmsg='ERROR -- error -- ERROR -- error'
2155
2156      CALL table_characteristics(datname, tabname, linedataset, Nvariations,         &
2157        Nrows, Ncols)
2158
2159      IF (Nrows /= Nsoiltypes) THEN
2160        PRINT *,errmsg
2161        PRINT '  ' // TRIM(fname) // ': Number of soil-types: ',Nsoiltypes,          &
2162          ' does not concide with the number of values: ', Nrows,                    &
2163          ' from the data-set "' // TRIM(datname) // '" !!!'
2164        STOP
2165      END IF
2166
2167      IF (ALLOCATED(values)) DEALLOCATE(values)
2168      ALLOCATE(values(Nvariations,Nrows,Ncols))
2169
2170      IF (ALLOCATED(namevalues)) DEALLOCATE(namevalues)
2171      ALLOCATE(namevalues(Nrows))
2172
2173      CALL read_table(datname, tabname, linedataset, Nvariations, Nrows, Ncols,      &
2174        values, namevalues)
2175
2176      IF (ALLOCATED(soiltypes)) DEALLOCATE(soiltypes)
2177      ALLOCATE(soiltypes(Nrows))
2178
2179! On WRF3.3 SOILPARM.TBL
2180      SELECT CASE (TRIM(datname))
2181
2182        CASE ('STAS')
2183          NNOsoiltypes = 1
2184          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2185          ALLOCATE(NOsoiltypes(NNOsoiltypes))
2186
2187          NOsoiltypes=(/ 'WATER' /)
2188        CASE ('STAS-RUC')
2189          NNOsoiltypes = 1
2190          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2191          ALLOCATE(NOsoiltypes( NNOsoiltypes))
2192
2193          NOsoiltypes=(/ 'WATER' /)
2194      END SELECT
2195
2196      soilTOTmass = 0.
2197      soilTOThumidity = 0.
2198
2199      DO iz=1,Nsoillayers
2200        IF (iz == 1 ) THEN
2201          dz = soillayers(iz)
2202        ELSE
2203          dz = soillayers(iz) - soillayers(iz - 1)
2204        END IF
2205        IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN
2206          DO ix=1,Nrows
2207            CALL interpolate_lin(soillayers(1), soiltypest(ix),                      &
2208              soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz), soiltypes(ix))
2209          END DO
2210        ELSEIF (iz == 1) THEN
2211          soiltypes = soiltypest
2212        ELSE
2213          soiltypes = soiltypesb
2214        END IF
2215!        PRINT *,'lev: ',iz,' dz: ',dz,' soil types: ' ,soiltypes
2216
2217        CALL compute_soil_mass(Nvariations,Nrows,Ncols,values,namevalues,soiltypes,  &
2218          NNOsoiltypes,NOsoiltypes,smois(iz),1.,1.,dz,soilmass,soilhumidity)
2219        soilTOTmass = soilTOTmass + soilmass
2220        soilTOThumidity = soilTOThumidity + soilhumidity
2221      END DO
2222
2223      DEALLOCATE(values)
2224      DEALLOCATE(namevalues)
2225      DEALLOCATE(soiltypes)
2226      DEALLOCATE(NOsoiltypes)
2227
2228  END SUBROUTINE compute_TOTsoil_mass_water
2229
2230  SUBROUTINE compute_TOTsoil_mass_water_values(datname, Nvariations, Nrows, Ncols,   &
2231    values, namevalues, ip, jp, Nsoiltypes, soiltypest, soiltypesb, Nsoillayers,     &
2232    smois, soillayers, sizeX, sizeY, soilTOTmass, soilTOThumidity)
2233! Subroutine to compute total soil mass and water from a given grid point using the
2234!   values from a given data-set of SOILPARM.TBL
2235
2236      IMPLICIT NONE
2237
2238      CHARACTER(LEN=50), INTENT(IN)                      :: datname
2239      INTEGER, INTENT(IN)                                :: Nvariations, Nrows, Ncols
2240      REAL, DIMENSION(Nvariations, Nrows, Ncols), INTENT(IN) :: values
2241      CHARACTER(LEN=50), DIMENSION(Nrows), INTENT(IN)    :: namevalues
2242      INTEGER, INTENT(IN)                                :: ip, jp
2243      INTEGER, INTENT(IN)                                :: Nsoiltypes, Nsoillayers
2244      REAL, DIMENSION(Nsoiltypes), INTENT(IN)            :: soiltypest, soiltypesb
2245      REAL, DIMENSION(Nsoillayers), INTENT(IN)           :: smois, soillayers
2246      REAL, INTENT(IN)                                   :: sizeX, sizeY
2247      REAL, INTENT(OUT)                                  :: soilTOTmass,             &
2248         soilTOThumidity
2249
2250! Local
2251      INTEGER                                            :: ix,iy,iz
2252      CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:)       :: NOsoiltypes
2253      INTEGER                                            :: linedataset
2254      INTEGER                                            :: NNOsoiltypes
2255      REAL, ALLOCATABLE, DIMENSION(:)                    :: soiltypes
2256      REAL                                               :: dz,soilmass,soilhumidity
2257      CHARACTER(LEN=50)                                  :: fname, errmsg
2258
2259!!!!!!! Variables
2260! datname: name of the data set to use
2261! values: mmatrix with the name of the values: variation (Nvariations, g.e. seasons),
2262!   table row (Nrows),  table column (Ncols)
2263! linedataset: line inside data-set where values start
2264! namevalues: string assigned to each table row
2265! ip, jp: coordinates of the grid point
2266! Nsoiltypes: number of soil types
2267! soiltypes[t/b]: percentage of the soil types at a given grid point on top/bottom
2268! soiltypes: percentage of the soil types at a given grid point (lineray interpolated)
2269! NOsoiltypes: vector with the names of soil types which are not land
2270
2271      fname='"compute_TOTsoil_mass_water_values"'
2272      errmsg='ERROR -- error -- ERROR -- error'
2273
2274      IF (Nrows /= Nsoiltypes) THEN
2275        PRINT *,errmsg
2276        PRINT '  ' // TRIM(fname) // ': Number of soil-types: ',Nsoiltypes,          &
2277          ' does not concide with the number of values: ', Nrows,                    &
2278          ' from the data-set "' // TRIM(datname) // '" !!!'
2279        STOP
2280      END IF
2281
2282      IF (ALLOCATED(soiltypes)) DEALLOCATE(soiltypes)
2283      ALLOCATE(soiltypes(Nrows))
2284
2285! On WRF3.3 SOILPARM.TBL
2286      SELECT CASE (TRIM(datname))
2287
2288        CASE ('STAS')
2289          NNOsoiltypes = 1
2290          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2291          ALLOCATE(NOsoiltypes(NNOsoiltypes))
2292
2293          NOsoiltypes=(/ 'WATER' /)
2294        CASE ('STAS-RUC')
2295          NNOsoiltypes = 1
2296          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2297          ALLOCATE(NOsoiltypes( NNOsoiltypes))
2298
2299          NOsoiltypes=(/ 'WATER' /)
2300      END SELECT
2301
2302      soilTOTmass = 0.
2303      soilTOThumidity = 0.
2304
2305      DO iz=1,Nsoillayers
2306        IF (iz == 1 ) THEN
2307          dz = soillayers(iz)
2308        ELSE
2309          dz = soillayers(iz) - soillayers(iz - 1)
2310        END IF
2311        IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN
2312          DO ix=1,Nrows
2313            CALL interpolate_lin(soillayers(1), soiltypest(ix),                      &
2314              soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz), soiltypes(ix))
2315          END DO
2316        ELSEIF (iz == 1) THEN
2317          soiltypes = soiltypest
2318        ELSE
2319          soiltypes = soiltypesb
2320        END IF
2321!        PRINT *,'lev: ',iz,' dz: ',dz,' soil types: ' ,soiltypes
2322
2323        CALL compute_soil_mass(Nvariations,Nrows,Ncols,values,namevalues,soiltypes,  &
2324          NNOsoiltypes,NOsoiltypes,smois(iz),sizeX,sizeY,dz,soilmass,soilhumidity)
2325        soilTOTmass = soilTOTmass + soilmass
2326        soilTOThumidity = soilTOThumidity + soilhumidity
2327      END DO
2328
2329      IF ( (soilTOTmass < 0.) .OR. (soilTOThumidity < 0.) ) THEN
2330        PRINT *,errmsg
2331        PRINT *,'  ' // TRIM(fname) // ' at point ',ip, ', ', jp,                    &
2332          ' negative total soil mass: ',soilTOTmass, ' or soil water content',       &
2333          soilTOThumidity
2334        PRINT *,'  Soil_name density soil_types_________'
2335        DO ix=1, Nrows
2336           PRINT *,ix, TRIM(namevalues(ix)), values(1,ix,2), soiltypes(ix)
2337        END DO
2338        PRINT *,'  depth_layer partial_mass soil+moisure partial_moisture__________'
2339        DO iz=1,Nsoillayers
2340          IF (iz == 1 ) THEN
2341            dz = soillayers(iz)
2342          ELSE
2343            dz = soillayers(iz) - soillayers(iz - 1)
2344          END IF
2345          IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN
2346            DO ix=1,Nrows
2347              CALL interpolate_lin(soillayers(1), soiltypest(ix),                    &
2348                soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz),soiltypes(ix))
2349            END DO
2350          ELSEIF (iz == 1) THEN
2351            soiltypes = soiltypest
2352          ELSE
2353            soiltypes = soiltypesb
2354          END IF
2355          PRINT *,'  ',iz, dz, values(1,iz,2)*sizeX*sizeY*dz*1000., smois(iz),       &
2356             smois(iz)*values(1,iz,2)*sizeX*sizeY*dz*1000.
2357        END DO
2358        STOP
2359      END IF
2360
2361
2362      DEALLOCATE(soiltypes)
2363      DEALLOCATE(NOsoiltypes)
2364
2365  END SUBROUTINE compute_TOTsoil_mass_water_values
2366
2367! cat Registry_out.LMDZ | awk '{print  "        pwrf_grid%"tolower($3)" "$4" = p"tolower($3)" "$4}' >> get_lmdz_out.f90
2368! cat get_lmdz_out.f90n' ',' '{print $3}' | tr '\n
2369  SUBROUTINE get_lmdz_out(ddx, ddy, ddz, bdy, dlmdz, Nks, Nz0, nqtot, ivap, iliq,    &
2370    iflag_pbl, iflag_con, iflag_thermals, iflag_wake, read_climoz, itap,             &
2371    pdtphys, paprs, pplay, pphi, pphis, presnivs, qx,                                &
2372    d_u, d_v, d_t, d_qx, d_ps,                                                       &
2373    pwrf_grid)
2374! Subroutine to get all the LMDZ output variables into the WRF Registry structure
2375
2376! WRF modules
2377    USE module_domain_type
2378
2379! LMDZ modules
2380    USE indice_sol_mod
2381    USE phys_state_var_mod
2382    USE phys_output_var_mod
2383!    USE pbl_surface_mod
2384    USE phys_local_var_mod
2385    USE comgeomphy
2386
2387! WRF_LMDZ modules
2388    USE output_lmdz_NOmodule
2389
2390! netcdf
2391    USE netcdf, ONLY: nf90_fill_real
2392
2393    IMPLICIT NONE
2394
2395    TYPE(domain) , TARGET                                :: pwrf_grid
2396    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
2397      dlmdz, Nks, Nz0, nqtot, ivap, iliq, iflag_pbl, iflag_con, iflag_thermals,      &
2398      iflag_wake, read_climoz, itap
2399    REAL, INTENT(IN)                                     :: pdtphys
2400    REAL, DIMENSION(dlmdz), INTENT(IN)                   :: pphis, d_ps
2401    REAL, DIMENSION(ddz), INTENT(IN)                     :: presnivs
2402    REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN)             :: paprs
2403    REAL, DIMENSION(dlmdz,ddz), INTENT(IN)               :: pplay, pphi, d_u, d_v, d_t
2404    REAL, DIMENSION(dlmdz,ddz,nqtot), INTENT(IN)         :: qx, d_qx
2405
2406!    EXTERNAL suphel
2407
2408! Local
2409    INTEGER                                              :: dx, dy, ix, iy, ixy, iz
2410    INTEGER                                              :: i,j,k,l,n
2411    INTEGER                                              :: nsrf
2412    REAL, DIMENSION(dlmdz,ddz)                           :: var3d
2413    REAL, DIMENSION(dlmdz,ddz+1)                         :: var3dstagg
2414    REAL                                                 :: RDAY, RG, RD, RMO3
2415    CHARACTER(LEN=50)                                    :: fname, errmsg
2416    CHARACTER(LEN=4)                                     :: bb2
2417    REAL, DIMENSION(dlmdz)                               :: zx_tmp_fi2d ! variable temporaire grille physique
2418    REAL, DIMENSION(dlmdz,ddz)                           :: zx_tmp_fi3d ! variable temporaire pour champs 3D
2419    REAL, DIMENSION(dlmdz,ddz+1)                         :: zx_tmp_fi3d1 !variable temporaire pour champs 3D (kelvp1)
2420    REAL(KIND=8), DIMENSION(dlmdz,ddz)                   :: zx_tmp2_fi3d ! variable temporaire pour champs 3D
2421    REAL, PARAMETER                                      :: dobson_u = 2.1415e-05
2422
2423!    klon=dldmz
2424!    klev=ddz
2425
2426
2427!! To think about / discuss with Frederic
2428
2429    INCLUDE "declare_STDlev.h"
2430!    INCLUDE "calcul_STDlev.h"
2431
2432!!!!!!! Variables
2433! pwrf_grid: WRF grid structure
2434
2435! L. Fita, LMD January 2014. It should work using suphel, but not?
2436!    CALL suphel
2437
2438    fname='"get_lmdz_out"'
2439    errmsg='ERROR -- error -- ERROR -- error'
2440
2441    RDAY = 86400.
2442    RG = 9.81
2443    RD = 287.
2444    RMO3=47.9942
2445
2446    PRINT *,'  Lluis in ' //TRIM(fname)//'!'
2447
2448! Remove variable assignation already done in get_lmdz? These variables are not used
2449!    in restart/boundaries
2450
2451    DO ix=1,ddx
2452      DO iy=1,ddy
2453        ixy=ddx*(iy-1)+ix
2454
2455        pwrf_grid%lages_ter(ix,iy) = agesno(ixy,1)
2456        pwrf_grid%lages_lic(ix,iy) = agesno(ixy,2)
2457        pwrf_grid%lages_oce(ix,iy) = agesno(ixy,3)
2458        pwrf_grid%lages_sic(ix,iy) = agesno(ixy,4)
2459!!        pwrf_grid%lahyb(misc??) = pahyb misc
2460        pwrf_grid%laire(ix,iy) = airephy(ixy)
2461        pwrf_grid%laireter(ix,iy) = paire_ter(ixy)
2462        pwrf_grid%lalb1(ix,iy) = albsol1(ixy)
2463        pwrf_grid%lalb2(ix,iy) = albsol2(ixy)
2464        pwrf_grid%lalbe_ter(ix,iy) = falb1(ixy,1)
2465        pwrf_grid%lalbe_lic(ix,iy) = falb1(ixy,2)
2466        pwrf_grid%lalbe_oce(ix,iy) = falb1(ixy,3)
2467        pwrf_grid%lalbe_sic(ix,iy) = falb1(ixy,4)
2468        pwrf_grid%lale(ix,iy) = ale(ixy)
2469        pwrf_grid%lale_bl(ix,iy) = ale_bl(ixy)
2470        pwrf_grid%lale_wk(ix,iy) = ale_wake(ixy)
2471        pwrf_grid%lalp(ix,iy) = alp(ixy)
2472        pwrf_grid%lalp_bl(ix,iy) = alp_bl(ixy)
2473        pwrf_grid%lalp_wk(ix,iy) = alp_wake(ixy)
2474!!        pwrf_grid%lalt(misc??) = palt misc
2475!!        pwrf_grid%lbhyb(misc??) = pbhyb misc
2476        pwrf_grid%lbils(ix,iy) = bils(ixy)
2477        pwrf_grid%lbils_diss(ix,iy) = bils_diss(ixy)
2478        pwrf_grid%lbils_ec(ix,iy) = bils_ec(ixy)
2479        pwrf_grid%lbils_enthalp(ix,iy) = bils_enthalp(ixy)
2480        pwrf_grid%lbils_kinetic(ix,iy) = bils_kinetic(ixy)
2481        pwrf_grid%lbils_latent(ix,iy) = bils_latent(ixy)
2482        pwrf_grid%lbils_tke(ix,iy) = bils_tke(ixy)
2483        pwrf_grid%lcape(ix,iy) = cape(ixy)
2484        pwrf_grid%lcape_max(ix,iy) = cape(ixy)
2485        pwrf_grid%lcdrh(ix,iy) = cdragh(ixy)
2486        pwrf_grid%lcdrm(ix,iy) = cdragm(ixy)
2487        pwrf_grid%lcldh(ix,iy) = cldh(ixy)
2488        pwrf_grid%lcldl(ix,iy) = cldl(ixy)
2489        pwrf_grid%lcldm(ix,iy) = cldm(ixy)
2490        pwrf_grid%lcldq(ix,iy) = cldq(ixy)
2491        pwrf_grid%lcldt(ix,iy) = cldt(ixy)
2492        pwrf_grid%lcontfracatm(ix,iy) = pctsrf(ixy,is_ter)+pctsrf(ixy,is_lic)
2493        pwrf_grid%lcontfracor(ix,iy) = pctsrf(ixy,is_ter)
2494!NOTknown!        pwrf_grid%lcumpb(ix,iy) = cumpb(ixy)
2495!NOTknown!        pwrf_grid%lcumrn(ix,iy) = cumrn(ixy)
2496        pwrf_grid%ldthmin(ix,iy) = dthmin(ixy)
2497        pwrf_grid%ldtsvdft(ix,iy) = d_ts(ixy,is_ter)
2498        pwrf_grid%ldtsvdfo(ix,iy) = d_ts(ixy,is_oce)
2499        pwrf_grid%ldtsvdfg(ix,iy) = d_ts(ixy,is_lic)
2500        pwrf_grid%ldtsvdfi(ix,iy) = d_ts(ixy,is_sic)
2501        pwrf_grid%levap(ix,iy) = fevap(ixy,1)
2502        pwrf_grid%levap_ter(ix,iy) = fevap(ixy,1)
2503        pwrf_grid%levap_lic(ix,iy) = fevap(ixy,2)
2504        pwrf_grid%levap_oce(ix,iy) = fevap(ixy,3)
2505        pwrf_grid%levap_sic(ix,iy) = fevap(ixy,4)
2506        pwrf_grid%levappot_ter(ix,iy) = evap_pot(ixy,1)
2507        pwrf_grid%levappot_lic(ix,iy) = evap_pot(ixy,2)
2508        pwrf_grid%levappot_oce(ix,iy) = evap_pot(ixy,3)
2509        pwrf_grid%levappot_sic(ix,iy) = evap_pot(ixy,4)
2510        pwrf_grid%lfbase(ix,iy) = ema_cbmf(ixy)
2511        pwrf_grid%lfder(ix,iy) = fder(ixy)
2512        pwrf_grid%lffonte(ix,iy) = zxffonte(ixy)
2513        pwrf_grid%lflat(ix,iy) = zxfluxlat(ixy)
2514        pwrf_grid%lflw_ter(ix,iy) = fsollw(ixy,1)
2515        pwrf_grid%lflw_lic(ix,iy) = fsollw(ixy,2)
2516        pwrf_grid%lflw_oce(ix,iy) = fsollw(ixy,3)
2517        pwrf_grid%lflw_sic(ix,iy) = fsollw(ixy,4)
2518        pwrf_grid%lfqcalving(ix,iy) = zxfqcalving(ixy)
2519        pwrf_grid%lfqfonte(ix,iy) = zxfqfonte(ixy)
2520        pwrf_grid%lfract_ter(ix,iy) = pctsrf(ixy,1)
2521        pwrf_grid%lfract_lic(ix,iy) = pctsrf(ixy,2)
2522        pwrf_grid%lfract_oce(ix,iy) = pctsrf(ixy,3)
2523        pwrf_grid%lfract_sic(ix,iy) = pctsrf(ixy,4)
2524        pwrf_grid%lfsnow(ix,iy) = zfra_o(ixy)
2525        pwrf_grid%lfsw_ter(ix,iy) = fsolw(ixy,1)
2526        pwrf_grid%lfsw_lic(ix,iy) = fsolw(ixy,2)
2527        pwrf_grid%lfsw_oce(ix,iy) = fsolw(ixy,3)
2528        pwrf_grid%lfsw_sic(ix,iy) = fsolw(ixy,4)
2529!NOtnow!        pwrf_grid%lftime_con(ix,iy) = float(itau_con)/float(itap)
2530        pwrf_grid%lftime_th(ix,iy) = 0.
2531        pwrf_grid%liwp(ix,iy) = fiwp(ixy)
2532!!        pwrf_grid%llat j = pllat j
2533        pwrf_grid%llat_ter(ix,iy) = fluxlat(ixy,1)
2534        pwrf_grid%llat_lic(ix,iy) = fluxlat(ixy,2)
2535        pwrf_grid%llat_oce(ix,iy) = fluxlat(ixy,3)
2536        pwrf_grid%llat_sic(ix,iy) = fluxlat(ixy,4)
2537        pwrf_grid%llmaxth(ix,iy) = lmax_th(ixy)
2538!!        pwrf_grid%llon i = pllon i
2539        pwrf_grid%llwdn200(ix,iy) = lwdn200(ixy)
2540        pwrf_grid%llwdn200clr(ix,iy) = lwdn200clr(ixy)
2541        pwrf_grid%llwdnsfc(ix,iy) = sollwdown(ixy)
2542        pwrf_grid%llwdnsfcclr(ix,iy) = sollwdownclr(ixy)
2543        pwrf_grid%llwdownor(ix,iy) = sollwdown(ixy)
2544        pwrf_grid%llwp(ix,iy) = flwp(ixy)
2545        pwrf_grid%llwup200(ix,iy) = lwup200(ixy)
2546        pwrf_grid%llwup200clr(ix,iy) = lwup200clr(ixy)
2547        pwrf_grid%llwupsfc(ix,iy) = sollwdown(ixy)-sollw(ixy)
2548        pwrf_grid%llwupsfcclr(ix,iy) = -1.*lwdn0(ixy,1)-sollw0(ixy)
2549        pwrf_grid%lmsnow(ix,iy) = snow_o(ixy)
2550        pwrf_grid%lndayrain(ix,iy) = nday_rain(ixy)
2551        pwrf_grid%lnettop(ix,iy) = topsw(ixy)-toplw(ixy)
2552        pwrf_grid%lpbase(ix,iy) = ema_pcb(ixy)
2553        pwrf_grid%lphis(ix,iy) = pphis(ixy)
2554        pwrf_grid%lplcl(ix,iy) = plcl(ixy)
2555        pwrf_grid%lplfc(ix,iy) = plfc(ixy)
2556        pwrf_grid%lpluc(ix,iy) = rain_con(ixy) + snow_con(ixy)
2557        pwrf_grid%lplul(ix,iy) = rain_lsc(ixy) + snow_lsc(ixy)
2558        pwrf_grid%lplulst(ix,iy) = plul_st(ixy)
2559        pwrf_grid%lplulth(ix,iy) = plul_th(ixy)
2560        pwrf_grid%lpourc_ter(ix,iy) = pctsrf(ixy,1)*100.
2561        pwrf_grid%lpourc_lic(ix,iy) = pctsrf(ixy,2)*100.
2562        pwrf_grid%lpourc_oce(ix,iy) = pctsrf(ixy,3)*100.
2563        pwrf_grid%lpourc_sic(ix,iy) = pctsrf(ixy,4)*100.
2564        pwrf_grid%lprecip(ix,iy) = rain_fall(ixy) + snow_fall(ixy)
2565!!        pwrf_grid%lpresnivs k = plpresnivs k
2566        pwrf_grid%lprw(ix,iy) = prw(ixy)
2567        pwrf_grid%lpsol(ix,iy) = paprs(ixy,1)
2568        pwrf_grid%lptop(ix,iy) = ema_pct(ixy)
2569        pwrf_grid%lq2m(ix,iy) = zq2m(ixy)
2570        pwrf_grid%lqsat2m(ix,iy) = qsat2m(ixy)
2571        pwrf_grid%lqsol(ix,iy) = qsol(ixy)
2572        pwrf_grid%lqsurf(ix,iy) = zxqsurf(ixy)
2573        pwrf_grid%lradsol(ix,iy) = radsol(ixy)
2574        pwrf_grid%lrh2m(ix,iy) = MIN(100.,rh2m(ixy)*100.)
2575        pwrf_grid%lrh2m_max(ix,iy) = MIN(100.,rh2m(ixy)*100.)
2576        pwrf_grid%lrh2m_min(ix,iy) = MIN(100.,rh2m(ixy)*100.)
2577        pwrf_grid%lrugs(ix,iy) = zxrugs(ixy)
2578        pwrf_grid%lrugs_ter(ix,iy) = frugs(ixy,1)
2579        pwrf_grid%lrugs_lic(ix,iy) = frugs(ixy,2)
2580        pwrf_grid%lrugs_oce(ix,iy) = frugs(ixy,3)
2581        pwrf_grid%lrugs_sic(ix,iy) = frugs(ixy,4)
2582        pwrf_grid%lsens(ix,iy) = -1.*sens(ixy)
2583        pwrf_grid%lsicf(ix,iy) = pctsrf(ixy,is_sic)
2584        pwrf_grid%ls_lcl(ix,iy) = s_lcl(ixy)
2585        pwrf_grid%lslp(ix,iy) = slp(ixy)
2586        pwrf_grid%lsnow(ix,iy) = snow_fall(ixy)
2587        pwrf_grid%lsnowl(ix,iy) = snow_lsc(ixy)
2588        pwrf_grid%lsoll(ix,iy) = sollw(ixy)
2589        pwrf_grid%lsoll0(ix,iy) = sollw0(ixy)
2590        pwrf_grid%lsollwdown(ix,iy) = sollwdown(ixy)
2591        pwrf_grid%lsols(ix,iy) = solsw(ixy)
2592        pwrf_grid%lsols0(ix,iy) = solsw0(ixy)
2593        pwrf_grid%ls_pblh(ix,iy) = s_pblh(ixy)
2594        pwrf_grid%ls_pblt(ix,iy) = s_pblt(ixy)
2595        pwrf_grid%ls_therm(ix,iy) = s_therm(ixy)
2596        pwrf_grid%lswdn200(ix,iy) = swdn200(ixy)
2597        pwrf_grid%lswdn200clr(ix,iy) = swdn200clr(ixy)
2598        pwrf_grid%lswdnsfc(ix,iy) = swdn(ixy,1)
2599        pwrf_grid%lswdnsfcclr(ix,iy) = swdn0(ixy,1)
2600        pwrf_grid%lswdntoa(ix,iy) = swdn(ixy,klev+1)
2601        pwrf_grid%lswdntoaclr(ix,iy) = swdn0(ixy,klev+1)
2602        pwrf_grid%lswdownor(ix,iy) = solsw(ixy)/(1.-albsol1(ixy))
2603        pwrf_grid%lswnetor(ix,iy) = fsolsw(ixy,is_ter)
2604        pwrf_grid%lswup200(ix,iy) = swup200(ixy)
2605        pwrf_grid%lswup200clr(ix,iy) = swup200clr(ixy)
2606        pwrf_grid%lswupsfc(ix,iy) = swup(ixy,1)
2607        pwrf_grid%lswupsfcclr(ix,iy) = swup0(ixy,1)
2608        pwrf_grid%lswuptoa(ix,iy) = swup(ixy,klev+1)
2609        pwrf_grid%lswuptoaclr(ix,iy) = swup0(ixy,klev+1)
2610        pwrf_grid%lt2m(ix,iy) = zt2m(ixy)
2611        pwrf_grid%lt2m_max(ix,iy) = zt2m(ixy)
2612        pwrf_grid%lt2m_min(ix,iy) = zt2m(ixy)
2613        pwrf_grid%lt2m_ter(ix,iy) = t2m(ixy,1)
2614        pwrf_grid%lt2m_lic(ix,iy) = t2m(ixy,2)
2615        pwrf_grid%lt2m_oce(ix,iy) = t2m(ixy,3)
2616        pwrf_grid%lt2m_sic(ix,iy) = t2m(ixy,4)
2617        pwrf_grid%ltaux(ix,iy) = 0.
2618        pwrf_grid%ltauy(ix,iy) = 0.
2619        DO iz=1,4
2620          pwrf_grid%ltaux(ix,iy) = pwrf_grid%ltauy(ix,iy) + pctsrf(ixy,iz)*          &
2621            fluxu(ixy,1,iz)
2622          pwrf_grid%ltauy(ix,iy) = pwrf_grid%ltauy(ix,iy) + pctsrf(ixy,iz)*          &
2623            fluxv(ixy,1,iz)
2624        END DO
2625        pwrf_grid%ltaux_ter(ix,iy) = fluxu(ixy,1,1)
2626        pwrf_grid%ltaux_lic(ix,iy) = fluxu(ixy,1,2)
2627        pwrf_grid%ltaux_oce(ix,iy) = fluxu(ixy,1,3)
2628        pwrf_grid%ltaux_sic(ix,iy) = fluxu(ixy,1,4)
2629        pwrf_grid%ltauy_ter(ix,iy) = fluxv(ixy,1,1)
2630        pwrf_grid%ltauy_lic(ix,iy) = fluxv(ixy,1,2)
2631        pwrf_grid%ltauy_oce(ix,iy) = fluxv(ixy,1,3)
2632        pwrf_grid%ltauy_sic(ix,iy) = fluxv(ixy,1,4)
2633!!        pwrf_grid%ltime - = pltime -
2634!!        pwrf_grid%ltime_bounds {lbnds} = pltime_bounds {lbnds}
2635!!        pwrf_grid%lti86400 - = plti86400 -
2636        IF ( (pctsrf(ixy,is_oce).GT.epsfra).OR.(pctsrf(ixy,is_sic).GT.epsfra) ) THEN
2637          pwrf_grid%lt_oce_sic(ix,iy) = (ftsol(ixy, is_oce) * pctsrf(ixy,is_oce)+    &
2638          ftsol(ixy, is_sic)*pctsrf(ixy,is_sic))/(pctsrf(ixy,is_oce)+                &
2639          pctsrf(ixy,is_sic))
2640        ELSE
2641          pwrf_grid%lt_oce_sic(ix,iy) = 273.15
2642        END IF
2643!!        pwrf_grid%lt_op_00001800(misc??) = pt_op_00001800 misc
2644!!        pwrf_grid%lt_op_00001800_bnds(misc??) = pt_op_00001800_bnds misc
2645        pwrf_grid%ltopl(ix,iy) = toplw(ixy)
2646        pwrf_grid%ltopl0(ix,iy) = toplw0(ixy)
2647        pwrf_grid%ltops(ix,iy) = topsw(ixy)
2648        pwrf_grid%ltops0(ix,iy) = topsw0(ixy)
2649        pwrf_grid%ltpot(ix,iy) = tpot(ixy)
2650        pwrf_grid%ltpote(ix,iy) = tpote(ixy)
2651        pwrf_grid%ltsol(ix,iy) = zxtsol(ixy)
2652        pwrf_grid%ltsol_ter(ix,iy) = ftsol(ixy,1)
2653        pwrf_grid%ltsol_lic(ix,iy) = ftsol(ixy,2)
2654        pwrf_grid%ltsol_oce(ix,iy) = ftsol(ixy,3)
2655        pwrf_grid%ltsol_sic(ix,iy) = ftsol(ixy,4)
2656        pwrf_grid%lu10m(ix,iy) = zu10m(ixy)
2657        pwrf_grid%lu10m_ter(ix,iy) = u10m(ixy,1)
2658        pwrf_grid%lu10m_lic(ix,iy) = u10m(ixy,2)
2659        pwrf_grid%lu10m_oce(ix,iy) = u10m(ixy,3)
2660        pwrf_grid%lu10m_sic(ix,iy) = u10m(ixy,4)
2661        pwrf_grid%lue(ix,iy) = ue(ixy)
2662        pwrf_grid%luq(ix,iy) = uq(ixy)
2663        pwrf_grid%lustar(ix,iy) = zustar(ixy)
2664        pwrf_grid%lustar_ter(ix,iy) = ustar(ixy,1)
2665        pwrf_grid%lustar_lic(ix,iy) = ustar(ixy,2)
2666        pwrf_grid%lustar_oce(ix,iy) = ustar(ixy,3)
2667        pwrf_grid%lustar_sic(ix,iy) = ustar(ixy,4)
2668        pwrf_grid%lv10m(ix,iy) = zv10m(ixy)
2669        pwrf_grid%lv10m_ter(ix,iy) = v10m(ixy,1)
2670        pwrf_grid%lv10m_lic(ix,iy) = v10m(ixy,2)
2671        pwrf_grid%lv10m_oce(ix,iy) = v10m(ixy,3)
2672        pwrf_grid%lv10m_sic(ix,iy) = v10m(ixy,4)
2673        pwrf_grid%lve(ix,iy) = ve(ixy)
2674        pwrf_grid%lvq(ix,iy) = vq(ixy)
2675        pwrf_grid%lwake_h(ix,iy) = wake_h(ixy)
2676        pwrf_grid%lwake_s(ix,iy) = wake_s(ixy)
2677        pwrf_grid%lwape(ix,iy) = wake_pe(ixy)
2678        pwrf_grid%lwbeff(ix,iy) = wbeff(ixy)
2679        pwrf_grid%lwbilo_ter(ix,iy) = wfbilo(ixy,1)
2680        pwrf_grid%lwbilo_lic(ix,iy) = wfbilo(ixy,2)
2681        pwrf_grid%lwbilo_oce(ix,iy) = wfbilo(ixy,3)
2682        pwrf_grid%lwbilo_sic(ix,iy) = wfbilo(ixy,4)
2683        pwrf_grid%lwbils_ter(ix,iy) = wfbils(ixy,1)
2684        pwrf_grid%lwbils_lic(ix,iy) = wfbils(ixy,2)
2685        pwrf_grid%lwbils_oce(ix,iy) = wfbils(ixy,3)
2686        pwrf_grid%lwbils_sic(ix,iy) = wfbils(ixy,4)
2687        pwrf_grid%lweakinv(ix,iy) = weak_inversion(ixy)
2688        pwrf_grid%lwind10m(ix,iy) = SQRT(zu10m(ixy)*zu10m(ixy)+zv10m(ixy)*zv10m(ixy))
2689        pwrf_grid%lwind10max(ix,iy) = SQRT(zu10m(ixy)*zu10m(ixy)+zv10m(ixy)*zv10m(ixy))
2690        pwrf_grid%lzmax_th(ix,iy) = zmax_th(ixy)
2691! Output at the standard levels
2692!!
2693!        DO iz=1, nlevSTD
2694!          bb2=clevSTD(iz)
2695 
2696!          IF (TRIM(bb2) == "850") THEN
2697!            pwrf_grid%lq850(ix,iy) = q850(ixy)
2698!            pwrf_grid%lt850(ix,iy) = t850(ixy)
2699!            pwrf_grid%lu850(ix,iy) = u850(ixy)
2700!            pwrf_grid%lv850(ix,iy) = v850(ixy)
2701!            pwrf_grid%lw850(ix,iy) = w850(ixy)
2702!            pwrf_grid%lz850(ix,iy) = z850(ixy)
2703!          ELSE IF (TRIM(bb2) == "700") THEN
2704!            pwrf_grid%lq700(ix,iy) = q700(ixy)
2705!            pwrf_grid%lt700(ix,iy) = t700(ixy)
2706!            pwrf_grid%lu700(ix,iy) = u700(ixy)
2707!            pwrf_grid%lv700(ix,iy) = v700(ixy)
2708!            pwrf_grid%lw700(ix,iy) = w700(ixy)
2709!            pwrf_grid%lz700(ix,iy) = z700(ixy)
2710!          ELSE IF (TRIM(bb2) == "500") THEN
2711!            pwrf_grid%lq500(ix,iy) = q500(ixy)
2712!            pwrf_grid%lt500(ix,iy) = t500(ixy)
2713!            pwrf_grid%lu500(ix,iy) = u500(ixy)
2714!            pwrf_grid%lv500(ix,iy) = v500(ixy)
2715!            pwrf_grid%lw500(ix,iy) = w500(ixy)
2716!            pwrf_grid%lz500(ix,iy) = z500(ixy)
2717!           ELSE IF (TRIM(bb2) == "200")THEN
2718!            pwrf_grid%lq200(ix,iy) = q200(ixy)
2719!            pwrf_grid%lt200(ix,iy) = t200(ixy)
2720!            pwrf_grid%lu200(ix,iy) = u200(ixy)
2721!            pwrf_grid%lv200(ix,iy) = v200(ixy)
2722!            pwrf_grid%lw200(ix,iy) = w200(ixy)
2723!            pwrf_grid%lz200(ix,iy) = z200(ixy)
2724!          ELSE IF (TRIM(bb2) == "100") THEN
2725!            pwrf_grid%lq100(ix,iy) = q100(ixy)
2726!            pwrf_grid%lt100(ix,iy) = t100(ixy)
2727!            pwrf_grid%lu100(ix,iy) = u100(ixy)
2728!            pwrf_grid%lv100(ix,iy) = v100(ixy)
2729!            pwrf_grid%lw100(ix,iy) = w100(ixy)
2730!            pwrf_grid%lz100(ix,iy) = z100(ixy)
2731!          ELSE IF (TRIM(bb2) == "50") THEN
2732!            pwrf_grid%lq50(ix,iy) = q50(ixy)
2733!            pwrf_grid%lt50(ix,iy) = t50(ixy)
2734!            pwrf_grid%lu50(ix,iy) = u50(ixy)
2735!            pwrf_grid%lv50(ix,iy) = v50(ixy)
2736!            pwrf_grid%lw50(ix,iy) = w50(ixy)
2737!            pwrf_grid%lz50(ix,iy) = z50(ixy)
2738!          ELSE IF (TRIM(bb2) == "10") THEN
2739!            pwrf_grid%lq10(ix,iy) = q10(ixy)
2740!            pwrf_grid%lt10(ix,iy) = t10(ixy)
2741!            pwrf_grid%lu10(ix,iy) = u10(ixy)
2742!            pwrf_grid%lv10(ix,iy) = v10(ixy)
2743!            pwrf_grid%lw10(ix,iy) = w10(ixy)
2744!            pwrf_grid%lz10(ix,iy) = z10(ixy)
2745!          END IF
2746!        END DO
2747
2748      END DO
2749    END DO
2750
2751    dx=ddx+2*bdy
2752    dy=ddy+2*bdy
2753
2754    CALL vect_mat_zstagg(fraca, dx, dy, ddz+1, bdy, pwrf_grid%la_th)
2755    CALL vect_mat(beta_prec, dx, dy, ddz, bdy, pwrf_grid%lbeta_prec)
2756    var3d=upwd + dnwd+ dnwd0
2757    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldmc)
2758    CALL vect_mat(dnwd, dx, dy, ddz, bdy, pwrf_grid%ldnwd)
2759    CALL vect_mat(dnwd0, dx, dy, ddz, bdy, pwrf_grid%ldnwd0)
2760! Originally is var3d = d_q_ajsb(1:ddx*ddy,1:ddz)/pdtphys
2761    var3d = d_q_ajsb/pdtphys
2762    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqajs)
2763    var3d = d_q_con/pdtphys
2764    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqcon)
2765    CALL vect_mat(d_q_dyn, dx, dy, ddz, bdy, pwrf_grid%ldqdyn)
2766    var3d = d_q_eva/pdtphys
2767    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqeva)
2768    var3d = d_q_lsc/pdtphys
2769    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlsc)
2770    var3d = d_q_lscst/pdtphys
2771    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlscst)
2772    var3d = d_q_lscth/pdtphys
2773    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlscth)
2774! L. Fita, LMD January 2014. d_qx is the tendency of moisture due to physiq!
2775!   MOISTtend in the WRF_LMDZ coupling
2776    CALL vect_mat(d_qx(:,:,ivap), dx, dy, ddz, bdy, pwrf_grid%ldqphy)
2777    var3d = d_q_ajs/pdtphys - d_q_ajsb/pdtphys
2778    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqthe)
2779    var3d = d_q_vdf/pdtphys
2780    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqvdf)
2781    var3d = d_q_wake/pdtphys
2782    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqwak)
2783    var3d = d_t_ajsb/pdtphys
2784    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtajs)
2785    var3d = d_t_con/pdtphys
2786    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtcon)
2787    var3d = d_t_diss/pdtphys
2788    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtdis)
2789    CALL vect_mat(d_t_dyn, dx, dy, ddz, bdy, pwrf_grid%ldtdyn)
2790    var3d = d_t_ec/pdtphys
2791    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtec)
2792    var3d = d_t_eva/pdtphys
2793    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldteva)
2794    CALL vect_mat(detr_therm, dx, dy, ddz, bdy, pwrf_grid%ld_th)
2795    CALL vect_mat(cldemi, dx, dy, ddz, bdy, pwrf_grid%lcldemi)
2796    CALL vect_mat(cldtau, dx, dy, ddz, bdy, pwrf_grid%lcldtau)
2797    CALL vect_mat(clwcon, dx, dy, ddz, bdy, pwrf_grid%lclwcon)
2798    var3d = d_t_lif/pdtphys
2799    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlif)
2800    var3d = d_t_lsc/pdtphys
2801    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlsc)
2802    var3d = (d_t_lsc+d_t_eva)/pdtphys
2803    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlschr)
2804    var3d = d_t_lscst/pdtphys
2805    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlscst)
2806    var3d = d_t_lscth/pdtphys
2807    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlscth)
2808    var3d=-1.*cool0/RDAY
2809    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlw0)
2810    var3d = -1.*cool/RDAY
2811    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlwr)
2812    var3d=d_t_oro/pdtphys
2813    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtoro)
2814    CALL vect_mat(d_t, dx, dy, ddz, bdy, pwrf_grid%ldtphy)
2815    var3d=heat0/RDAY
2816    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtsw0)
2817    var3d=heat/RDAY
2818    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtswr)
2819    var3d = d_t_ajs/pdtphys - d_t_ajsb/pdtphys
2820    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtthe)
2821    var3d=d_t_vdf/pdtphys
2822    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtvdf)
2823    var3d=d_t_wake/pdtphys
2824    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtwak)
2825    var3d=d_u_con/pdtphys
2826    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lducon)
2827    CALL vect_mat(d_u_dyn, dx, dy, ddz, bdy, pwrf_grid%ldudyn)
2828    var3d=d_u_lif/pdtphys
2829    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldulif)
2830    var3d=d_u_oro/pdtphys
2831    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lduoro)
2832    var3d=d_u_vdf/pdtphys
2833    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lduvdf)
2834    var3d=d_v_con/pdtphys
2835    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvcon)
2836    CALL vect_mat(d_v_dyn, dx, dy, ddz, bdy, pwrf_grid%ldvdyn)
2837    var3d=d_v_lif/pdtphys
2838    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvlif)
2839    var3d=d_v_oro/pdtphys
2840    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvoro)
2841    var3d=d_v_vdf/pdtphys
2842    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvvdf)
2843!NOTfound!    CALL vect_mat(ec550aer, dx, dy, ddz, bdy, pwrf_grid%lec550aer)
2844    CALL vect_mat(entr_therm, dx, dy, ddz, bdy, pwrf_grid%le_th)
2845    CALL vect_mat(coefm(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%levu)
2846    CALL vect_mat(fl, dx, dy, ddz, bdy, pwrf_grid%lfl)
2847    CALL vect_mat(zphi, dx, dy, ddz, bdy, pwrf_grid%lgeop)
2848    var3d = q_seri + ql_seri
2849    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lh2o)
2850    CALL vect_mat(fiwc, dx, dy, ddz, bdy, pwrf_grid%liwcon)
2851    CALL vect_mat(coefh(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%lkz)
2852    CALL vect_mat(coefh(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%lkz_max)
2853!    CALL vect_mat(lambda_th, dx, dy, ddz, bdy, pwrf_grid%llambda_th)
2854    CALL vect_mat(flwc, dx, dy, ddz, bdy, pwrf_grid%llwcon)
2855!NOTknown!    CALL vect_mat(ma, dx, dy, ddz, bdy, pwrf_grid%lma)
2856    CALL vect_mat(zmasse, dx, dy, ddz, bdy, pwrf_grid%lmass)
2857!NOTknown!    CALL vect_mat(mc, dx, dy, ddz, bdy, pwrf_grid%lmc)
2858!NOTknown!    CALL vect_mat(mcd, dx, dy, ddz, bdy, pwrf_grid%lmcd)
2859    CALL vect_mat(ql_seri, dx, dy, ddz, bdy, pwrf_grid%loliq)
2860    CALL vect_mat(q_seri, dx, dy, ddz, bdy, pwrf_grid%lovap)
2861! L. Fita, LMD January 2014. d_qx is the tendency of moisture due to physiq!
2862!   MOISTtend in the WRF_LMDZ coupling
2863    CALL vect_mat(qx(:,:,ivap), dx, dy, ddz, bdy, pwrf_grid%lovapinit)
2864!Notnow!    var3d = wo(:,:,1)*dobson_u*1e3 / zmasse / rmo3 * rmd)
2865!Notnow!    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lozone)
2866    CALL vect_mat_zstagg(paprs, dx, dy, ddz+1, bdy, pwrf_grid%lpaprs)
2867!NOTknown!    CALL vect_mat(pb, dx, dy, ddz, bdy, pwrf_grid%lpb)
2868    CALL vect_mat_zstagg(pmflxs, dx, dy, ddz+1, bdy, pwrf_grid%lpr_con_i)
2869    CALL vect_mat_zstagg(pmflxr, dx, dy, ddz+1, bdy, pwrf_grid%lpr_con_l)
2870    CALL vect_mat(pplay, dx, dy, ddz, bdy, pwrf_grid%lpres)
2871    CALL vect_mat_zstagg(psfl, dx, dy, ddz+1, bdy, pwrf_grid%lpr_lsc_i)
2872    CALL vect_mat_zstagg(prfl, dx, dy, ddz+1, bdy, pwrf_grid%lpr_lsc_l)
2873    var3d = 0.
2874    WHERE (ptconv) var3d = 1.
2875    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lptconv)
2876    CALL vect_mat(zqasc, dx, dy, ddz, bdy, pwrf_grid%lq_th)
2877    CALL vect_mat(ratqs, dx, dy, ddz, bdy, pwrf_grid%lratqs)
2878    CALL vect_mat(re, dx, dy, ddz, bdy, pwrf_grid%lre)
2879    CALL vect_mat(ref_ice, dx, dy, ddz, bdy, pwrf_grid%lref_ice)
2880    CALL vect_mat(ref_liq, dx, dy, ddz, bdy, pwrf_grid%lref_liq)
2881    CALL vect_mat(zx_rh, dx, dy, ddz, bdy, pwrf_grid%lrhum)
2882    CALL vect_mat(lwdn, dx, dy, ddz, bdy, pwrf_grid%lrld)
2883    CALL vect_mat(lwdn0, dx, dy, ddz, bdy, pwrf_grid%lrldcs)
2884    CALL vect_mat(lwup, dx, dy, ddz, bdy, pwrf_grid%lrlu)
2885    CALL vect_mat(lwup0, dx, dy, ddz, bdy, pwrf_grid%lrlucs)
2886!NOTknown!    CALL vect_mat(rn, dx, dy, ddz, bdy, pwrf_grid%lrn)
2887    CALL vect_mat(cldfra, dx, dy, ddz, bdy, pwrf_grid%lrneb)
2888    CALL vect_mat(rnebcon, dx, dy, ddz, bdy, pwrf_grid%lrnebcon)
2889    CALL vect_mat(rneb, dx, dy, ddz, bdy, pwrf_grid%lrnebls)
2890    CALL vect_mat(swdn, dx, dy, ddz, bdy, pwrf_grid%lrsd)
2891    CALL vect_mat(swdn0, dx, dy, ddz, bdy, pwrf_grid%lrsdcs)
2892    CALL vect_mat(swup, dx, dy, ddz, bdy, pwrf_grid%lrsu)
2893    CALL vect_mat(swup0, dx, dy, ddz, bdy, pwrf_grid%lrsucs)
2894    CALL vect_mat(fluxt(:,:,1), dx, dy, ddz, bdy, pwrf_grid%lsens_ter)
2895    CALL vect_mat(fluxt(:,:,2), dx, dy, ddz, bdy, pwrf_grid%lsens_lic)
2896    CALL vect_mat(fluxt(:,:,3), dx, dy, ddz, bdy, pwrf_grid%lsens_oce)
2897    CALL vect_mat(fluxt(:,:,4), dx, dy, ddz, bdy, pwrf_grid%lsens_sic)
2898    CALL vect_mat(t_seri, dx, dy, ddz, bdy, pwrf_grid%ltemp)
2899    CALL vect_mat(theta, dx, dy, ddz, bdy, pwrf_grid%ltheta)
2900    var3d=0.
2901    DO nsrf=1,nbsrf
2902       DO iz=1,ddz+1
2903         var3dstagg(:,iz)=var3dstagg(:,iz)+pctsrf(:,nsrf)*pbl_tke(:,iz,nsrf)
2904       ENDDO
2905    ENDDO
2906    CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%ltke)
2907    CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%ltke_max)
2908    CALL vect_mat(pbl_tke(:,:,1), dx, dy, ddz+1, bdy, pwrf_grid%ltke_ter)
2909    CALL vect_mat(pbl_tke(:,:,2), dx, dy, ddz+1, bdy, pwrf_grid%ltke_lic)
2910    CALL vect_mat(pbl_tke(:,:,3), dx, dy, ddz+1, bdy, pwrf_grid%ltke_oce)
2911    CALL vect_mat(pbl_tke(:,:,4), dx, dy, ddz+1, bdy, pwrf_grid%ltke_sic)
2912!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,1), dx, dy, ddz+1, bdy, pwrf_grid%ltke_ter_max)
2913!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,2), dx, dy, ddz+1, bdy, pwrf_grid%ltke_lic_max)
2914!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,3), dx, dy, ddz+1, bdy, pwrf_grid%ltke_oce_max)
2915!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,4), dx, dy, ddz+1, bdy, pwrf_grid%ltke_sic_max)
2916    var3d = d_qx(:,:,ivap)+d_q_dyn
2917    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhus)
2918    IF (iflag_thermals == 1) THEN
2919      var3d = d_q_con/pdtphys
2920    ELSE IF ((iflag_thermals > 1) .AND. (iflag_wake == 1)) THEN
2921      var3d = d_q_con/pdtphys + d_q_ajs/pdtphys + d_q_wake/pdtphys
2922    END IF
2923    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhusc)
2924    var3d = d_q_lsc/pdtphys + d_q_eva/pdtphys
2925    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhusscpbl)
2926    var3d = d_t + d_t_dyn
2927    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnt)
2928    IF (iflag_thermals == 1) THEN
2929      var3d = d_t_con/pdtphys + d_t_ajsb/pdtphys
2930    ELSE IF ((iflag_thermals  > 1) .AND. (iflag_wake == 1)) THEN
2931      var3d=d_t_con/pdtphys + d_t_ajs/pdtphys + d_t_wake/pdtphys
2932    END IF
2933    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntc)
2934    var3d = heat/RDAY - cool/RDAY
2935    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntr)
2936    var3d = (d_t_lsc + d_t_eva + d_t_vdf)/pdtphys
2937    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntscpbl)
2938    CALL vect_mat(upwd, dx, dy, ddz, bdy, pwrf_grid%lupwd)
2939    CALL vect_mat(u_seri, dx, dy, ddz, bdy, pwrf_grid%lvitu)
2940    CALL vect_mat(v_seri, dx, dy, ddz, bdy, pwrf_grid%lvitv)
2941    CALL vect_mat(omega, dx, dy, ddz, bdy, pwrf_grid%lvitw)
2942    CALL vect_mat_zstagg(Vprecip, dx, dy, ddz+1, bdy, pwrf_grid%lvprecip)
2943    CALL vect_mat(wake_deltaq, dx, dy, ddz, bdy, pwrf_grid%lwake_deltaq)
2944    CALL vect_mat(wake_deltat, dx, dy, ddz, bdy, pwrf_grid%lwake_deltat)
2945    CALL vect_mat(wake_omg, dx, dy, ddz, bdy, pwrf_grid%lwake_omg)
2946    CALL vect_mat(wdtrainA, dx, dy, ddz, bdy, pwrf_grid%lwdtraina)
2947    CALL vect_mat(wdtrainM, dx, dy, ddz, bdy, pwrf_grid%lwdtrainm)
2948! L. Fita, LMD, January 2014. pphis is the given value of the pressure at the surface
2949    var3dstagg(:,1) = pphis(:)/RG
2950    DO iz=1,ddz
2951      var3dstagg(:,iz+1)= var3dstagg(:,iz)-(t_seri(:,iz)*RD*(paprs(:,iz+1)-          &
2952       paprs(:,iz)) )/ (pplay(:,iz)*RG)
2953    ENDDO
2954    CALL vect_mat(var3dstagg, dx, dy, ddz, bdy, pwrf_grid%lzfull)
2955    var3dstagg(:,1) = pphis(:)/RG-( (t_seri(:,1)+zxtsol)/2.*RD*(pplay(:,1)-          &
2956      paprs(:,1)) )/ ((paprs(:,1)+pplay(:,1))/2.*RG)
2957    DO iz=1, ddz-1
2958      var3dstagg(:,iz+1)= var3dstagg(:,iz)-( (t_seri(:,iz)+t_seri(:,iz+1))/          &
2959        2.*RD*(pplay(:,iz+1)-pplay(:,iz)))/( paprs(:,iz)*RG)
2960    ENDDO
2961    CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%lzhalf)
2962
2963    RETURN
2964
2965  END SUBROUTINE get_lmdz_out
2966
2967  SUBROUTINE put_lmdz_in_WRFout(is, ie, js, je, ks, ke, ddx, ddy, ddz, bdy, dlmdz,xp,&
2968    yp, lmdp, Nks, Nz0, nqtot, ivap, iliq, iflag_pbl, iflag_con, iflag_thermals,     &
2969    iflag_wake, itap, pdtphys, paprs, pplay, pphi, pphis, presnivs, qx, d_u, d_v,    &
2970    d_t, d_qx, d_ps, orchidlev, tsoilorchid, nslayers, wslay, wzs, wdzs, wnest_pos,  &
2971    wq2, wt2, wth2, wpsfc, wu10, wv10, wresm, wzetatop, wtslb, wsmois,               &
2972    wsh2o, wsmcrel, wsfroff, wudroff, wgrdflx, wacgrdflx, wsnow, wsnowh,             &
2973    wcanwat, wsst, wsstsk, wlai, wh_diabatic, wsina, wtsk, wtiso,                    &
2974    wmax_msftx, wmax_msfty, wrainc, wraincv, wrainsh, wrainshv, wrainnc, wrainncv,   &
2975    wsnownc, wsnowncv, wgraupelnc, wgraupelncv, whailnc, whailncv, wstepave_count,   &
2976    wcldfra, wswdown, wglw, wswnorm, wolr, wemiss, wtmn, wust, wpblh, whfx,          &
2977    wqfx, wlh, wachfx, waclhf, wsnowc, wsr, wpotevp, wsnopcx, wsoiltb)
2978!    wsave_topo_from_real, wavg_fuel_frac, wuah, wvah, wseed1, wseed2)
2979! Subroutine to get LMDZ output and get into WRF standard output variables
2980!   (as they appear in the wrfout)
2981
2982! LMDZ modules
2983    USE indice_sol_mod
2984    USE phys_state_var_mod
2985    USE phys_output_var_mod
2986!    USE pbl_surface_mod
2987    USE phys_local_var_mod
2988    USE comgeomphy
2989
2990! WRF_LMDZ modules
2991    USE output_lmdz_NOmodule
2992    USE lmdz_wrf_variables_mod
2993
2994    IMPLICIT NONE
2995
2996    INTEGER, INTENT(IN)                                  :: is, ie, js, je, ks, ke,  &
2997      ddx, ddy, ddz, bdy, dlmdz, xp, yp, lmdp, Nks, Nz0, nqtot, ivap, iliq,          &
2998      iflag_pbl, iflag_con, iflag_thermals, iflag_wake, itap, nslayers
2999    REAL, INTENT(IN)                                     :: pdtphys
3000    REAL, DIMENSION(dlmdz), INTENT(IN)                   :: pphis, d_ps
3001    REAL, DIMENSION(ddz), INTENT(IN)                     :: presnivs
3002    REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN)             :: paprs
3003    REAL, DIMENSION(dlmdz,ddz), INTENT(IN)               :: pplay, pphi, d_u, d_v, d_t
3004    REAL, DIMENSION(dlmdz,ddz,nqtot), INTENT(IN)         :: qx, d_qx
3005    REAL, DIMENSION(Nz0), INTENT(IN)                     :: orchidlev
3006    REAL, DIMENSION(dlmdz,Nz0,Nks), INTENT(IN)           :: tsoilorchid
3007    INTEGER, INTENT(INOUT)                               :: wstepave_count
3008    REAL, INTENT(INOUT)                                  :: wresm, wzetatop, wtiso,  &
3009      wmax_msftx, wmax_msfty
3010    REAL, DIMENSION(1:nslayers), INTENT(INOUT)           :: wzs, wdzs
3011    REAL, DIMENSION(1:nslayers), INTENT(IN)              :: wslay
3012    REAL, DIMENSION(is:ie,js:je), INTENT(INOUT)          :: wnest_pos, wq2, wt2,     &
3013      wth2, wpsfc, wu10, wv10, wsfroff, wudroff, wgrdflx, wacgrdflx, wsnow, wsnowh,  &
3014      wcanwat, wsst, wsstsk, wlai, wsina, wtsk, wrainc, wraincv, wrainsh, wrainshv,  &
3015      wrainnc, wrainncv, wsnownc, wsnowncv, wgraupelnc, wgraupelncv, whailnc,        &
3016      whailncv, wswdown, wglw, wswnorm, wolr, wemiss, wtmn, wust, wpblh, whfx,       &
3017      wqfx, wlh, wachfx, waclhf, wsnowc, wsr, wpotevp, wsnopcx, wsoiltb
3018    REAL, DIMENSION(is:ie,1:nslayers,js:je), INTENT(INOUT) :: wtslb, wsmois, wsh2o,  &
3019      wsmcrel
3020    REAL, DIMENSION(is:ie,ks:ke,js:je), INTENT(INOUT)    :: wh_diabatic, wcldfra
3021!    EXTERNAL suphel
3022
3023! Local
3024    INTEGER                                              :: dx, dy, ix, iy, ixy, iz
3025    INTEGER                                              :: nsrf
3026    REAL                                                 :: varixy
3027    REAL, DIMENSION(Nz0)                                 :: yvals1
3028    REAL                                                 :: RDAY, RG, RD
3029    CHARACTER(LEN=50)                                    :: fname
3030
3031!!!!!!! Variables
3032! pwrf_grid: WRF grid structure
3033
3034    fname = "'put_lmdz_in_WRFout'"
3035
3036! L. Fita, LMD January 2014. It should work using suphel, but not?
3037!    CALL suphel
3038
3039    RDAY = 86400.
3040    RG = 9.81
3041    RD = 287.
3042
3043    wdzs(1)=wslay(nslayers)
3044    DO iz=1,nslayers-1
3045      wdzs(iz+1) = wslay(nslayers-iz+1) - wslay(nslayers-iz)
3046    END DO
3047!    wmax_msftx =
3048!    wmax_msfty =
3049!    wresm =
3050!    wstepave_count =
3051!    wtiso =
3052!    wzetatop =
3053    wzs = wslay
3054
3055    DO ix=1,ddx
3056      DO iy=1,ddy
3057        ixy=ddx*(iy-1)+ix
3058        wacgrdflx(ix,iy) = wacgrdflx(ix,iy) + 0.
3059        wachfx(ix,iy) = wachfx(ix,iy) + sens(ixy)
3060        waclhf(ix,iy) = waclhf(ix,iy) + zxfluxlat(ixy)
3061        DO iz=1,ddz
3062          wcldfra(ix,iz,iy) = cldfra(ixy,iz)
3063        END DO
3064        wcanwat(ix,iy) = 0.
3065        wglw(ix,iy) = sollwdown(ixy)
3066! These are precipitation fluxes !!
3067!        wgraupelnc(ix,iy) = wgraupelnc(ix,iy) + prmflxs(ixy) + prfs(ixy,iz)
3068!        wgraupelncv(ix,iy) = prmflxs(ixy) + prfs(ixy,iz)
3069        wgrdflx(ix,iy) = 0.
3070        whailnc(ix,iy) = whailnc(ix,iy) + 0.
3071        whailncv(ix,iy) = 0.
3072        wh_diabatic(ix,ks:ke,iy) = 0.
3073        whfx(ix,iy) = sens(ixy)
3074        wlai(ix,iy) = 0.
3075        wlh(ix,iy) = zxfluxlat(ixy)
3076        wnest_pos(ix,iy) = 0.
3077!!        wnoahres(ix,iy) = 0.
3078        wolr(ix,iy) = toplw(ixy)
3079        wpblh(ix,iy) = s_pblh(ixy)
3080! Should be something from evappot_srf ?
3081        wpotevp(ix,iy) = wpotevp(ix,iy) + 0.
3082        wpsfc(ix,iy) = paprs(ixy,1)
3083        wq2(ix,iy) = zq2m(ixy)
3084        wqfx(ix,iy) = 0.
3085        wrainc(ix,iy) = wrainc(ix,iy) + rain_fall(ixy)*pdtphys
3086        wraincv(ix,iy) = rain_fall(ixy)*pdtphys
3087        wrainnc(ix,iy) = wrainnc(ix,iy) + rain_lsc(ixy)*pdtphys
3088        wrainncv(ix,iy) = rain_lsc(ixy)*pdtphys
3089        wrainsh(ix,iy) = wrainsh(ix,iy) + 0.
3090        wrainshv(ix,iy) = 0.
3091! L. Fita, LMD. February 2014
3092! LMDZ run_off_lic is for continental ice. We could try to get run_off_ter, but....
3093!        wsfroff(ix,iy) = 0.
3094        DO iz=1,Nks
3095          yvals1(iz) = tsoilorchid(ixy,iz,1)*pctsrf(ixy,is_ter) +                    &
3096            tsoilorchid(ixy,iz,2)*pctsrf(ixy,is_lic) +                               &
3097            tsoilorchid(ixy,iz,3)*pctsrf(ixy,is_oce) +                               &
3098            tsoilorchid(ixy,iz,4)*pctsrf(ixy,is_sic)
3099        END DO
3100        CALL interpolate_layers(Nz0,orchidlev,yvals1,nslayers,wslay,wtslb(ix,:,iy))
3101! No way ???
3102!          wsh2o(is:ie,iz,js:je) = wsh2o(is:ie,iz,js:je) +
3103!          wsmcrel(is:ie,iz,js:je) = 0.
3104
3105!        wsinalpha(ix,iy) =
3106        wsina(ix,iy) = 0.
3107        wsnopcx(ix,iy) = wsnopcx(ix,iy) + 0.
3108!!        wsnowh(ix,iy) = snowght(ixy)
3109        wsnownc(ix,iy) = wsnownc(ix,iy) + snow_fall(ixy)*pdtphys
3110        wsnowncv(ix,iy) = snow_fall(ixy)*pdtphys
3111        wsoiltb(ix,iy) = wtslb(ix,nslayers,iy)
3112        IF (rain_fall(ixy) /= 0.) THEN
3113          wsr(ix,iy) = ( snow_fall(ixy) ) / ( rain_fall(ixy))
3114        ELSE
3115          wsr(ix,iy) = 0.
3116        END IF
3117        wsstsk(ix,iy) = ftsol(ixy,is_oce)
3118        wsst(ix,iy) = sst(ixy)
3119        wswdown(ix,iy) = swdn(ixy,1)*pctsrf(ixy,is_ter) +                            &
3120          swdn(ixy,2)*pctsrf(ixy,is_lic) + swdn(ixy,3)*pctsrf(ixy,is_oce) +          &
3121          swdn(ixy,4)*pctsrf(ixy,is_sic)
3122        wswnorm(ix,iy) = 0.
3123        wt2(ix,iy) = zt2m(ixy)
3124! This is not exact!!
3125        CALL temp_to_theta1(zt2m(ixy), paprs(ixy,1), wth2(ix,iy))
3126        wtsk(ix,iy) = ftsol(ixy, is_ter)*pctsrf(ixy,is_ter) + ftsol(ixy, is_lic)*    &
3127          pctsrf(ixy,is_lic) + ftsol(ixy, is_oce)*pctsrf(ixy,is_oce) +               &
3128          ftsol(ixy, is_sic)*pctsrf(ixy,is_sic)
3129        wu10(ix,iy) = zu10m(ixy)
3130        wudroff(ix,iy) = wudroff(ix,iy) + 0.
3131        wv10(ix,iy) = zv10m(ixy)
3132      END DO
3133    END DO
3134
3135  END SUBROUTINE put_lmdz_in_WRFout
3136
3137END MODULE wrf_lmdz_mod
3138
Note: See TracBrowser for help on using the repository browser.