source: lmdz_wrf/trunk/WRFV3/lmdz/wrf_lmdz_mod.F90 @ 1939

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

Removing checking printings
File size: 119.4 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    DO iy=1,ddy
932      DO ix=1,ddx
933        ixy=ddx*(iy-1)+ix
934! It does not change on time: pzmask
935! It might change due to ice dynamics?
936        pwrf_grid%lter(ix,iy) = ppctsrf(ixy,is_ter)
937        pwrf_grid%llic(ix,iy) = ppctsrf(ixy,is_lic)
938        pwrf_grid%loce(ix,iy) = ppctsrf(ixy,is_oce)
939        pwrf_grid%lsic(ix,iy) = ppctsrf(ixy,is_sic)
940        DO iz=1,Nks
941          pwrf_grid%ltksoil(ix,iz,iy) = pftsol(ixy,iz)
942        END DO
943        DO iz=1,NzO
944          pwrf_grid%lotter(ix,iz,iy) = pftsoil(ixy,iz,1)
945          pwrf_grid%lotlic(ix,iz,iy) = pftsoil(ixy,iz,2)
946          pwrf_grid%lotoce(ix,iz,iy) = pftsoil(ixy,iz,3)
947          pwrf_grid%lotsic(ix,iz,iy) = pftsoil(ixy,iz,4)
948        END DO
949        DO iz=1,Nks
950          pwrf_grid%lqksoil(ix,iz,iy) = pqsurf(ixy,iz)
951        END DO
952        pwrf_grid%lwsol(ix,iy) = pqsol(ixy)
953        DO iz=1,Nks
954          pwrf_grid%lalbksoil(ix,iz,iy) = pfalb1(ixy,iz)
955          pwrf_grid%llwalbksoil(ix,iz,iy) = pfalb2(ixy,iz)
956          pwrf_grid%levapksoil(ix,iz,iy) = pevap(ixy,iz)
957          pwrf_grid%lsnowksoil(ix,iz,iy) = psnow(ixy,iz)
958        END DO
959        pwrf_grid%lrads(ix,iy) = pradsol(ixy)
960        pwrf_grid%lsolsw(ix,iy) = psolsw(ixy)
961        pwrf_grid%lsollw(ix,iy) = psollw(ixy)
962        pwrf_grid%lfder(ix,iy) = pfder(ixy)
963        pwrf_grid%lrain(ix,iy) = prain_fall(ixy)
964        pwrf_grid%lsnow(ix,iy) = psnow_fall(ixy)
965        DO iz=1,Nks
966          pwrf_grid%lrugksoil(ix,iz,iy) = prugos(ixy,iz)
967          pwrf_grid%lagesnoksoil(ix,iz,iy) = pagesno(ixy,iz)
968        END DO
969! These variables do not change on time!
970!!
971!    pzmea, pzstd, pzsig, pzgam, pzthe, pzpic, pzval, prugsrel
972        pwrf_grid%lrugsea(ix,iy) = prugos(ixy,is_oce)
973!        pwrf_grid%lrunofflic(ix,iy)= prun_off_lic(ixy)
974        pwrf_grid%lzmax0(ix,iy) = pzmax0(ixy)
975        pwrf_grid%lf0(ix,iy) = pf0(ixy)
976        pwrf_grid%lwake_s(ix,iy) = pwake_s(ixy)
977        pwrf_grid%lwake_cstar(ix,iy) = pwake_cstar(ixy)
978        pwrf_grid%lwake_pe(ix,iy) = pwake_pe(ixy)
979        pwrf_grid%lwake_fip(ix,iy) = pwake_fip(ixy)
980        pwrf_grid%lnat(ix,iy) = pnat(ixy)
981        pwrf_grid%sst(ix,iy) = psst(ixy)
982!        pwrf_grid%lbils(ix,iy) = pbils(ixy)
983        pwrf_grid%albedo(ix,iy) = palbedo(ixy)
984!        pwrf_grid%lrug(ix,iy) = prug(ixy)
985      END DO
986    END DO
987
988    dx=ddx+2*bdy
989    dy=ddy+2*bdy
990
991    CALL vect_mat(pclwcon, dx, dy, ddz, bdy, pwrf_grid%lclwcon)
992    CALL vect_mat(prnebcon, dx, dy, ddz, bdy, pwrf_grid%lrnebcon)
993    CALL vect_mat(pratqs, dx, dy, ddz, bdy, pwrf_grid%lratqs)
994    CALL vect_mat(pema_work1, dx, dy, ddz, bdy, pwrf_grid%lema_work1)
995    CALL vect_mat(pema_work2, dx, dy, ddz, bdy, pwrf_grid%lema_work2)
996    CALL vect_mat(pwake_deltat, dx, dy, ddz, bdy, pwrf_grid%lwake_deltat)
997    CALL vect_mat(pwake_deltaq, dx, dy, ddz, bdy, pwrf_grid%lwake_deltaq)
998    CALL vect_mat(pfm_therm, dx, dy, ddz+1, bdy, pwrf_grid%lfm_therm)
999    CALL vect_mat(pentr_therm, dx, dy, ddz, bdy, pwrf_grid%lentr_therm)
1000    CALL vect_mat(pdetr_therm, dx, dy, ddz, bdy, pwrf_grid%ldetr_therm)
1001
1002    RETURN
1003
1004  END SUBROUTINE get_lmdz
1005
1006  SUBROUTINE get_lowbdy(ddx, ddy, ddz, bdy, dlmdz, Nks, Nz0, pkglo, pwrf_grid,       &
1007    ppctsrf, pftsol, prugos, palbedo, psst, ptsurf_lim, pz0_lim, palb_lim, prug_glo, &
1008    ppct_glo)
1009! Subroutine to load LMDZ limit variables with values from lowbdy WRF
1010
1011! WRF modules
1012!    USE module_model_constants
1013    USE module_domain_type
1014
1015! LMDZ modules
1016    USE indice_sol_mod
1017
1018    IMPLICIT NONE
1019
1020    TYPE(domain) , TARGET                                :: pwrf_grid
1021    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
1022      dlmdz, Nks, Nz0, pkglo
1023    REAL, DIMENSION(dlmdz), INTENT(INOUT)                :: psst, palbedo,           &
1024      ptsurf_lim, pz0_lim, palb_lim
1025! What to do with these variables? ,pnat,pbils
1026    REAL, DIMENSION(dlmdz,Nks), INTENT(INOUT)            :: ppctsrf, pftsol, prugos
1027    REAL, DIMENSION(pkglo), INTENT(OUT)                  :: prug_glo
1028    REAL, DIMENSION(pkglo,Nks), INTENT(OUT)              :: ppct_glo
1029
1030! Local
1031    INTEGER                                              :: ix,iy,iz,ixy,dx,dy,ixx,iyy
1032    CHARACTER(LEN=50)                                    :: LMDZvarmethod,fname,errmsg
1033
1034    fname='get_lowbdy'
1035    errmsg='ERROR -- error -- ERROR -- error'
1036
1037! Filling all 2D variables
1038!!
1039    DO iy=1,ddy
1040      DO ix=1,ddx
1041        ixy=ddx*(iy-1)+ix
1042        ppctsrf(ixy,is_ter) = pwrf_grid%lter(ix,iy)
1043        ppctsrf(ixy,is_lic) = pwrf_grid%llic(ix,iy)
1044        IF ( pwrf_grid%xice(ix,iy) /= pwrf_grid%lsic(ix,iy) ) THEN
1045! No ocean/frozen ocean fraction before
1046          IF (pwrf_grid%loce(ix,iy) == 0.) pwrf_grid%lrugksoil(ix,3,iy) = 0.0001
1047          IF (pwrf_grid%lsic(ix,iy) == 0.) pwrf_grid%lrugksoil(ix,4,iy) = 0.001
1048          pwrf_grid%loce(ix,iy) = 1. - pwrf_grid%lsic(ix,iy)
1049          pwrf_grid%lsic(ix,iy) = pwrf_grid%lsic(ix,iy)
1050          ppctsrf(ixy,is_oce) = pwrf_grid%loce(ix,iy)
1051          ppctsrf(ixy,is_sic) = pwrf_grid%lsic(ix,iy)
1052! Recomputing total roughness
1053!!
1054          pwrf_grid%lrug(ix,iy)=pwrf_grid%lter(ix,iy)*pwrf_grid%lrugksoil(ix,1,iy) +&
1055            pwrf_grid%llic(ix,iy)*pwrf_grid%lrugksoil(ix,2,iy) +                    &
1056            pwrf_grid%loce(ix,iy)*pwrf_grid%lrugksoil(ix,3,iy) +                    &
1057            pwrf_grid%lsic(ix,iy)*pwrf_grid%lrugksoil(ix,4,iy)
1058        END IF
1059! Full SST grid points must have a SST value from the wrfinput (WRF 1/0 landmask)
1060!   but that LMDZ points with fractions of land/mask could not. Using the closest one
1061!   if it is an isolated grid point, using TSK instead
1062        IF (pwrf_grid%sst(ix,iy) < 200.15) THEN
1063          pwrf_grid%ltksoil(ix,3,iy) = -9.
1064          IF ( (ix > 1).AND.(ix < ddx) .AND. (iy > 1).AND.(iy < ddy) ) THEN
1065            DO ixx=-1,1
1066              DO iyy=-1,1
1067!               PRINT *,ixx,iyy,wrf_grid%sst(ix+ixx,iy+iyy)
1068                IF (pwrf_grid%sst(ix+ixx,iy+iyy) > 200.15) THEN
1069                  pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%sst(ix+ixx,iy+iyy)
1070                  EXIT
1071                END IF
1072              END DO
1073              IF ( pwrf_grid%ltksoil(ix,3,iy) > 200.15) EXIT
1074            END DO
1075            IF ( pwrf_grid%ltksoil(ix,3,iy) == -9.) pwrf_grid%ltksoil(ix,3,iy) =     &
1076              pwrf_grid%tsk(ix,iy)
1077            ELSE
1078              pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%tsk(ix,iy)
1079            END IF
1080        ELSE
1081          pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%sst(ix,iy)
1082        END IF
1083        IF ( (pwrf_grid%ltksoil(ix,3,iy) < 200.15) .OR.                              &
1084          (pwrf_grid%ltksoil(ix,3,iy) > 370.) ) THEN
1085          PRINT *,TRIM(errmsg)
1086!          WRITE(wrf_err_message,*) '  ' // TRIM(fname) //                            &
1087!            ' SST inconsistency at point: ', ix, ', ', iy,' with loce: ',            &
1088!            wrf_grid%loce(ix,iy),' has a tsoil(oce)= ',                              &
1089!            wrf_grid%ltksoil(ix,3,iy),' sst: ',wrf_grid%sst(ix,iy),' tsk: ',         &
1090!            wrf_grid%tsk(ix,iy)
1091!          CALL wrf_error_fatal(TRIM(wrf_err_message))
1092          PRINT *, '  ' // TRIM(fname) //                                            &
1093          ' SST inconsistency at point: ', ix, ', ', iy,' with loce: ',              &
1094            pwrf_grid%loce(ix,iy),' has a tsoil(oce)= ',                             &
1095            pwrf_grid%ltksoil(ix,3,iy),' sst: ', pwrf_grid%sst(ix,iy),' tsk: ',      &
1096            pwrf_grid%tsk(ix,iy)
1097        END IF
1098        psst(ixy) = pwrf_grid%ltksoil(ix,3,iy)
1099        palbedo(ixy) = pwrf_grid%albbck(ix,iy)
1100        DO iz=1,Nks
1101          prugos(ixy,iz) = pwrf_grid%lrugksoil(ix,iz,iy)
1102        END DO
1103
1104! In ocean_forced_mod is used as tsurf_lim = sst(knindex(1:knon)). Let's complicate it
1105!   a little bit...
1106        ptsurf_lim(ixy)= psst(ixy)
1107! In surf_land_bucket_mod is used as z0_new=rugos(knindex(1:knon),1) [z0_new == z0_limit]
1108        pz0_lim(ixy) = prugos(ixy,1)
1109        palb_lim(ixy) = palbedo(ixy)
1110
1111        prug_glo(ixy) = prugos(ixy,is_ter)
1112        ppct_glo(ixy,:) = ppctsrf(ixy,:)
1113
1114      END DO
1115    END DO
1116
1117    RETURN
1118
1119  END SUBROUTINE get_lowbdy
1120
1121  SUBROUTINE load_lmdz(pwrf_grid,ddx,ddy,ddz,bdy,dlmdz,Nks,NzO,prlon,prlat,pzmasq,   &
1122    ppctsrf,pftsol,pftsoil,pqsurf,pqsol,pfalb1,pfalb2,pevap,psnow,pradsol,psolsw,    &
1123    psollw,pfder,prain_fall,psnow_fall,pagesno,pzmea,pzstd,pzgam,pzthe,pzpic,pzval,  &
1124    prugoro,prugos,prun_off_lic,pzmax0,pf0,pwake_s,pwake_cstar,pwake_pe,pwake_fip,   &
1125    psst,palbedo,pt_ancien,pq_ancien,pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,    &
1126    pema_work1,pema_work2,pwake_deltat,pwake_deltaq,pfm_therm,pentr_therm,           &
1127    pdetr_therm,pnat)
1128! Subroutine to load LMDZ variables with values from WRF
1129
1130! WRF modules
1131!    USE module_model_constants
1132    USE module_domain_type
1133
1134! LMDZ modules
1135    USE indice_sol_mod
1136
1137    IMPLICIT NONE
1138
1139    TYPE(domain) , TARGET                                :: pwrf_grid
1140    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
1141      dlmdz, Nks, NzO
1142    REAL, DIMENSION(dlmdz), INTENT(OUT)                  :: prlon, prlat, pzmasq,    &
1143      pqsol, pradsol, psolsw, psollw, pfder, prain_fall, psnow_fall, pzmea, pzstd,   &
1144      pzgam, pzthe, pzpic, pzval, prugoro, prun_off_lic, pzmax0, pf0, pwake_s,       &
1145      pwake_cstar, pwake_pe, pwake_fip, psst, palbedo,pnat
1146    REAL, DIMENSION(dlmdz,Nks), INTENT(OUT)              :: ppctsrf,pftsol,pqsurf,   &
1147      pfalb1,pfalb2,pevap,psnow,prugos,pagesno
1148    REAL, DIMENSION(dlmdz,ddz), INTENT(OUT)              :: pt_ancien,pq_ancien,     &
1149      pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,pema_work1,pema_work2,pwake_deltat,&
1150      pwake_deltaq,pentr_therm,pdetr_therm
1151    REAL, DIMENSION(dlmdz,NzO,Nks), INTENT(OUT)          :: pftsoil
1152    REAL, DIMENSION(dlmdz,ddz+1), INTENT(OUT)            :: pfm_therm
1153
1154
1155! Local
1156    INTEGER                                              :: ix,iy,iz,ixy,dx,dy
1157    CHARACTER(LEN=50)                                    :: LMDZvarmethod,fname
1158
1159    fname='load_lmdz'
1160
1161! Filling all 2D variables
1162!!
1163    PRINT *,'  ' // TRIM(fname) // '; LUBOUND(grid%lmsk): ',LBOUND(pwrf_grid%lmsk), UBOUND(pwrf_grid%lmsk)
1164    PRINT *,'  ' // TRIM(fname) // '; ddx ddy: ',ddy,ddx
1165
1166    DO iy=1,ddy
1167      DO ix=1,ddx
1168        ixy=ddx*(iy-1)+ix
1169        prlon(ixy) = pwrf_grid%xlong(ix,iy)
1170        prlat(ixy) = pwrf_grid%xlat(ix,iy)
1171        pzmasq(ixy) = pwrf_grid%lmsk(ix,iy)
1172        ppctsrf(ixy,is_ter) = pwrf_grid%lter(ix,iy)
1173        ppctsrf(ixy,is_lic) = pwrf_grid%llic(ix,iy)
1174        ppctsrf(ixy,is_oce) = pwrf_grid%loce(ix,iy)
1175        ppctsrf(ixy,is_sic) = pwrf_grid%lsic(ix,iy)
1176        DO iz=1,Nks
1177          pftsol(ixy,iz) = pwrf_grid%ltksoil(ix,iz,iy)
1178        END DO
1179        DO iz=1,NzO
1180          pftsoil(ixy,iz,1) = pwrf_grid%lotter(ix,iz,iy)
1181          pftsoil(ixy,iz,2) = pwrf_grid%lotlic(ix,iz,iy)
1182          pftsoil(ixy,iz,3) = pwrf_grid%lotoce(ix,iz,iy)
1183          pftsoil(ixy,iz,4) = pwrf_grid%lotsic(ix,iz,iy)
1184        END DO
1185        DO iz=1,Nks
1186          pqsurf(ixy,iz) = pwrf_grid%lqksoil(ix,iz,iy)
1187        END DO
1188        pqsol(ixy) = pwrf_grid%lwsol(ix,iy)
1189        DO iz=1,Nks
1190          pfalb1(ixy,iz) = pwrf_grid%lalbksoil(ix,iz,iy)
1191          pfalb2(ixy,iz) = pwrf_grid%llwalbksoil(ix,iz,iy)
1192          pevap(ixy,iz) = pwrf_grid%levapksoil(ix,iz,iy)
1193          psnow(ixy,iz) = pwrf_grid%lsnowksoil(ix,iz,iy)
1194        END DO
1195! No way to get back the independent radiative values (unless using LMDZ 'hist*'
1196!   files). Use it in the LMDZ's added way also from grid (no need to initialize it,
1197!   no radiation on t=0)
1198        pradsol(ixy) = pwrf_grid%lrads(ix,iy)
1199        psolsw(ixy) = pwrf_grid%lsolsw(ix,iy)
1200        psollw(ixy) = pwrf_grid%lsollw(ix,iy)
1201        pfder(ixy) = pwrf_grid%lfder(ix,iy)
1202! LMDZ does not provide a way (in the restart) to distingusih between convective,
1203!    non-convective. Thus Use it as a whole
1204        prain_fall(ixy) = pwrf_grid%lrain(ix,iy)
1205        psnow_fall(ixy) = pwrf_grid%lsnow(ix,iy)
1206        DO iz=1,Nks
1207          prugos(ixy,iz) = pwrf_grid%lrugksoil(ix,iz,iy)
1208          pagesno(ixy,iz) = pwrf_grid%lagesnoksoil(ix,iz,iy)
1209        END DO
1210        pzmea(ixy) = pwrf_grid%lzmea(ix,iy)
1211        pzstd(ixy) = pwrf_grid%lzstd(ix,iy)
1212        pzgam(ixy) = pwrf_grid%lzgam(ix,iy)
1213        pzthe(ixy) = pwrf_grid%lzthe(ix,iy)
1214        pzpic(ixy) = pwrf_grid%lzpic(ix,iy)
1215        pzval(ixy) = pwrf_grid%lzval(ix,iy)
1216        prugoro(ixy) = pwrf_grid%lzrugsrel(ix,iy)
1217!        prugos(ixy,is_oce) = pwrf_grid%lrugsea(ix,iy)
1218        prun_off_lic(ixy) = pwrf_grid%lrunofflic(ix,iy)
1219        pzmax0(ixy) = pwrf_grid%lzmax0(ix,iy)
1220        pf0(ixy) = pwrf_grid%lf0(ix,iy)
1221        pwake_s(ixy) = pwrf_grid%lwake_s(ix,iy)
1222        pwake_cstar(ixy) = pwrf_grid%lwake_cstar(ix,iy)
1223        pwake_pe(ixy) = pwrf_grid%lwake_pe(ix,iy)
1224        pwake_fip(ixy) = pwrf_grid%lwake_fip(ix,iy)
1225!! limit.nc
1226! nat == pctsrf, not needed ?
1227        pnat(ixy) = pwrf_grid%lnat(ix,iy)
1228        psst(ixy) = pwrf_grid%sst(ix,iy)
1229! not used, not defined....
1230!        pbils(ixy) = pwrf_grid%lbils(ix,iy)
1231        palbedo(ixy) = pwrf_grid%albedo(ix,iy)
1232! already done in previous steps
1233!        prugos(ixy) = pwrf_grid%lrug(ix,iy)
1234
1235      END DO
1236    END DO
1237
1238    dx=ddx+2*bdy
1239    dy=ddy+2*bdy
1240
1241    CALL mat_vect(pwrf_grid%t_2, dx, dy, ddz, bdy, pt_ancien)
1242! New variable added in the Registry qv_2
1243    CALL mat_vect(pwrf_grid%qv_2, dx, dy, ddz, bdy, pq_ancien)
1244    CALL mat_vect(pwrf_grid%u_2, dx, dy, ddz, bdy, pu_ancien)
1245    CALL mat_vect(pwrf_grid%v_2, dx, dy, ddz, bdy, pv_ancien)
1246    CALL mat_vect(pwrf_grid%lclwcon, dx, dy, ddz, bdy, pclwcon)
1247    CALL mat_vect(pwrf_grid%lrnebcon, dx, dy, ddz, bdy, prnebcon)
1248    CALL mat_vect(pwrf_grid%lratqs, dx, dy, ddz, bdy, pratqs)
1249    CALL mat_vect(pwrf_grid%lema_work1, dx, dy, ddz, bdy, pema_work1)
1250    CALL mat_vect(pwrf_grid%lema_work2, dx, dy, ddz, bdy, pema_work2)
1251    CALL mat_vect(pwrf_grid%lwake_deltat, dx, dy, ddz, bdy, pwake_deltat)
1252    CALL mat_vect(pwrf_grid%lwake_deltaq, dx, dy, ddz, bdy, pwake_deltaq)
1253    CALL mat_vect(pwrf_grid%lfm_therm, dx, dy, ddz+1, bdy, pfm_therm)
1254    CALL mat_vect(pwrf_grid%lentr_therm, dx, dy, ddz, bdy, pentr_therm)
1255    CALL mat_vect(pwrf_grid%ldetr_therm, dx, dy, ddz, bdy, pdetr_therm)
1256
1257    RETURN
1258
1259  END SUBROUTINE load_lmdz
1260
1261  SUBROUTINE vect_mat(vect, dx, dy, dz, wbdy, mat)
1262! Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
1263
1264    IMPLICIT NONE
1265
1266    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1267    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(IN) :: vect
1268    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1), INTENT(OUT) :: mat
1269
1270! Local
1271    INTEGER                                              :: i,j,k
1272    INTEGER                                              :: ddx, ddy
1273    INTEGER                                              :: lp, il, jl
1274
1275! Variables
1276! d[x/y/z]: dimensions
1277! mat: matrix to transform
1278! vect: vector to produce
1279    mat=0.
1280
1281    ddx=dx-2*wbdy
1282    ddy=dy-2*wbdy
1283!    lp=321
1284!    jl=INT(lp/ddx)+1
1285!    il=lp-ddx*(jl-1)
1286
1287    DO k=1,dz
1288      DO j=1,ddy
1289        DO i=1,ddx
1290          mat(i,k,j)=vect(ddx*(j-1)+i,k)
1291        END DO
1292      END DO
1293    END DO
1294
1295  END SUBROUTINE vect_mat
1296
1297
1298  SUBROUTINE vect_mat_zstagg(vect, dx, dy, dz, wbdy, mat)
1299! Subroutine to transform a LMDZ physic 1D vector to a 3D matrix
1300
1301    IMPLICIT NONE
1302
1303    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1304    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(IN) :: vect
1305    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz,1-wbdy:dy-wbdy+1), INTENT(OUT) :: mat
1306
1307! Local
1308    INTEGER                                              :: i,j,k
1309    INTEGER                                              :: ddx, ddy
1310
1311! Variables
1312! d[x/y/z]: dimensions
1313! mat: matrix to transform
1314! vect: vector to produce
1315
1316    ddx=dx-2*wbdy
1317    ddy=dy-2*wbdy
1318
1319    mat=0.
1320
1321    DO k=1,dz
1322      DO j=1,ddy
1323        DO i=1,ddx
1324          mat(i,k,j)=vect(ddx*(j-1)+i,k)
1325        END DO
1326      END DO
1327    END DO
1328
1329  END SUBROUTINE vect_mat_zstagg
1330
1331
1332  SUBROUTINE mat_vect(mat, dx, dy, dz, wbdy, vect)
1333! Subroutine to transform a 3D matrix to a LMDZ physic 1D vector
1334
1335    IMPLICIT NONE
1336
1337    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1338    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1), INTENT(IN) :: mat
1339    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(OUT) ::                      &
1340      vect
1341
1342! Local
1343    INTEGER                                              :: i,j,k
1344    INTEGER                                              :: ddx,ddy
1345
1346! Variables
1347! d[x/y/z]: dimensions
1348! mat: matrix to transform
1349! vect: vector to produce
1350
1351! Removing HALO
1352    ddx=dx-2*wbdy
1353    ddy=dy-2*wbdy
1354
1355    vect=0.
1356
1357    DO k=1,dz
1358      DO j=1,ddy
1359        DO i=1,ddx
1360          vect(ddx*(j-1)+i,k)= mat(i,k,j)
1361        END DO
1362      END DO
1363    END DO
1364
1365  END SUBROUTINE mat_vect
1366
1367  SUBROUTINE mat_vect_zstagg(mat, dx, dy, dz, wbdy, vect)
1368! Subroutine to transform a 3D matrix z-staggered to a LMDZ physic 1D vector
1369
1370    IMPLICIT NONE
1371
1372    INTEGER, INTENT(IN)                                  :: dx, dy, dz, wbdy
1373    REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz,1-wbdy:dy-wbdy+1), INTENT(IN) :: mat
1374    REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(OUT) ::                      &
1375      vect
1376
1377! Local
1378    INTEGER                                              :: i,j,k
1379    INTEGER                                              :: ddx,ddy
1380
1381! Variables
1382! d[x/y/z]: dimensions
1383! mat: matrix to transform
1384! vect: vector to produce
1385
1386! Removing HALO
1387    ddx=dx-2*wbdy
1388    ddy=dy-2*wbdy
1389
1390    DO k=1,dz
1391      DO j=1,ddy
1392        DO i=1,ddx
1393          vect(ddx*(j-1)+i,k)= mat(i,k,j)
1394        END DO
1395      END DO
1396    END DO
1397
1398  END SUBROUTINE mat_vect_zstagg
1399
1400  SUBROUTINE WRFlanduse_LMDZKsoil(lndcn,is,ie,js,je,dx,dy,ncat,lndu,mask,kt,kl,ko,ks)
1401! Subroutine to retrieve LMDZ Ksoil classification using WRF landuse
1402
1403    IMPLICIT NONE
1404
1405    CHARACTER(LEN=50), INTENT(IN)                        :: lndcn
1406    INTEGER, INTENT(IN)                                  :: is,ie,js,je,dx,dy,ncat
1407    REAL, DIMENSION(is:ie,ncat,js:je), INTENT(IN)        :: lndu
1408    REAL, DIMENSION(is:ie,js:je), INTENT(IN)             :: mask
1409    REAL, DIMENSION(is:ie,js:je), INTENT(OUT)            :: kt,kl,ko,ks
1410
1411! Local
1412    INTEGER                                              :: i,j,k,ij,Nnones,Ncattot
1413    CHARACTER(LEN=50)                                    :: fname, errmsg, wrnmsg
1414    REAL, DIMENSION(is:ie,js:je)                         :: totlndu
1415    INTEGER, ALLOCATABLE, DIMENSION(:)                   :: ter_wk, lic_wk, oce_wk,  &
1416      sic_wk
1417    LOGICAL, ALLOCATABLE, DIMENSION(:,:)                 :: wk
1418    INTEGER                                              :: Nter_wk, Nlic_wk,        &
1419        Noce_wk, Nsic_wk, klks_equal
1420    REAL*8, DIMENSION(is:ie,js:je)                       :: Ksoilval
1421
1422
1423!!!!!!!! Variables
1424! lndcn: land-use data-set (same as in VEGPARM.TBL)
1425! dx,dy: horizontal dimension of the domain
1426! ncat: number of categories of the data-set
1427! lndu: WRF landuse matrix
1428! mask: land-sea mask from WRF
1429! k[t/l/o/s]: LMDZ ter, lic, oce, sic matrices
1430! Ncattot: theoretical number of laduses for a given data-base
1431! [ter/lic/oce/sic]_wk: vectors with the indices of the landuse for LMDZ Ksoil
1432! N[ter/lic/oce/sic]_wk: number of indices of the landuse for LMDZ Ksoil
1433! wk: boolean matrix for each Ksoil indicating 'true' if the WRF landuse is for Ksoil
1434
1435    fname="'WRFlanduse_LMDZKsoil'"
1436    errmsg='ERROR -- error -- ERROR -- error'
1437    wrnmsg='WARNING -- warning -- WARNING -- warning'
1438   
1439! Checking if at all the grid points SUM(landuse) == 1.
1440!!
1441    PRINT *,'  ' // TRIM(fname) // " using '" // TRIM(lndcn) // "' WRF land types"
1442
1443    totlndu = SUM(lndu, 2)
1444
1445    Nnones=0
1446    DO i=1,dx
1447      DO j=1,dy
1448        IF (totlndu(i,j) /= 1. ) Nnones = Nnones+1
1449      END DO
1450    END DO
1451
1452    IF (Nnones /= 0) THEN
1453      PRINT *,TRIM(errmsg)
1454      PRINT *,'    ' // TRIM(fname) // ': ', Nnones,' grid points do not have a ' // &
1455        'TOTlanduse == 1. from WRF !!!!!'
1456      DO i=1,dx
1457        DO j=1,dy
1458          IF (totlndu(i,j) /= 1.) PRINT *,'      ',i,', ',j,' total land use: ',     &
1459            totlndu(i,j),' indiv. values :', lndu(i,:,j)
1460        END DO
1461      END DO
1462    END IF
1463
1464    SELECT CASE (lndcn)
1465      CASE ('OLD')
1466        Ncattot = 13
1467        Nter_wk = 11
1468        Nlic_wk = 1
1469        Noce_wk = 1
1470        Nsic_wk = 1
1471
1472        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1473        ALLOCATE(ter_wk(Nter_wk))
1474        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1475        ALLOCATE(lic_wk(Nlic_wk))
1476        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1477        ALLOCATE(oce_wk(Noce_wk))
1478        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1479        ALLOCATE(sic_wk(Nsic_wk))
1480
1481! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1482!!
1483!  1,'Urban land'                    2,'Agriculture'
1484!  3,'Range-grassland'               4,'Deciduous forest'
1485!  5,'Coniferous forest'             6,'Mixed forest and wet land'
1486!  7,'Water'                         8,'Marsh or wet land'
1487!  9,'Desert'                       10,'Tundra'
1488! 11,'Permanent ice'                12,'Tropical or subtropical forest'
1489! 13,'Savannah'
1490
1491        ter_wk = (/ 1,2,3,4,5,6,8,9,10,12,13 /)
1492        lic_wk = (/ 11 /)
1493        oce_wk = (/ 7 /)
1494        sic_wk = (/ 11 /)
1495
1496      CASE ('USGS')
1497        Ncattot = 33
1498        Nter_wk = 22
1499        Nlic_wk = 1
1500        Noce_wk = 1
1501        Nsic_wk = 1
1502
1503        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1504        ALLOCATE(ter_wk(Nter_wk))
1505        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1506        ALLOCATE(lic_wk(Nlic_wk))
1507        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1508        ALLOCATE(oce_wk(Noce_wk))
1509        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1510        ALLOCATE(sic_wk(Nsic_wk))
1511
1512! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1513!!
1514! 1,'Urban and Built-Up Land'        2,'Dryland Cropland and Pasture'
1515! 3,'Irrigated Cropland and Pasture' 4,'Mixed Dryland/Irrigated Cropland and Pasture'
1516! 5,'Cropland/Grassland Mosaic'      6,'Cropland/Woodland Mosaic'
1517! 7,'Grassland'                      8,'Shrubland'
1518! 9,'Mixed Shrubland/Grassland'     10,'Savanna'
1519!11,'Deciduous Broadleaf Forest'    12,'Deciduous Needleleaf Forest'
1520!13,'Evergreen Broadleaf Forest'    14,'Evergreen Needleleaf Forest'
1521!15,'Mixed Forest'                  16,'Water Bodies'
1522!17,'Herbaceous Wetland'            18,'Wooded Wetland'
1523!19,'Barren or Sparsely Vegetated'  20,'Herbaceous Tundra'
1524!21,'Wooded Tundra'                 22,'Mixed Tundra'
1525!23,'Bare Ground Tundra'            24,'Snow or Ice'
1526!25,'Playa'                         26,'Lava'
1527!27,'White Sand'                    28,'Unassigned'
1528!29,'Unassigned'                    30,'Unassigned'
1529!31,'Low Intensity Residential '    32,'High Intensity Residential'
1530!33,'Industrial or Commercial'     
1531
1532! Correct values
1533!        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,  &
1534!          27,31,32,33 /)
1535! How ever, USGS by default has 24 so...
1536        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,18,19,20,21,22,23 /)
1537        lic_wk = (/ 24 /)
1538        oce_wk = (/ 16 /)
1539        sic_wk = (/ 24 /)
1540
1541      CASE ('MODIFIED_IGBP_MODIS_NOAH')
1542        Ncattot = 33
1543        Nter_wk = 28
1544        Nlic_wk = 1
1545        Noce_wk = 1
1546        Nsic_wk = 1
1547
1548        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1549        ALLOCATE(ter_wk(Nter_wk))
1550        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1551        ALLOCATE(lic_wk(Nlic_wk))
1552        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1553        ALLOCATE(oce_wk(Noce_wk))
1554        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1555        ALLOCATE(sic_wk(Nsic_wk))
1556
1557! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1558!!
1559!  1,'Evergreen Needleleaf Forest'             2,'Evergreen Broadleaf Forest'
1560!  3,'Deciduous Needleleaf Forest'             4,'Deciduous Broadleaf Forest'
1561!  5,'Mixed Forests'                           6,'Closed Shrublands'
1562!  7,'Open Shrublands'                         8,'Woody Savannas'
1563!  9,'Savannas'                               10,'Grasslands'
1564! 11,'Permanent wetlands'                     12,'Croplands'
1565! 13,'Urban and Built-Up'                     14,'cropland/natural vegetation mosaic'
1566! 15,'Snow and Ice'                           16,'Barren or Sparsely Vegetated'
1567! 17,'Water'                                  18,'Wooded Tundra'
1568! 19,'Mixed Tundra'                           20,'Barren Tundra'
1569! 21,'Unassigned'                             22,'Unassigned'
1570! 23,'Unassigned'                             24,'Unassigned'
1571! 25,'Unassigned'                             26,'Unassigned'
1572! 27,'Unassigned'                             28,'Unassigned'
1573! 29,'Unassigned'                             30,'Unassigned'
1574! 31,'Low Intensity Residential '             32,'High Intensity Residential'
1575! 33,'Industrial or Commercial'
1576
1577        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,19,20,31,32,33 /)
1578        lic_wk = (/ 15 /)
1579        oce_wk = (/ 17 /)
1580        sic_wk = (/ 15 /)
1581
1582      CASE ('SiB')
1583        Ncattot = 16
1584        Nter_wk = 14
1585        Nlic_wk = 1
1586        Noce_wk = 1
1587        Nsic_wk = 1
1588
1589        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1590        ALLOCATE(ter_wk(Nter_wk))
1591        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1592        ALLOCATE(lic_wk(Nlic_wk))
1593        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1594        ALLOCATE(oce_wk(Noce_wk))
1595        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1596        ALLOCATE(sic_wk(Nsic_wk))
1597
1598! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1599!!
1600!  1,'Evergreen Broadleaf Trees'        2,'Broadleaf Deciduous Trees'
1601!  3,'Deciduous and Evergreen Trees'    4,'Evergreen Needleleaf Trees'
1602!  5,'Deciduous Needleleaf Trees'       6,'Ground Cover with Trees and Shrubs'
1603!  7,'Groundcover Only'              8,'Broadleaf Shrubs with Perennial Ground Cover'
1604!  9,'Broadleaf Shrubs with Bare Soil' 10,'Groundcover with Dwarf Trees and Shrubs'
1605! 11,'Bare Soil'                       12,'Agriculture or C3 Grassland'
1606! 13,'Persistent Wetland'              14,'Dry Coastal Complexes'
1607! 15,'Water'                           16,'Ice Cap and Glacier'
1608
1609        ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14 /)
1610        lic_wk = (/ 16 /)
1611        oce_wk = (/ 15 /)
1612        sic_wk = (/ 16 /)
1613
1614      CASE ('LW12')
1615        Ncattot = 3
1616        Nter_wk = 1
1617        Nlic_wk = 1
1618        Noce_wk = 1
1619        Nsic_wk = 1
1620
1621        IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk)
1622        ALLOCATE(ter_wk(Nter_wk))
1623        IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk)
1624        ALLOCATE(lic_wk(Nlic_wk))
1625        IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk)
1626        ALLOCATE(oce_wk(Noce_wk))
1627        IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk)
1628        ALLOCATE(sic_wk(Nsic_wk))
1629
1630! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3)
1631!!
1632!  1,'Land'                              2,'Water'
1633!  3,'Snow or Ice'
1634
1635        ter_wk = (/ 1 /)
1636        lic_wk = (/ 3 /)
1637        oce_wk = (/ 2 /)
1638        sic_wk = (/ 3 /)
1639
1640      CASE DEFAULT
1641        PRINT *,TRIM(errmsg)
1642        PRINT *,'  ' // TRIM(fname) // ": landuse data-base '" // TRIM(lndcn) //     &
1643          "' not ready !!!"
1644        STOP
1645
1646    END SELECT
1647
1648! Checking the correct number of categories
1649    IF (Ncattot /= ncat) THEN
1650      PRINT *,TRIM(wrnmsg)
1651      PRINT *,'  ' // TRIM(fname) // ": '" // TRIM(lndcn) // "' landuse data-base "  &
1652        // ' has ',Ncattot,' land-uses types ',ncat,' passed !!!!'
1653!      STOP
1654    END IF
1655
1656! Mimic WRF landuse with 4 LMDZ Ksoils
1657!!
1658    IF (ALLOCATED(wk)) DEALLOCATE(wk)
1659!    ALLOCATE(wk(4,ncattot))
1660    ALLOCATE(wk(4,ncat))
1661
1662    wk = .FALSE.
1663
1664    DO i=1,Nter_wk
1665      wk(1,ter_wk(i)) = .TRUE.
1666    END DO
1667    DO i=1,Nlic_wk
1668      wk(2,lic_wk(i)) = .TRUE.
1669    END DO
1670    DO i=1,Noce_wk
1671      wk(3,oce_wk(i)) = .TRUE.
1672    END DO
1673    DO i=1,Nsic_wk
1674      wk(4,sic_wk(i)) = .TRUE.
1675    END DO
1676
1677    klks_equal = 0
1678    DO i=1,dx
1679      DO j=1,dy
1680! Using double precission numbers to avoid truncation errors...       
1681        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(1,:))
1682        kt(i,j) = Ksoilval(1,1)
1683        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(2,:))
1684        kl(i,j) = Ksoilval(1,1)
1685        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(3,:))
1686        ko(i,j) = Ksoilval(1,1)
1687        Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(4,:))
1688        ks(i,j) = Ksoilval(1,1)
1689
1690        IF (kl(i,j) == ks(i,j)) klks_equal = klks_equal+1
1691      END DO
1692    END DO
1693
1694! Assigning lic/sic to all that ice/snow according to the percentage of ter/oce
1695!   Is there any other better way? Satellites, do not know what is below the ice ?!
1696
1697    IF (klks_equal == dx*dy) THEN
1698!   Assuming that kl=ks=ice
1699      DO i=1,dx
1700        DO j=1,dy
1701          IF ((kt(i,j)+ko(i,j)) /= 0.) THEN
1702            Ksoilval(1,1) = kl(i,j)*kt(i,j)/(kt(i,j)+ko(i,j))
1703            kl(i,j) = Ksoilval(1,1)
1704            Ksoilval(1,1) = ks(i,j)*ko(i,j)/(kt(i,j)+ko(i,j))
1705            ks(i,j) = Ksoilval(1,1)
1706          ELSE
1707            Ksoilval(1,1) = kl(i,j)/(kl(i,j)+ks(i,j))
1708            Ksoilval(1,2) = ks(i,j)/(kl(i,j)+ks(i,j))
1709            kl(i,j) = Ksoilval(1,1)
1710            ks(i,j) = Ksoilval(1,2)
1711          END IF
1712! Introducing WRF land-sea mask
1713          IF (mask(i,j) == 1.) THEN
1714            kl(i,j)=kl(i,j)+ks(i,j)
1715            ks(i,j)=0.
1716          ELSE
1717            ks(i,j)=kl(i,j)+ks(i,j)
1718            kl(i,j)=0.
1719          END IF
1720        END DO
1721      END DO
1722    ELSE
1723      PRINT *,'  ' // TRIM(fname) // ': land-use data provided different lic & sic !'
1724    END IF
1725
1726! Checking if at all the grid points kt+kl+ko+ks == 1.
1727!!
1728    totlndu = kt+kl+ko+ks
1729
1730    Nnones=0
1731    DO i=1,dx
1732      DO j=1,dy
1733        IF (totlndu(i,j) /= 1. ) Nnones = Nnones+1
1734        IF (kt(i,j) + ko(i,j) == 0.)THEN
1735          PRINT *,TRIM(wrnmsg)
1736          PRINT *,'  ' // TRIM(fname) // ': point without land or water category !!'
1737          PRINT *,'i j WRF_ter ter ______________'
1738          DO k=1,Nter_wk
1739            PRINT *,i,j,lndu(i,ter_wk(k),j), kt(i,j)
1740          END DO
1741          PRINT *,'i j WRF_oce oce ______________'
1742          DO k=1,Noce_wk
1743            PRINT *,i,j,lndu(i,oce_wk(k),j), ko(i,j)
1744          END DO
1745          PRINT *,'  wrf%landusef: ', lndu(i,:,j)
1746!          STOP
1747        END IF
1748      END DO
1749    END DO
1750
1751    ij=1
1752    IF (Nnones /= 0) THEN
1753      PRINT *,TRIM(errmsg)
1754      PRINT *,'  ' // TRIM(fname) // ': ', Nnones,' grid points do not have a ' //   &
1755        'TOTlanduse == 1.!!!!!'
1756      PRINT *,'pt i j ter+lic+oce+sic : ter lic oce sic ______________'
1757      DO i=1,dx
1758        DO j=1,dy
1759          IF (totlndu(i,j) /= 1.) THEN
1760            PRINT *,ij,'; ',i,j,totlndu(i,j),' : ',kt(i,j),kl(i,j),ko(i,j),ks(i,j)
1761            ij=ij+1
1762          END IF
1763        END DO
1764      END DO
1765    END IF
1766
1767    DEALLOCATE(ter_wk)
1768    DEALLOCATE(lic_wk)
1769    DEALLOCATE(oce_wk)
1770    DEALLOCATE(sic_wk)
1771    DEALLOCATE(wk)
1772
1773    RETURN
1774
1775  END SUBROUTINE WRFlanduse_LMDZKsoil
1776
1777  SUBROUTINE interpolate_layers(NvalsA,posA,valsA,NvalsB,posB,valsB)
1778! Subroutine to interpolate values between different layers (constant values along the depth of a layer)
1779
1780    IMPLICIT NONE
1781
1782    INTEGER, INTENT(IN)                                  :: NvalsA,NvalsB
1783    REAL, DIMENSION(NvalsA), INTENT(IN)                  :: posA,valsA
1784    REAL, DIMENSION(NvalsB), INTENT(IN)                  :: posB
1785    REAL, DIMENSION(NvalsB), INTENT(OUT)                 :: valsB
1786
1787! Local
1788    INTEGER                                              :: i,k
1789    REAL                                                 :: xA,yA,xB,yB
1790
1791!!!!!!! Variables
1792! A: original values used to interpolate
1793! B: values to obtain
1794! Nvals[A/B]: number of values
1795! pos[A/B]: depths of each layer
1796! vals[A/B]: values of each layer
1797
1798!!!!!!! Functions/Subroutines
1799! interpolate_lin: subroutine to linearly interpolate
1800
1801    i=1
1802    DO k=1,NvalsB
1803      IF (posB(k) > posA(i)) i = i+1
1804      IF (i-1 < 1) THEN
1805        xA = posA(i)
1806        yA = valsA(i)
1807        xB = posA(i+1)
1808        yB = valsA(i+1)
1809      ELSEIF ( i > NvalsA) THEN
1810        xA = posA(i-2)
1811        yA = valsA(i-2)
1812        xB = posA(i-1)
1813        yB = valsA(i-1)
1814      ELSE
1815        xA = posA(i-1)
1816        yA = valsA(i-1)
1817        xB = posA(i)
1818        yB = valsA(i)
1819      END IF
1820      IF (posB(k) <= xB .AND. i-1 > 1) THEN
1821        CALL interpolate_lin(xA,yA,xB,yB,posB(k),valsB(k))
1822      ELSEIF (posB(k) <= xA .AND. i-1 < 1) THEN
1823        valsB(k) = yA + (yB - yA)*(posB(k) - xA)/(xB - xA)
1824      ELSE
1825        valsB(k) = yB + (yB - yA)*(posB(k) - xB)/(xB - xA)
1826      END IF
1827!      PRINT *,'i ',i,'  Interpolating for: ',posB(k),' with xA: ', xA,' yA: ',yA,' & ',   &
1828!       ' xB: ', xB,' yB: ',yB,' --> ',valsB(k)
1829    END DO
1830
1831  END SUBROUTINE interpolate_layers
1832
1833  SUBROUTINE interpolate_lin(x1,y1,x2,y2,x0,y)
1834! Program to interpolate values y=a+b*x with: (from clWRF modifications)
1835!       a=y1
1836!       b=(y2-y1)/(x2-x1)
1837!       x=(x2-x0)
1838
1839    IMPLICIT NONE
1840
1841    REAL, INTENT (IN)                                    :: x1,x2,x0
1842!    REAL(r8), INTENT (IN)                                :: y1,y2
1843!    REAL(r8), INTENT (OUT)                               :: y
1844!    REAL(r8)                                             :: a,b,x
1845    REAL, INTENT (IN)                                   :: y1,y2
1846    REAL, INTENT (OUT)                                  :: y
1847    REAL                                                :: a,b,x
1848
1849    a=y1
1850    b=(y2-y1)/(x2-x1)
1851
1852    IF (x0 .GE. x1) THEN
1853      x=x0-x1
1854    ELSE
1855      x=x1-x0
1856      b=-b
1857    ENDIF
1858
1859    y=a+b*x
1860   
1861!    PRINT *,'y=a+b*x'
1862!    IF (b .LT. 0.) THEN
1863!      PRINT *,'y=',y1,' ',b,'*x'
1864!    ELSE
1865!      PRINT *,'y=',y1,'+',b,'*x'
1866!    ENDIF
1867!    PRINT *,'y(',x,')=',y
1868
1869  END SUBROUTINE interpolate_lin
1870
1871  SUBROUTINE table_characteristics(datasetn,tablen,datasetline,Nvar,Nr,Nc)
1872! Subroutine to read values from a WRF .TBL
1873
1874    IMPLICIT NONE
1875
1876    CHARACTER(LEN=50), INTENT(IN)                        :: datasetn
1877    CHARACTER(LEN=100), INTENT(IN)                       :: tablen
1878    INTEGER, INTENT(OUT)                                 :: datasetline, Nvar, Nr, Nc
1879
1880! Local
1881    INTEGER                                              :: il, iline, ic
1882    CHARACTER(LEN=50)                                    :: fname, errmsg
1883    CHARACTER(LEN=200)                                   :: lineval
1884    LOGICAL                                              :: existsfile, is_used
1885    INTEGER                                              :: funit, Llineval, icoma
1886
1887!!!!!! Variables
1888! datasetn: name of the data set to use
1889! tablen: name of the table to read
1890! Nvar: number of variations (g.e. seasons)
1891! Nr: table row number
1892! Nc: table column number
1893! datasetline: line in the table where the dataset starts
1894
1895    fname='table_characteristics'
1896    errmsg='ERROR -- error -- ERROR -- error'
1897
1898    INQUIRE(FILE=TRIM(tablen), EXIST=existsfile)
1899
1900    IF (existsfile) THEN
1901      DO funit=10,100
1902        INQUIRE(UNIT=funit, OPENED=is_used)
1903        IF (.NOT. is_used) EXIT
1904      END DO
1905
1906    ELSE
1907      PRINT *,TRIM(errmsg)
1908      PRINT *,TRIM(fname),': file "'// TRIM(tablen) // '" does not exist !!!'
1909      STOP
1910    END IF
1911
1912    OPEN(UNIT=funit, FILE=TRIM(tablen), STATUS='OLD')
1913
1914! Looking for the place where data-set starts
1915!!
1916    datasetline = -1
1917    iline=1
1918    DO il=1,100000
1919      READ(funit,*,END=100)lineval
1920      IF (TRIM(lineval) == TRIM(datasetn)) datasetline = il
1921    END DO
1922
1923100 CONTINUE
1924    IF (datasetline == -1) THEN
1925      PRINT *,TRIM(errmsg)
1926      PRINT *,TRIM(fname),': in file "'// TRIM(tablen) // '" dataset: "'//           &
1927        TRIM(datasetn) // '" does not exist !!!'
1928      STOP
1929    END IF
1930
1931! Looking for number of values
1932!!
1933    REWIND(UNIT=funit)
1934    DO ic=1,datasetline
1935      READ(funit,*)lineval
1936    END DO
1937    READ(funit,*)Nr,Nvar
1938
1939    IF (Nvar > 1) READ(funit,*)lineval
1940
1941    READ(funit,'(200A)')lineval
1942    Llineval=LEN(TRIM(lineval))
1943
1944!    PRINT *,'  lineval: ******' // TRIM(lineval) // '******', Llineval
1945
1946    Nc = 0
1947    ic = 1
1948    DO WHILE (ic <= Llineval)
1949!      PRINT *,lineval(ic:Llineval)
1950      icoma = SCAN(lineval(ic:Llineval),',')
1951      IF (icoma > 1) THEN
1952        Nc = Nc + 1
1953!        PRINT *,'  Nc: ',Nc,' ic ',ic,' icoma: ',icoma
1954      ic = ic + icoma
1955      ELSE
1956        ic = ic + Llineval
1957      END IF
1958    END DO
1959    PRINT *,'  dataset :"' // TRIM(datasetn) // '" at line ',datasetline,            &
1960      ' in table :"'// TRIM(tablen) // '" has N variations: ',Nvar,                  &
1961      ' Number of types: ',Nr, ' Number of values: ',Nc
1962
1963    CLOSE(funit)
1964
1965    RETURN
1966
1967  END SUBROUTINE table_characteristics
1968
1969  SUBROUTINE read_table(datasetn,tablen,datasetline,Nvar,Nr,Nc,vals,valns)
1970! Subroutine to read values from a WRF .TBL
1971
1972    IMPLICIT NONE
1973
1974    CHARACTER(LEN=50), INTENT(IN)                        :: datasetn
1975    CHARACTER(LEN=100), INTENT(IN)                       :: tablen
1976    INTEGER, INTENT(IN)                                  :: datasetline, Nvar, Nr, Nc
1977    REAL, DIMENSION(Nvar,Nr,Nc), INTENT(OUT)             :: vals
1978    CHARACTER(LEN=50), DIMENSION(Nr), INTENT(OUT)        :: valns
1979
1980
1981! Local
1982    INTEGER                                              :: iv, ir, ic
1983    CHARACTER(LEN=50)                                    :: fname, errmsg, line
1984    LOGICAL                                              :: existsfile, is_used
1985    INTEGER                                              :: funit
1986
1987!!!!!!! Variables
1988! datasetn: name of the data set to use
1989! tablen: name of the table to read
1990! Nvar: number of variations (g.e. seasons)
1991! Nr: table row number
1992! Nc: table column number
1993! datasetline: line in the table where the dataset starts
1994! vals: values of the data-set
1995! valns: assigned name of each value in the dataset
1996
1997    fname='read_table'
1998    errmsg='ERROR -- error -- ERROR -- error'
1999
2000    INQUIRE(FILE=TRIM(tablen), EXIST=existsfile)
2001
2002    IF (existsfile) THEN
2003      DO funit=10,100
2004        INQUIRE(UNIT=funit, OPENED=is_used)
2005        IF (.NOT. is_used) EXIT
2006      END DO
2007
2008    ELSE
2009      PRINT *,TRIM(errmsg)
2010        PRINT *,TRIM(fname),': file "'// TRIM(tablen) // '" does not exist !!!'
2011        STOP
2012    END IF
2013
2014    OPEN(UNIT=funit, FILE=TRIM(tablen), STATUS='OLD')
2015
2016    DO ic=1,datasetline+1
2017      READ(funit,*)line
2018    END DO
2019
2020    DO iv=1, Nvar
2021      IF (Nvar > 1) READ(funit,*)line
2022      DO ir=1, Nr
2023        READ(funit,*)(vals(iv,ir,ic), ic=1,Nc),valns(ir)
2024      END DO
2025    END DO
2026
2027!    PRINT *,'  values dataset :"' // TRIM(datasetn) // '" in table :"' //TRIM(tablen)
2028!    DO iv=1,Nvar
2029!      PRINT *,'   variation: ',iv
2030!      DO ir=1,Nr
2031!        PRINT *,'    ',(vals(iv,ir,ic), ic=1,Nc),valns(ir)
2032!      END DO
2033!    END DO
2034
2035    CLOSE(funit)
2036
2037    RETURN
2038
2039  END SUBROUTINE read_table
2040
2041
2042  SUBROUTINE compute_soil_mass(Nvar,Nr,Nc,vals,valns,soilt,Nnost,NOsoilt,smois,ddx,  &
2043    ddy,ddz,soilm,soilq)
2044! Subroutine to compute the total soil mass (kg) from a point with different
2045!   percentages of soil types
2046
2047    IMPLICIT NONE
2048
2049    INTEGER, INTENT(IN)                                  :: Nvar, Nr, Nc, Nnost
2050    REAL, DIMENSION(Nvar,Nr,Nc), INTENT(IN)              :: vals
2051    CHARACTER(LEN=50), DIMENSION(Nr), INTENT(IN)         :: valns
2052    REAL, DIMENSION(Nr), INTENT(IN)                      :: soilt
2053    CHARACTER(LEN=50), DIMENSION(Nnost), INTENT(IN)      :: NOsoilt
2054    REAL, INTENT(IN)                                     :: smois, ddx, ddy, ddz
2055    REAL, INTENT(OUT)                                    :: soilm,soilq
2056
2057! Local
2058    INTEGER                                              :: iv, ir, ic
2059    CHARACTER(LEN=50)                                    :: fname, errmsg
2060    LOGICAL                                              :: is_land
2061
2062!!!!!!! Variables
2063! Nvar: number of variations (g.e. seasons)
2064! Nr: table row number
2065! Nc: table column number
2066! vals: values of the data-set (col2: density [g cm-3])
2067! valns: assigned name of each value in the dataset
2068! soilt: vector with the percentages of soil types
2069! Nnost: number of soil types which are not land
2070! NOsoilt: labels of the non-land soil types
2071! smois: soil moisture [m3 m-3]
2072! ddx,ddy,ddz: volume of the grid cell [m]
2073! soilm: total soil mass [kg]
2074! soilq: total soil water [kg kg-1]
2075
2076    fname='read_table'
2077    errmsg='ERROR -- error -- ERROR -- error'
2078
2079    soilm = 0.
2080    soilq = 0.
2081
2082!    PRINT *,'  soilmoisture: ',smois
2083    DO ir=1,Nr
2084      is_land = .TRUE.
2085      IF (soilt(ir) /= 0.) THEN
2086        DO iv=1, Nnost
2087          IF (TRIM(valns(ir)) == TRIM(NOsoilt(iv))) is_land = .FALSE.
2088        END DO
2089        IF (is_land) THEN
2090!          PRINT *,'  soil moisture range: ', vals(1,ir,3), vals(1,ir,5)
2091          soilm = soilm + soilt(ir)*vals(1,ir,2)*ddx*ddy*ddz*1000.
2092          soilq = soilq + soilt(ir)*vals(1,ir,2)*smois*ddx*ddy*ddz*1000.
2093!          PRINT *,'  Percentage of soil "' // TRIM(valns(ir)) //'": ',soilt(ir)*100.,&
2094!            ' % density: ',vals(1,ir,2),' partial mass: ',                           &
2095!            soilt(ir)*vals(1,ir,2)*ddx*ddy*ddz*1000. , ' partial water: ',           &
2096!            soilt(ir)*vals(1,ir,2)*smois*ddx*ddy*ddz*1000.
2097        END IF
2098      END IF
2099    END DO
2100! Remember 1 kg H20 = 1 l = 1mm / m2
2101!!  PRINT *,'  TOTAL. soil mass: ',soilm,' [kg] soil water: ',soilq,' [kg]',           &
2102!!    ' equiv. water depth in soil ',soilq/(ddx*ddy*1000.),' [m]'
2103
2104  END SUBROUTINE compute_soil_mass
2105
2106  SUBROUTINE compute_TOTsoil_mass_water(datname, tabname, Nsoiltypes, soiltypest,    &
2107    soiltypesb, Nsoillayers, smois, soillayers, sizeX, sizeY, soilTOTmass,           &
2108    soilTOThumidity)
2109! Subroutine to compute total soil mass and water from a given grid point using a
2110!   given data-set from SOILPARM.TBL
2111
2112      IMPLICIT NONE
2113
2114      CHARACTER(LEN=50), INTENT(IN)                      :: datname
2115      CHARACTER(LEN=100), INTENT(IN)                     :: tabname
2116      INTEGER, INTENT(IN)                                :: Nsoiltypes, Nsoillayers
2117      REAL, DIMENSION(Nsoiltypes), INTENT(IN)            :: soiltypest, soiltypesb
2118      REAL, DIMENSION(Nsoillayers), INTENT(IN)           :: smois, soillayers
2119      REAL, INTENT(IN)                                   :: sizeX, sizeY
2120      REAL, INTENT(OUT)                                  :: soilTOTmass,             &
2121         soilTOThumidity
2122
2123! Local
2124      INTEGER                                            :: ix,iy,iz
2125      REAL, ALLOCATABLE, DIMENSION(:,:,:)                :: values
2126      CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:)       :: namevalues,NOsoiltypes
2127      INTEGER                                            :: linedataset,Nvariations, &
2128        Nrows, Ncols
2129      INTEGER                                            :: NNOsoiltypes
2130      REAL, ALLOCATABLE, DIMENSION(:)                    :: soiltypes
2131      REAL                                               :: dz,soilmass,soilhumidity
2132      CHARACTER(LEN=50)                                  :: fname, errmsg
2133
2134!!!!!!! Variables
2135! datname: name of the data set to use
2136! tabname: name of the table to read
2137! values: mmatrix with the name of the values: variation (Nvariations, g.e. seasons),
2138!   table row (Nrows),  table column (Ncols)
2139! linedataset: line inside data-set where values start
2140! namevalues: string assigned to each table row
2141! Nsoiltypes: number of soil types
2142! soiltypes[t/b]: percentage of the soil types at a given grid point on top/bottom
2143! soiltypes: percentage of the soil types at a given grid point (lineray interpolated)
2144! NOsoiltypes: vector with the names of soil types which are not land
2145
2146      fname='"compute_TOTsoil_mass_water"'
2147      errmsg='ERROR -- error -- ERROR -- error'
2148
2149      CALL table_characteristics(datname, tabname, linedataset, Nvariations,         &
2150        Nrows, Ncols)
2151
2152      IF (Nrows /= Nsoiltypes) THEN
2153        PRINT *,errmsg
2154        PRINT '  ' // TRIM(fname) // ': Number of soil-types: ',Nsoiltypes,          &
2155          ' does not concide with the number of values: ', Nrows,                    &
2156          ' from the data-set "' // TRIM(datname) // '" !!!'
2157        STOP
2158      END IF
2159
2160      IF (ALLOCATED(values)) DEALLOCATE(values)
2161      ALLOCATE(values(Nvariations,Nrows,Ncols))
2162
2163      IF (ALLOCATED(namevalues)) DEALLOCATE(namevalues)
2164      ALLOCATE(namevalues(Nrows))
2165
2166      CALL read_table(datname, tabname, linedataset, Nvariations, Nrows, Ncols,      &
2167        values, namevalues)
2168
2169      IF (ALLOCATED(soiltypes)) DEALLOCATE(soiltypes)
2170      ALLOCATE(soiltypes(Nrows))
2171
2172! On WRF3.3 SOILPARM.TBL
2173      SELECT CASE (TRIM(datname))
2174
2175        CASE ('STAS')
2176          NNOsoiltypes = 1
2177          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2178          ALLOCATE(NOsoiltypes(NNOsoiltypes))
2179
2180          NOsoiltypes=(/ 'WATER' /)
2181        CASE ('STAS-RUC')
2182          NNOsoiltypes = 1
2183          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2184          ALLOCATE(NOsoiltypes( NNOsoiltypes))
2185
2186          NOsoiltypes=(/ 'WATER' /)
2187      END SELECT
2188
2189      soilTOTmass = 0.
2190      soilTOThumidity = 0.
2191
2192      DO iz=1,Nsoillayers
2193        IF (iz == 1 ) THEN
2194          dz = soillayers(iz)
2195        ELSE
2196          dz = soillayers(iz) - soillayers(iz - 1)
2197        END IF
2198        IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN
2199          DO ix=1,Nrows
2200            CALL interpolate_lin(soillayers(1), soiltypest(ix),                      &
2201              soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz), soiltypes(ix))
2202          END DO
2203        ELSEIF (iz == 1) THEN
2204          soiltypes = soiltypest
2205        ELSE
2206          soiltypes = soiltypesb
2207        END IF
2208!        PRINT *,'lev: ',iz,' dz: ',dz,' soil types: ' ,soiltypes
2209
2210        CALL compute_soil_mass(Nvariations,Nrows,Ncols,values,namevalues,soiltypes,  &
2211          NNOsoiltypes,NOsoiltypes,smois(iz),1.,1.,dz,soilmass,soilhumidity)
2212        soilTOTmass = soilTOTmass + soilmass
2213        soilTOThumidity = soilTOThumidity + soilhumidity
2214      END DO
2215
2216      DEALLOCATE(values)
2217      DEALLOCATE(namevalues)
2218      DEALLOCATE(soiltypes)
2219      DEALLOCATE(NOsoiltypes)
2220
2221  END SUBROUTINE compute_TOTsoil_mass_water
2222
2223  SUBROUTINE compute_TOTsoil_mass_water_values(datname, Nvariations, Nrows, Ncols,   &
2224    values, namevalues, ip, jp, Nsoiltypes, soiltypest, soiltypesb, Nsoillayers,     &
2225    smois, soillayers, sizeX, sizeY, soilTOTmass, soilTOThumidity)
2226! Subroutine to compute total soil mass and water from a given grid point using the
2227!   values from a given data-set of SOILPARM.TBL
2228
2229      IMPLICIT NONE
2230
2231      CHARACTER(LEN=50), INTENT(IN)                      :: datname
2232      INTEGER, INTENT(IN)                                :: Nvariations, Nrows, Ncols
2233      REAL, DIMENSION(Nvariations, Nrows, Ncols), INTENT(IN) :: values
2234      CHARACTER(LEN=50), DIMENSION(Nrows), INTENT(IN)    :: namevalues
2235      INTEGER, INTENT(IN)                                :: ip, jp
2236      INTEGER, INTENT(IN)                                :: Nsoiltypes, Nsoillayers
2237      REAL, DIMENSION(Nsoiltypes), INTENT(IN)            :: soiltypest, soiltypesb
2238      REAL, DIMENSION(Nsoillayers), INTENT(IN)           :: smois, soillayers
2239      REAL, INTENT(IN)                                   :: sizeX, sizeY
2240      REAL, INTENT(OUT)                                  :: soilTOTmass,             &
2241         soilTOThumidity
2242
2243! Local
2244      INTEGER                                            :: ix,iy,iz
2245      CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:)       :: NOsoiltypes
2246      INTEGER                                            :: linedataset
2247      INTEGER                                            :: NNOsoiltypes
2248      REAL, ALLOCATABLE, DIMENSION(:)                    :: soiltypes
2249      REAL                                               :: dz,soilmass,soilhumidity
2250      CHARACTER(LEN=50)                                  :: fname, errmsg
2251
2252!!!!!!! Variables
2253! datname: name of the data set to use
2254! values: mmatrix with the name of the values: variation (Nvariations, g.e. seasons),
2255!   table row (Nrows),  table column (Ncols)
2256! linedataset: line inside data-set where values start
2257! namevalues: string assigned to each table row
2258! ip, jp: coordinates of the grid point
2259! Nsoiltypes: number of soil types
2260! soiltypes[t/b]: percentage of the soil types at a given grid point on top/bottom
2261! soiltypes: percentage of the soil types at a given grid point (lineray interpolated)
2262! NOsoiltypes: vector with the names of soil types which are not land
2263
2264      fname='"compute_TOTsoil_mass_water_values"'
2265      errmsg='ERROR -- error -- ERROR -- error'
2266
2267      IF (Nrows /= Nsoiltypes) THEN
2268        PRINT *,errmsg
2269        PRINT '  ' // TRIM(fname) // ': Number of soil-types: ',Nsoiltypes,          &
2270          ' does not concide with the number of values: ', Nrows,                    &
2271          ' from the data-set "' // TRIM(datname) // '" !!!'
2272        STOP
2273      END IF
2274
2275      IF (ALLOCATED(soiltypes)) DEALLOCATE(soiltypes)
2276      ALLOCATE(soiltypes(Nrows))
2277
2278! On WRF3.3 SOILPARM.TBL
2279      SELECT CASE (TRIM(datname))
2280
2281        CASE ('STAS')
2282          NNOsoiltypes = 1
2283          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2284          ALLOCATE(NOsoiltypes(NNOsoiltypes))
2285
2286          NOsoiltypes=(/ 'WATER' /)
2287        CASE ('STAS-RUC')
2288          NNOsoiltypes = 1
2289          IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes)
2290          ALLOCATE(NOsoiltypes( NNOsoiltypes))
2291
2292          NOsoiltypes=(/ 'WATER' /)
2293      END SELECT
2294
2295      soilTOTmass = 0.
2296      soilTOThumidity = 0.
2297
2298      DO iz=1,Nsoillayers
2299        IF (iz == 1 ) THEN
2300          dz = soillayers(iz)
2301        ELSE
2302          dz = soillayers(iz) - soillayers(iz - 1)
2303        END IF
2304        IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN
2305          DO ix=1,Nrows
2306            CALL interpolate_lin(soillayers(1), soiltypest(ix),                      &
2307              soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz), soiltypes(ix))
2308          END DO
2309        ELSEIF (iz == 1) THEN
2310          soiltypes = soiltypest
2311        ELSE
2312          soiltypes = soiltypesb
2313        END IF
2314!        PRINT *,'lev: ',iz,' dz: ',dz,' soil types: ' ,soiltypes
2315
2316        CALL compute_soil_mass(Nvariations,Nrows,Ncols,values,namevalues,soiltypes,  &
2317          NNOsoiltypes,NOsoiltypes,smois(iz),sizeX,sizeY,dz,soilmass,soilhumidity)
2318        soilTOTmass = soilTOTmass + soilmass
2319        soilTOThumidity = soilTOThumidity + soilhumidity
2320      END DO
2321
2322      IF ( (soilTOTmass < 0.) .OR. (soilTOThumidity < 0.) ) THEN
2323        PRINT *,errmsg
2324        PRINT *,'  ' // TRIM(fname) // ' at point ',ip, ', ', jp,                    &
2325          ' negative total soil mass: ',soilTOTmass, ' or soil water content',       &
2326          soilTOThumidity
2327        PRINT *,'  Soil_name density soil_types_________'
2328        DO ix=1, Nrows
2329           PRINT *,ix, TRIM(namevalues(ix)), values(1,ix,2), soiltypes(ix)
2330        END DO
2331        PRINT *,'  depth_layer partial_mass soil+moisure partial_moisture__________'
2332        DO iz=1,Nsoillayers
2333          IF (iz == 1 ) THEN
2334            dz = soillayers(iz)
2335          ELSE
2336            dz = soillayers(iz) - soillayers(iz - 1)
2337          END IF
2338          IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN
2339            DO ix=1,Nrows
2340              CALL interpolate_lin(soillayers(1), soiltypest(ix),                    &
2341                soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz),soiltypes(ix))
2342            END DO
2343          ELSEIF (iz == 1) THEN
2344            soiltypes = soiltypest
2345          ELSE
2346            soiltypes = soiltypesb
2347          END IF
2348          PRINT *,'  ',iz, dz, values(1,iz,2)*sizeX*sizeY*dz*1000., smois(iz),       &
2349             smois(iz)*values(1,iz,2)*sizeX*sizeY*dz*1000.
2350        END DO
2351        STOP
2352      END IF
2353
2354
2355      DEALLOCATE(soiltypes)
2356      DEALLOCATE(NOsoiltypes)
2357
2358  END SUBROUTINE compute_TOTsoil_mass_water_values
2359
2360! cat Registry_out.LMDZ | awk '{print  "        pwrf_grid%"tolower($3)" "$4" = p"tolower($3)" "$4}' >> get_lmdz_out.f90
2361! cat get_lmdz_out.f90n' ',' '{print $3}' | tr '\n
2362  SUBROUTINE get_lmdz_out(ddx, ddy, ddz, bdy, dlmdz, Nks, Nz0, nqtot, ivap, iliq,    &
2363    iflag_pbl, iflag_con, iflag_thermals, iflag_wake, read_climoz, itap,             &
2364    pdtphys, paprs, pplay, pphi, pphis, presnivs, qx,                                &
2365    d_u, d_v, d_t, d_qx, d_ps,                                                       &
2366    pwrf_grid)
2367! Subroutine to get all the LMDZ output variables into the WRF Registry structure
2368
2369! WRF modules
2370    USE module_domain_type
2371
2372! LMDZ modules
2373    USE indice_sol_mod
2374    USE phys_state_var_mod
2375    USE phys_output_var_mod
2376!    USE pbl_surface_mod
2377    USE phys_local_var_mod
2378    USE comgeomphy
2379
2380! WRF_LMDZ modules
2381    USE output_lmdz_NOmodule
2382
2383! netcdf
2384    USE netcdf, ONLY: nf90_fill_real
2385
2386    IMPLICIT NONE
2387
2388    TYPE(domain) , TARGET                                :: pwrf_grid
2389    INTEGER, INTENT(IN)                                  :: ddx, ddy, ddz, bdy,      &
2390      dlmdz, Nks, Nz0, nqtot, ivap, iliq, iflag_pbl, iflag_con, iflag_thermals,      &
2391      iflag_wake, read_climoz, itap
2392    REAL, INTENT(IN)                                     :: pdtphys
2393    REAL, DIMENSION(dlmdz), INTENT(IN)                   :: pphis, d_ps
2394    REAL, DIMENSION(ddz), INTENT(IN)                     :: presnivs
2395    REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN)             :: paprs
2396    REAL, DIMENSION(dlmdz,ddz), INTENT(IN)               :: pplay, pphi, d_u, d_v, d_t
2397    REAL, DIMENSION(dlmdz,ddz,nqtot), INTENT(IN)         :: qx, d_qx
2398
2399!    EXTERNAL suphel
2400
2401! Local
2402    INTEGER                                              :: dx, dy, ix, iy, ixy, iz
2403    INTEGER                                              :: i,j,k,l,n
2404    INTEGER                                              :: nsrf
2405    REAL, DIMENSION(dlmdz,ddz)                           :: var3d
2406    REAL, DIMENSION(dlmdz,ddz+1)                         :: var3dstagg
2407    REAL                                                 :: RDAY, RG, RD, RMO3
2408    CHARACTER(LEN=50)                                    :: fname, errmsg
2409    CHARACTER(LEN=4)                                     :: bb2
2410    REAL, DIMENSION(dlmdz)                               :: zx_tmp_fi2d ! variable temporaire grille physique
2411    REAL, DIMENSION(dlmdz,ddz)                           :: zx_tmp_fi3d ! variable temporaire pour champs 3D
2412    REAL, DIMENSION(dlmdz,ddz+1)                         :: zx_tmp_fi3d1 !variable temporaire pour champs 3D (kelvp1)
2413    REAL(KIND=8), DIMENSION(dlmdz,ddz)                   :: zx_tmp2_fi3d ! variable temporaire pour champs 3D
2414    REAL, PARAMETER                                      :: dobson_u = 2.1415e-05
2415
2416!    klon=dldmz
2417!    klev=ddz
2418
2419
2420!! To think about / discuss with Frederic
2421
2422    INCLUDE "declare_STDlev.h"
2423!    INCLUDE "calcul_STDlev.h"
2424
2425!!!!!!! Variables
2426! pwrf_grid: WRF grid structure
2427
2428! L. Fita, LMD January 2014. It should work using suphel, but not?
2429!    CALL suphel
2430
2431    fname='"get_lmdz_out"'
2432    errmsg='ERROR -- error -- ERROR -- error'
2433
2434    RDAY = 86400.
2435    RG = 9.81
2436    RD = 287.
2437    RMO3=47.9942
2438
2439! Remove variable assignation already done in get_lmdz? These variables are not used
2440!    in restart/boundaries
2441
2442    DO ix=1,ddx
2443      DO iy=1,ddy
2444        ixy=ddx*(iy-1)+ix
2445
2446        pwrf_grid%lages_ter(ix,iy) = agesno(ixy,1)
2447        pwrf_grid%lages_lic(ix,iy) = agesno(ixy,2)
2448        pwrf_grid%lages_oce(ix,iy) = agesno(ixy,3)
2449        pwrf_grid%lages_sic(ix,iy) = agesno(ixy,4)
2450!!        pwrf_grid%lahyb(misc??) = pahyb misc
2451        pwrf_grid%laire(ix,iy) = airephy(ixy)
2452        pwrf_grid%laireter(ix,iy) = paire_ter(ixy)
2453        pwrf_grid%lalb1(ix,iy) = albsol1(ixy)
2454        pwrf_grid%lalb2(ix,iy) = albsol2(ixy)
2455        pwrf_grid%lalbe_ter(ix,iy) = falb1(ixy,1)
2456        pwrf_grid%lalbe_lic(ix,iy) = falb1(ixy,2)
2457        pwrf_grid%lalbe_oce(ix,iy) = falb1(ixy,3)
2458        pwrf_grid%lalbe_sic(ix,iy) = falb1(ixy,4)
2459        pwrf_grid%lale(ix,iy) = ale(ixy)
2460        pwrf_grid%lale_bl(ix,iy) = ale_bl(ixy)
2461        pwrf_grid%lale_wk(ix,iy) = ale_wake(ixy)
2462        pwrf_grid%lalp(ix,iy) = alp(ixy)
2463        pwrf_grid%lalp_bl(ix,iy) = alp_bl(ixy)
2464        pwrf_grid%lalp_wk(ix,iy) = alp_wake(ixy)
2465!!        pwrf_grid%lalt(misc??) = palt misc
2466!!        pwrf_grid%lbhyb(misc??) = pbhyb misc
2467        pwrf_grid%lbils(ix,iy) = bils(ixy)
2468        pwrf_grid%lbils_diss(ix,iy) = bils_diss(ixy)
2469        pwrf_grid%lbils_ec(ix,iy) = bils_ec(ixy)
2470        pwrf_grid%lbils_enthalp(ix,iy) = bils_enthalp(ixy)
2471        pwrf_grid%lbils_kinetic(ix,iy) = bils_kinetic(ixy)
2472        pwrf_grid%lbils_latent(ix,iy) = bils_latent(ixy)
2473        pwrf_grid%lbils_tke(ix,iy) = bils_tke(ixy)
2474        pwrf_grid%lcape(ix,iy) = cape(ixy)
2475        pwrf_grid%lcape_max(ix,iy) = cape(ixy)
2476        pwrf_grid%lcdrh(ix,iy) = cdragh(ixy)
2477        pwrf_grid%lcdrm(ix,iy) = cdragm(ixy)
2478        pwrf_grid%lcldh(ix,iy) = cldh(ixy)
2479        pwrf_grid%lcldl(ix,iy) = cldl(ixy)
2480        pwrf_grid%lcldm(ix,iy) = cldm(ixy)
2481        pwrf_grid%lcldq(ix,iy) = cldq(ixy)
2482        pwrf_grid%lcldt(ix,iy) = cldt(ixy)
2483        pwrf_grid%lcontfracatm(ix,iy) = pctsrf(ixy,is_ter)+pctsrf(ixy,is_lic)
2484        pwrf_grid%lcontfracor(ix,iy) = pctsrf(ixy,is_ter)
2485!NOTknown!        pwrf_grid%lcumpb(ix,iy) = cumpb(ixy)
2486!NOTknown!        pwrf_grid%lcumrn(ix,iy) = cumrn(ixy)
2487        pwrf_grid%ldthmin(ix,iy) = dthmin(ixy)
2488        pwrf_grid%ldtsvdft(ix,iy) = d_ts(ixy,is_ter)
2489        pwrf_grid%ldtsvdfo(ix,iy) = d_ts(ixy,is_oce)
2490        pwrf_grid%ldtsvdfg(ix,iy) = d_ts(ixy,is_lic)
2491        pwrf_grid%ldtsvdfi(ix,iy) = d_ts(ixy,is_sic)
2492! L. Fita, October 2014. with fevap there is not evaporation!
2493!        pwrf_grid%levap(ix,iy) = fevap(ixy,1)
2494        pwrf_grid%levap(ix,iy) = evap(ixy)
2495        pwrf_grid%levap_ter(ix,iy) = fevap(ixy,1)
2496        pwrf_grid%levap_lic(ix,iy) = fevap(ixy,2)
2497        pwrf_grid%levap_oce(ix,iy) = fevap(ixy,3)
2498        pwrf_grid%levap_sic(ix,iy) = fevap(ixy,4)
2499        pwrf_grid%levappot_ter(ix,iy) = evap_pot(ixy,1)
2500        pwrf_grid%levappot_lic(ix,iy) = evap_pot(ixy,2)
2501        pwrf_grid%levappot_oce(ix,iy) = evap_pot(ixy,3)
2502        pwrf_grid%levappot_sic(ix,iy) = evap_pot(ixy,4)
2503        pwrf_grid%lfbase(ix,iy) = ema_cbmf(ixy)
2504        pwrf_grid%lfder(ix,iy) = fder(ixy)
2505        pwrf_grid%lffonte(ix,iy) = zxffonte(ixy)
2506        pwrf_grid%lflat(ix,iy) = zxfluxlat(ixy)
2507        pwrf_grid%lflw_ter(ix,iy) = fsollw(ixy,1)
2508        pwrf_grid%lflw_lic(ix,iy) = fsollw(ixy,2)
2509        pwrf_grid%lflw_oce(ix,iy) = fsollw(ixy,3)
2510        pwrf_grid%lflw_sic(ix,iy) = fsollw(ixy,4)
2511        pwrf_grid%lfqcalving(ix,iy) = zxfqcalving(ixy)
2512        pwrf_grid%lfqfonte(ix,iy) = zxfqfonte(ixy)
2513        pwrf_grid%lfract_ter(ix,iy) = pctsrf(ixy,1)
2514        pwrf_grid%lfract_lic(ix,iy) = pctsrf(ixy,2)
2515        pwrf_grid%lfract_oce(ix,iy) = pctsrf(ixy,3)
2516        pwrf_grid%lfract_sic(ix,iy) = pctsrf(ixy,4)
2517        pwrf_grid%lfsnow(ix,iy) = zfra_o(ixy)
2518        pwrf_grid%lfsw_ter(ix,iy) = fsolw(ixy,1)
2519        pwrf_grid%lfsw_lic(ix,iy) = fsolw(ixy,2)
2520        pwrf_grid%lfsw_oce(ix,iy) = fsolw(ixy,3)
2521        pwrf_grid%lfsw_sic(ix,iy) = fsolw(ixy,4)
2522!NOtnow!        pwrf_grid%lftime_con(ix,iy) = float(itau_con)/float(itap)
2523        pwrf_grid%lftime_th(ix,iy) = 0.
2524        pwrf_grid%liwp(ix,iy) = fiwp(ixy)
2525!!        pwrf_grid%llat j = pllat j
2526        pwrf_grid%llat_ter(ix,iy) = fluxlat(ixy,1)
2527        pwrf_grid%llat_lic(ix,iy) = fluxlat(ixy,2)
2528        pwrf_grid%llat_oce(ix,iy) = fluxlat(ixy,3)
2529        pwrf_grid%llat_sic(ix,iy) = fluxlat(ixy,4)
2530        pwrf_grid%llmaxth(ix,iy) = lmax_th(ixy)
2531!!        pwrf_grid%llon i = pllon i
2532        pwrf_grid%llwdn200(ix,iy) = lwdn200(ixy)
2533        pwrf_grid%llwdn200clr(ix,iy) = lwdn200clr(ixy)
2534        pwrf_grid%llwdnsfc(ix,iy) = sollwdown(ixy)
2535        pwrf_grid%llwdnsfcclr(ix,iy) = sollwdownclr(ixy)
2536        pwrf_grid%llwdownor(ix,iy) = sollwdown(ixy)
2537        pwrf_grid%llwp(ix,iy) = flwp(ixy)
2538        pwrf_grid%llwup200(ix,iy) = lwup200(ixy)
2539        pwrf_grid%llwup200clr(ix,iy) = lwup200clr(ixy)
2540        pwrf_grid%llwupsfc(ix,iy) = sollwdown(ixy)-sollw(ixy)
2541        pwrf_grid%llwupsfcclr(ix,iy) = -1.*lwdn0(ixy,1)-sollw0(ixy)
2542        pwrf_grid%lmsnow(ix,iy) = snow_o(ixy)
2543        pwrf_grid%lndayrain(ix,iy) = nday_rain(ixy)
2544        pwrf_grid%lnettop(ix,iy) = topsw(ixy)-toplw(ixy)
2545        pwrf_grid%lpbase(ix,iy) = ema_pcb(ixy)
2546        pwrf_grid%lphis(ix,iy) = pphis(ixy)
2547        pwrf_grid%lplcl(ix,iy) = plcl(ixy)
2548        pwrf_grid%lplfc(ix,iy) = plfc(ixy)
2549        pwrf_grid%lpluc(ix,iy) = rain_con(ixy) + snow_con(ixy)
2550        pwrf_grid%lplul(ix,iy) = rain_lsc(ixy) + snow_lsc(ixy)
2551        pwrf_grid%lplulst(ix,iy) = plul_st(ixy)
2552        pwrf_grid%lplulth(ix,iy) = plul_th(ixy)
2553        pwrf_grid%lpourc_ter(ix,iy) = pctsrf(ixy,1)*100.
2554        pwrf_grid%lpourc_lic(ix,iy) = pctsrf(ixy,2)*100.
2555        pwrf_grid%lpourc_oce(ix,iy) = pctsrf(ixy,3)*100.
2556        pwrf_grid%lpourc_sic(ix,iy) = pctsrf(ixy,4)*100.
2557        pwrf_grid%lprecip(ix,iy) = rain_fall(ixy) + snow_fall(ixy)
2558!!        pwrf_grid%lpresnivs k = plpresnivs k
2559        pwrf_grid%lprw(ix,iy) = prw(ixy)
2560        pwrf_grid%lpsol(ix,iy) = paprs(ixy,1)
2561        pwrf_grid%lptop(ix,iy) = ema_pct(ixy)
2562        pwrf_grid%lq2m(ix,iy) = zq2m(ixy)
2563        pwrf_grid%lqsat2m(ix,iy) = qsat2m(ixy)
2564        pwrf_grid%lqsol(ix,iy) = qsol(ixy)
2565        pwrf_grid%lqsurf(ix,iy) = zxqsurf(ixy)
2566        pwrf_grid%lradsol(ix,iy) = radsol(ixy)
2567        pwrf_grid%lrh2m(ix,iy) = MIN(100.,rh2m(ixy)*100.)
2568        pwrf_grid%lrh2m_max(ix,iy) = MIN(100.,rh2m(ixy)*100.)
2569        pwrf_grid%lrh2m_min(ix,iy) = MIN(100.,rh2m(ixy)*100.)
2570        pwrf_grid%lrugs(ix,iy) = zxrugs(ixy)
2571        pwrf_grid%lrugs_ter(ix,iy) = frugs(ixy,1)
2572        pwrf_grid%lrugs_lic(ix,iy) = frugs(ixy,2)
2573        pwrf_grid%lrugs_oce(ix,iy) = frugs(ixy,3)
2574        pwrf_grid%lrugs_sic(ix,iy) = frugs(ixy,4)
2575        pwrf_grid%lsens(ix,iy) = -1.*sens(ixy)
2576        pwrf_grid%lsicf(ix,iy) = pctsrf(ixy,is_sic)
2577        pwrf_grid%ls_lcl(ix,iy) = s_lcl(ixy)
2578        pwrf_grid%lslp(ix,iy) = slp(ixy)
2579        pwrf_grid%lsnow(ix,iy) = snow_fall(ixy)
2580        pwrf_grid%lsnowl(ix,iy) = snow_lsc(ixy)
2581        pwrf_grid%lsoll(ix,iy) = sollw(ixy)
2582        pwrf_grid%lsoll0(ix,iy) = sollw0(ixy)
2583        pwrf_grid%lsollwdown(ix,iy) = sollwdown(ixy)
2584        pwrf_grid%lsols(ix,iy) = solsw(ixy)
2585        pwrf_grid%lsols0(ix,iy) = solsw0(ixy)
2586        pwrf_grid%ls_pblh(ix,iy) = s_pblh(ixy)
2587        pwrf_grid%ls_pblt(ix,iy) = s_pblt(ixy)
2588        pwrf_grid%ls_therm(ix,iy) = s_therm(ixy)
2589        pwrf_grid%lswdn200(ix,iy) = swdn200(ixy)
2590        pwrf_grid%lswdn200clr(ix,iy) = swdn200clr(ixy)
2591        pwrf_grid%lswdnsfc(ix,iy) = swdn(ixy,1)
2592        pwrf_grid%lswdnsfcclr(ix,iy) = swdn0(ixy,1)
2593        pwrf_grid%lswdntoa(ix,iy) = swdn(ixy,klev+1)
2594        pwrf_grid%lswdntoaclr(ix,iy) = swdn0(ixy,klev+1)
2595        pwrf_grid%lswdownor(ix,iy) = solsw(ixy)/(1.-albsol1(ixy))
2596        pwrf_grid%lswnetor(ix,iy) = fsolsw(ixy,is_ter)
2597        pwrf_grid%lswup200(ix,iy) = swup200(ixy)
2598        pwrf_grid%lswup200clr(ix,iy) = swup200clr(ixy)
2599        pwrf_grid%lswupsfc(ix,iy) = swup(ixy,1)
2600        pwrf_grid%lswupsfcclr(ix,iy) = swup0(ixy,1)
2601        pwrf_grid%lswuptoa(ix,iy) = swup(ixy,klev+1)
2602        pwrf_grid%lswuptoaclr(ix,iy) = swup0(ixy,klev+1)
2603        pwrf_grid%lt2m(ix,iy) = zt2m(ixy)
2604        pwrf_grid%lt2m_max(ix,iy) = zt2m(ixy)
2605        pwrf_grid%lt2m_min(ix,iy) = zt2m(ixy)
2606        pwrf_grid%lt2m_ter(ix,iy) = t2m(ixy,1)
2607        pwrf_grid%lt2m_lic(ix,iy) = t2m(ixy,2)
2608        pwrf_grid%lt2m_oce(ix,iy) = t2m(ixy,3)
2609        pwrf_grid%lt2m_sic(ix,iy) = t2m(ixy,4)
2610        pwrf_grid%ltaux(ix,iy) = 0.
2611        pwrf_grid%ltauy(ix,iy) = 0.
2612        DO iz=1,4
2613          pwrf_grid%ltaux(ix,iy) = pwrf_grid%ltauy(ix,iy) + pctsrf(ixy,iz)*          &
2614            fluxu(ixy,1,iz)
2615          pwrf_grid%ltauy(ix,iy) = pwrf_grid%ltauy(ix,iy) + pctsrf(ixy,iz)*          &
2616            fluxv(ixy,1,iz)
2617        END DO
2618        pwrf_grid%ltaux_ter(ix,iy) = fluxu(ixy,1,1)
2619        pwrf_grid%ltaux_lic(ix,iy) = fluxu(ixy,1,2)
2620        pwrf_grid%ltaux_oce(ix,iy) = fluxu(ixy,1,3)
2621        pwrf_grid%ltaux_sic(ix,iy) = fluxu(ixy,1,4)
2622        pwrf_grid%ltauy_ter(ix,iy) = fluxv(ixy,1,1)
2623        pwrf_grid%ltauy_lic(ix,iy) = fluxv(ixy,1,2)
2624        pwrf_grid%ltauy_oce(ix,iy) = fluxv(ixy,1,3)
2625        pwrf_grid%ltauy_sic(ix,iy) = fluxv(ixy,1,4)
2626!!        pwrf_grid%ltime - = pltime -
2627!!        pwrf_grid%ltime_bounds {lbnds} = pltime_bounds {lbnds}
2628!!        pwrf_grid%lti86400 - = plti86400 -
2629        IF ( (pctsrf(ixy,is_oce).GT.epsfra).OR.(pctsrf(ixy,is_sic).GT.epsfra) ) THEN
2630          pwrf_grid%lt_oce_sic(ix,iy) = (ftsol(ixy, is_oce) * pctsrf(ixy,is_oce)+    &
2631          ftsol(ixy, is_sic)*pctsrf(ixy,is_sic))/(pctsrf(ixy,is_oce)+                &
2632          pctsrf(ixy,is_sic))
2633        ELSE
2634          pwrf_grid%lt_oce_sic(ix,iy) = 273.15
2635        END IF
2636!!        pwrf_grid%lt_op_00001800(misc??) = pt_op_00001800 misc
2637!!        pwrf_grid%lt_op_00001800_bnds(misc??) = pt_op_00001800_bnds misc
2638        pwrf_grid%ltopl(ix,iy) = toplw(ixy)
2639        pwrf_grid%ltopl0(ix,iy) = toplw0(ixy)
2640        pwrf_grid%ltops(ix,iy) = topsw(ixy)
2641        pwrf_grid%ltops0(ix,iy) = topsw0(ixy)
2642        pwrf_grid%ltpot(ix,iy) = tpot(ixy)
2643        pwrf_grid%ltpote(ix,iy) = tpote(ixy)
2644        pwrf_grid%ltsol(ix,iy) = zxtsol(ixy)
2645        pwrf_grid%ltsol_ter(ix,iy) = ftsol(ixy,1)
2646        pwrf_grid%ltsol_lic(ix,iy) = ftsol(ixy,2)
2647        pwrf_grid%ltsol_oce(ix,iy) = ftsol(ixy,3)
2648        pwrf_grid%ltsol_sic(ix,iy) = ftsol(ixy,4)
2649        pwrf_grid%lu10m(ix,iy) = zu10m(ixy)
2650        pwrf_grid%lu10m_ter(ix,iy) = u10m(ixy,1)
2651        pwrf_grid%lu10m_lic(ix,iy) = u10m(ixy,2)
2652        pwrf_grid%lu10m_oce(ix,iy) = u10m(ixy,3)
2653        pwrf_grid%lu10m_sic(ix,iy) = u10m(ixy,4)
2654        pwrf_grid%lue(ix,iy) = ue(ixy)
2655        pwrf_grid%luq(ix,iy) = uq(ixy)
2656        pwrf_grid%lustar(ix,iy) = zustar(ixy)
2657        pwrf_grid%lustar_ter(ix,iy) = ustar(ixy,1)
2658        pwrf_grid%lustar_lic(ix,iy) = ustar(ixy,2)
2659        pwrf_grid%lustar_oce(ix,iy) = ustar(ixy,3)
2660        pwrf_grid%lustar_sic(ix,iy) = ustar(ixy,4)
2661        pwrf_grid%lv10m(ix,iy) = zv10m(ixy)
2662        pwrf_grid%lv10m_ter(ix,iy) = v10m(ixy,1)
2663        pwrf_grid%lv10m_lic(ix,iy) = v10m(ixy,2)
2664        pwrf_grid%lv10m_oce(ix,iy) = v10m(ixy,3)
2665        pwrf_grid%lv10m_sic(ix,iy) = v10m(ixy,4)
2666        pwrf_grid%lve(ix,iy) = ve(ixy)
2667        pwrf_grid%lvq(ix,iy) = vq(ixy)
2668        pwrf_grid%lwake_h(ix,iy) = wake_h(ixy)
2669        pwrf_grid%lwake_s(ix,iy) = wake_s(ixy)
2670        pwrf_grid%lwape(ix,iy) = wake_pe(ixy)
2671        pwrf_grid%lwbeff(ix,iy) = wbeff(ixy)
2672        pwrf_grid%lwbilo_ter(ix,iy) = wfbilo(ixy,1)
2673        pwrf_grid%lwbilo_lic(ix,iy) = wfbilo(ixy,2)
2674        pwrf_grid%lwbilo_oce(ix,iy) = wfbilo(ixy,3)
2675        pwrf_grid%lwbilo_sic(ix,iy) = wfbilo(ixy,4)
2676        pwrf_grid%lwbils_ter(ix,iy) = wfbils(ixy,1)
2677        pwrf_grid%lwbils_lic(ix,iy) = wfbils(ixy,2)
2678        pwrf_grid%lwbils_oce(ix,iy) = wfbils(ixy,3)
2679        pwrf_grid%lwbils_sic(ix,iy) = wfbils(ixy,4)
2680        pwrf_grid%lweakinv(ix,iy) = weak_inversion(ixy)
2681        pwrf_grid%lwind10m(ix,iy) = SQRT(zu10m(ixy)*zu10m(ixy)+zv10m(ixy)*zv10m(ixy))
2682        pwrf_grid%lwind10max(ix,iy) = SQRT(zu10m(ixy)*zu10m(ixy)+zv10m(ixy)*zv10m(ixy))
2683        pwrf_grid%lzmax_th(ix,iy) = zmax_th(ixy)
2684! Output at the standard levels
2685!!
2686!        DO iz=1, nlevSTD
2687!          bb2=clevSTD(iz)
2688 
2689!          IF (TRIM(bb2) == "850") THEN
2690!            pwrf_grid%lq850(ix,iy) = q850(ixy)
2691!            pwrf_grid%lt850(ix,iy) = t850(ixy)
2692!            pwrf_grid%lu850(ix,iy) = u850(ixy)
2693!            pwrf_grid%lv850(ix,iy) = v850(ixy)
2694!            pwrf_grid%lw850(ix,iy) = w850(ixy)
2695!            pwrf_grid%lz850(ix,iy) = z850(ixy)
2696!          ELSE IF (TRIM(bb2) == "700") THEN
2697!            pwrf_grid%lq700(ix,iy) = q700(ixy)
2698!            pwrf_grid%lt700(ix,iy) = t700(ixy)
2699!            pwrf_grid%lu700(ix,iy) = u700(ixy)
2700!            pwrf_grid%lv700(ix,iy) = v700(ixy)
2701!            pwrf_grid%lw700(ix,iy) = w700(ixy)
2702!            pwrf_grid%lz700(ix,iy) = z700(ixy)
2703!          ELSE IF (TRIM(bb2) == "500") THEN
2704!            pwrf_grid%lq500(ix,iy) = q500(ixy)
2705!            pwrf_grid%lt500(ix,iy) = t500(ixy)
2706!            pwrf_grid%lu500(ix,iy) = u500(ixy)
2707!            pwrf_grid%lv500(ix,iy) = v500(ixy)
2708!            pwrf_grid%lw500(ix,iy) = w500(ixy)
2709!            pwrf_grid%lz500(ix,iy) = z500(ixy)
2710!           ELSE IF (TRIM(bb2) == "200")THEN
2711!            pwrf_grid%lq200(ix,iy) = q200(ixy)
2712!            pwrf_grid%lt200(ix,iy) = t200(ixy)
2713!            pwrf_grid%lu200(ix,iy) = u200(ixy)
2714!            pwrf_grid%lv200(ix,iy) = v200(ixy)
2715!            pwrf_grid%lw200(ix,iy) = w200(ixy)
2716!            pwrf_grid%lz200(ix,iy) = z200(ixy)
2717!          ELSE IF (TRIM(bb2) == "100") THEN
2718!            pwrf_grid%lq100(ix,iy) = q100(ixy)
2719!            pwrf_grid%lt100(ix,iy) = t100(ixy)
2720!            pwrf_grid%lu100(ix,iy) = u100(ixy)
2721!            pwrf_grid%lv100(ix,iy) = v100(ixy)
2722!            pwrf_grid%lw100(ix,iy) = w100(ixy)
2723!            pwrf_grid%lz100(ix,iy) = z100(ixy)
2724!          ELSE IF (TRIM(bb2) == "50") THEN
2725!            pwrf_grid%lq50(ix,iy) = q50(ixy)
2726!            pwrf_grid%lt50(ix,iy) = t50(ixy)
2727!            pwrf_grid%lu50(ix,iy) = u50(ixy)
2728!            pwrf_grid%lv50(ix,iy) = v50(ixy)
2729!            pwrf_grid%lw50(ix,iy) = w50(ixy)
2730!            pwrf_grid%lz50(ix,iy) = z50(ixy)
2731!          ELSE IF (TRIM(bb2) == "10") THEN
2732!            pwrf_grid%lq10(ix,iy) = q10(ixy)
2733!            pwrf_grid%lt10(ix,iy) = t10(ixy)
2734!            pwrf_grid%lu10(ix,iy) = u10(ixy)
2735!            pwrf_grid%lv10(ix,iy) = v10(ixy)
2736!            pwrf_grid%lw10(ix,iy) = w10(ixy)
2737!            pwrf_grid%lz10(ix,iy) = z10(ixy)
2738!          END IF
2739!        END DO
2740
2741      END DO
2742    END DO
2743
2744    dx=ddx+2*bdy
2745    dy=ddy+2*bdy
2746
2747    CALL vect_mat_zstagg(fraca, dx, dy, ddz+1, bdy, pwrf_grid%la_th)
2748    CALL vect_mat(beta_prec, dx, dy, ddz, bdy, pwrf_grid%lbeta_prec)
2749    var3d=upwd + dnwd+ dnwd0
2750    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldmc)
2751    CALL vect_mat(dnwd, dx, dy, ddz, bdy, pwrf_grid%ldnwd)
2752    CALL vect_mat(dnwd0, dx, dy, ddz, bdy, pwrf_grid%ldnwd0)
2753! Originally is var3d = d_q_ajsb(1:ddx*ddy,1:ddz)/pdtphys
2754    var3d = d_q_ajsb/pdtphys
2755    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqajs)
2756    var3d = d_q_con/pdtphys
2757    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqcon)
2758    CALL vect_mat(d_q_dyn, dx, dy, ddz, bdy, pwrf_grid%ldqdyn)
2759    var3d = d_q_eva/pdtphys
2760    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqeva)
2761    var3d = d_q_lsc/pdtphys
2762    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlsc)
2763    var3d = d_q_lscst/pdtphys
2764    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlscst)
2765    var3d = d_q_lscth/pdtphys
2766    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlscth)
2767! L. Fita, LMD January 2014. d_qx is the tendency of moisture due to physiq!
2768!   MOISTtend in the WRF_LMDZ coupling
2769    CALL vect_mat(d_qx(:,:,ivap), dx, dy, ddz, bdy, pwrf_grid%ldqphy)
2770    var3d = d_q_ajs/pdtphys - d_q_ajsb/pdtphys
2771    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqthe)
2772    var3d = d_q_vdf/pdtphys
2773    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqvdf)
2774    var3d = d_q_wake/pdtphys
2775    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqwak)
2776    var3d = d_t_ajsb/pdtphys
2777    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtajs)
2778    var3d = d_t_con/pdtphys
2779    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtcon)
2780    var3d = d_t_diss/pdtphys
2781    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtdis)
2782    CALL vect_mat(d_t_dyn, dx, dy, ddz, bdy, pwrf_grid%ldtdyn)
2783    var3d = d_t_ec/pdtphys
2784    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtec)
2785    var3d = d_t_eva/pdtphys
2786    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldteva)
2787    CALL vect_mat(detr_therm, dx, dy, ddz, bdy, pwrf_grid%ld_th)
2788    CALL vect_mat(cldemi, dx, dy, ddz, bdy, pwrf_grid%lcldemi)
2789    CALL vect_mat(cldtau, dx, dy, ddz, bdy, pwrf_grid%lcldtau)
2790    CALL vect_mat(clwcon, dx, dy, ddz, bdy, pwrf_grid%lclwcon)
2791    var3d = d_t_lif/pdtphys
2792    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlif)
2793    var3d = d_t_lsc/pdtphys
2794    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlsc)
2795    var3d = (d_t_lsc+d_t_eva)/pdtphys
2796    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlschr)
2797    var3d = d_t_lscst/pdtphys
2798    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlscst)
2799    var3d = d_t_lscth/pdtphys
2800    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlscth)
2801    var3d=-1.*cool0/RDAY
2802    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlw0)
2803    var3d = -1.*cool/RDAY
2804    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlwr)
2805    var3d=d_t_oro/pdtphys
2806    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtoro)
2807    CALL vect_mat(d_t, dx, dy, ddz, bdy, pwrf_grid%ldtphy)
2808    var3d=heat0/RDAY
2809    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtsw0)
2810    var3d=heat/RDAY
2811    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtswr)
2812    var3d = d_t_ajs/pdtphys - d_t_ajsb/pdtphys
2813    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtthe)
2814    var3d=d_t_vdf/pdtphys
2815    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtvdf)
2816    var3d=d_t_wake/pdtphys
2817    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtwak)
2818    var3d=d_u_con/pdtphys
2819    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lducon)
2820    CALL vect_mat(d_u_dyn, dx, dy, ddz, bdy, pwrf_grid%ldudyn)
2821    var3d=d_u_lif/pdtphys
2822    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldulif)
2823    var3d=d_u_oro/pdtphys
2824    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lduoro)
2825    var3d=d_u_vdf/pdtphys
2826    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lduvdf)
2827    var3d=d_v_con/pdtphys
2828    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvcon)
2829    CALL vect_mat(d_v_dyn, dx, dy, ddz, bdy, pwrf_grid%ldvdyn)
2830    var3d=d_v_lif/pdtphys
2831    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvlif)
2832    var3d=d_v_oro/pdtphys
2833    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvoro)
2834    var3d=d_v_vdf/pdtphys
2835    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvvdf)
2836!NOTfound!    CALL vect_mat(ec550aer, dx, dy, ddz, bdy, pwrf_grid%lec550aer)
2837    CALL vect_mat(entr_therm, dx, dy, ddz, bdy, pwrf_grid%le_th)
2838    CALL vect_mat(coefm(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%levu)
2839    CALL vect_mat(fl, dx, dy, ddz, bdy, pwrf_grid%lfl)
2840    CALL vect_mat(zphi, dx, dy, ddz, bdy, pwrf_grid%lgeop)
2841    var3d = q_seri + ql_seri
2842    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lh2o)
2843    CALL vect_mat(fiwc, dx, dy, ddz, bdy, pwrf_grid%liwcon)
2844    CALL vect_mat(coefh(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%lkz)
2845    CALL vect_mat(coefh(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%lkz_max)
2846!    CALL vect_mat(lambda_th, dx, dy, ddz, bdy, pwrf_grid%llambda_th)
2847    CALL vect_mat(flwc, dx, dy, ddz, bdy, pwrf_grid%llwcon)
2848!NOTknown!    CALL vect_mat(ma, dx, dy, ddz, bdy, pwrf_grid%lma)
2849    CALL vect_mat(zmasse, dx, dy, ddz, bdy, pwrf_grid%lmass)
2850!NOTknown!    CALL vect_mat(mc, dx, dy, ddz, bdy, pwrf_grid%lmc)
2851!NOTknown!    CALL vect_mat(mcd, dx, dy, ddz, bdy, pwrf_grid%lmcd)
2852    CALL vect_mat(ql_seri, dx, dy, ddz, bdy, pwrf_grid%loliq)
2853    CALL vect_mat(q_seri, dx, dy, ddz, bdy, pwrf_grid%lovap)
2854! L. Fita, LMD January 2014. d_qx is the tendency of moisture due to physiq!
2855!   MOISTtend in the WRF_LMDZ coupling
2856    CALL vect_mat(qx(:,:,ivap), dx, dy, ddz, bdy, pwrf_grid%lovapinit)
2857!Notnow!    var3d = wo(:,:,1)*dobson_u*1e3 / zmasse / rmo3 * rmd)
2858!Notnow!    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lozone)
2859    CALL vect_mat_zstagg(paprs, dx, dy, ddz+1, bdy, pwrf_grid%lpaprs)
2860!NOTknown!    CALL vect_mat(pb, dx, dy, ddz, bdy, pwrf_grid%lpb)
2861    CALL vect_mat_zstagg(pmflxs, dx, dy, ddz+1, bdy, pwrf_grid%lpr_con_i)
2862    CALL vect_mat_zstagg(pmflxr, dx, dy, ddz+1, bdy, pwrf_grid%lpr_con_l)
2863    CALL vect_mat(pplay, dx, dy, ddz, bdy, pwrf_grid%lpres)
2864    CALL vect_mat_zstagg(psfl, dx, dy, ddz+1, bdy, pwrf_grid%lpr_lsc_i)
2865    CALL vect_mat_zstagg(prfl, dx, dy, ddz+1, bdy, pwrf_grid%lpr_lsc_l)
2866    var3d = 0.
2867    WHERE (ptconv) var3d = 1.
2868    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lptconv)
2869    CALL vect_mat(zqasc, dx, dy, ddz, bdy, pwrf_grid%lq_th)
2870    CALL vect_mat(ratqs, dx, dy, ddz, bdy, pwrf_grid%lratqs)
2871    CALL vect_mat(re, dx, dy, ddz, bdy, pwrf_grid%lre)
2872    CALL vect_mat(ref_ice, dx, dy, ddz, bdy, pwrf_grid%lref_ice)
2873    CALL vect_mat(ref_liq, dx, dy, ddz, bdy, pwrf_grid%lref_liq)
2874    CALL vect_mat(zx_rh, dx, dy, ddz, bdy, pwrf_grid%lrhum)
2875    CALL vect_mat(lwdn, dx, dy, ddz, bdy, pwrf_grid%lrld)
2876    CALL vect_mat(lwdn0, dx, dy, ddz, bdy, pwrf_grid%lrldcs)
2877    CALL vect_mat(lwup, dx, dy, ddz, bdy, pwrf_grid%lrlu)
2878    CALL vect_mat(lwup0, dx, dy, ddz, bdy, pwrf_grid%lrlucs)
2879!NOTknown!    CALL vect_mat(rn, dx, dy, ddz, bdy, pwrf_grid%lrn)
2880    CALL vect_mat(cldfra, dx, dy, ddz, bdy, pwrf_grid%lrneb)
2881    CALL vect_mat(rnebcon, dx, dy, ddz, bdy, pwrf_grid%lrnebcon)
2882    CALL vect_mat(rneb, dx, dy, ddz, bdy, pwrf_grid%lrnebls)
2883    CALL vect_mat(swdn, dx, dy, ddz, bdy, pwrf_grid%lrsd)
2884    CALL vect_mat(swdn0, dx, dy, ddz, bdy, pwrf_grid%lrsdcs)
2885    CALL vect_mat(swup, dx, dy, ddz, bdy, pwrf_grid%lrsu)
2886    CALL vect_mat(swup0, dx, dy, ddz, bdy, pwrf_grid%lrsucs)
2887    CALL vect_mat(fluxt(:,:,1), dx, dy, ddz, bdy, pwrf_grid%lsens_ter)
2888    CALL vect_mat(fluxt(:,:,2), dx, dy, ddz, bdy, pwrf_grid%lsens_lic)
2889    CALL vect_mat(fluxt(:,:,3), dx, dy, ddz, bdy, pwrf_grid%lsens_oce)
2890    CALL vect_mat(fluxt(:,:,4), dx, dy, ddz, bdy, pwrf_grid%lsens_sic)
2891    CALL vect_mat(t_seri, dx, dy, ddz, bdy, pwrf_grid%ltemp)
2892    CALL vect_mat(theta, dx, dy, ddz, bdy, pwrf_grid%ltheta)
2893    var3d=0.
2894    DO nsrf=1,nbsrf
2895       DO iz=1,ddz+1
2896         var3dstagg(:,iz)=var3dstagg(:,iz)+pctsrf(:,nsrf)*pbl_tke(:,iz,nsrf)
2897       ENDDO
2898    ENDDO
2899    CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%ltke)
2900    CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%ltke_max)
2901    CALL vect_mat(pbl_tke(:,:,1), dx, dy, ddz+1, bdy, pwrf_grid%ltke_ter)
2902    CALL vect_mat(pbl_tke(:,:,2), dx, dy, ddz+1, bdy, pwrf_grid%ltke_lic)
2903    CALL vect_mat(pbl_tke(:,:,3), dx, dy, ddz+1, bdy, pwrf_grid%ltke_oce)
2904    CALL vect_mat(pbl_tke(:,:,4), dx, dy, ddz+1, bdy, pwrf_grid%ltke_sic)
2905!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,1), dx, dy, ddz+1, bdy, pwrf_grid%ltke_ter_max)
2906!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,2), dx, dy, ddz+1, bdy, pwrf_grid%ltke_lic_max)
2907!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,3), dx, dy, ddz+1, bdy, pwrf_grid%ltke_oce_max)
2908!NOTinRegisty!    CALL vect_mat(pbl_tke(:,:,4), dx, dy, ddz+1, bdy, pwrf_grid%ltke_sic_max)
2909    var3d = d_qx(:,:,ivap)+d_q_dyn
2910    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhus)
2911    IF (iflag_thermals == 1) THEN
2912      var3d = d_q_con/pdtphys
2913    ELSE IF ((iflag_thermals > 1) .AND. (iflag_wake == 1)) THEN
2914      var3d = d_q_con/pdtphys + d_q_ajs/pdtphys + d_q_wake/pdtphys
2915    END IF
2916    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhusc)
2917    var3d = d_q_lsc/pdtphys + d_q_eva/pdtphys
2918    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhusscpbl)
2919    var3d = d_t + d_t_dyn
2920    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnt)
2921    IF (iflag_thermals == 1) THEN
2922      var3d = d_t_con/pdtphys + d_t_ajsb/pdtphys
2923    ELSE IF ((iflag_thermals  > 1) .AND. (iflag_wake == 1)) THEN
2924      var3d=d_t_con/pdtphys + d_t_ajs/pdtphys + d_t_wake/pdtphys
2925    END IF
2926    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntc)
2927    var3d = heat/RDAY - cool/RDAY
2928    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntr)
2929    var3d = (d_t_lsc + d_t_eva + d_t_vdf)/pdtphys
2930    CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntscpbl)
2931    CALL vect_mat(upwd, dx, dy, ddz, bdy, pwrf_grid%lupwd)
2932    CALL vect_mat(u_seri, dx, dy, ddz, bdy, pwrf_grid%lvitu)
2933    CALL vect_mat(v_seri, dx, dy, ddz, bdy, pwrf_grid%lvitv)
2934    CALL vect_mat(omega, dx, dy, ddz, bdy, pwrf_grid%lvitw)
2935    CALL vect_mat_zstagg(Vprecip, dx, dy, ddz+1, bdy, pwrf_grid%lvprecip)
2936    CALL vect_mat(wake_deltaq, dx, dy, ddz, bdy, pwrf_grid%lwake_deltaq)
2937    CALL vect_mat(wake_deltat, dx, dy, ddz, bdy, pwrf_grid%lwake_deltat)
2938    CALL vect_mat(wake_omg, dx, dy, ddz, bdy, pwrf_grid%lwake_omg)
2939    CALL vect_mat(wdtrainA, dx, dy, ddz, bdy, pwrf_grid%lwdtraina)
2940    CALL vect_mat(wdtrainM, dx, dy, ddz, bdy, pwrf_grid%lwdtrainm)
2941! L. Fita, LMD, January 2014. pphis is the given value of the pressure at the surface
2942    var3dstagg(:,1) = pphis(:)/RG
2943    DO iz=1,ddz
2944      var3dstagg(:,iz+1)= var3dstagg(:,iz)-(t_seri(:,iz)*RD*(paprs(:,iz+1)-          &
2945       paprs(:,iz)) )/ (pplay(:,iz)*RG)
2946    ENDDO
2947    CALL vect_mat(var3dstagg, dx, dy, ddz, bdy, pwrf_grid%lzfull)
2948    var3dstagg(:,1) = pphis(:)/RG-( (t_seri(:,1)+zxtsol)/2.*RD*(pplay(:,1)-          &
2949      paprs(:,1)) )/ ((paprs(:,1)+pplay(:,1))/2.*RG)
2950    DO iz=1, ddz-1
2951      var3dstagg(:,iz+1)= var3dstagg(:,iz)-( (t_seri(:,iz)+t_seri(:,iz+1))/          &
2952        2.*RD*(pplay(:,iz+1)-pplay(:,iz)))/( paprs(:,iz)*RG)
2953    ENDDO
2954    CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%lzhalf)
2955
2956    RETURN
2957
2958  END SUBROUTINE get_lmdz_out
2959
2960  SUBROUTINE put_lmdz_in_WRFout(is, ie, js, je, ks, ke, ddx, ddy, ddz, bdy, dlmdz,xp,&
2961    yp, lmdp, Nks, Nz0, nqtot, ivap, iliq, iflag_pbl, iflag_con, iflag_thermals,     &
2962    iflag_wake, itap, pdtphys, paprs, pplay, pphi, pphis, presnivs, qx, d_u, d_v,    &
2963    d_t, d_qx, d_ps, orchidlev, tsoilorchid, nslayers, wslay, wzs, wdzs, wnest_pos,  &
2964    wq2, wt2, wth2, wpsfc, wu10, wv10, wresm, wzetatop, wtslb, wsmois,               &
2965    wsh2o, wsmcrel, wsfroff, wudroff, wgrdflx, wacgrdflx, wsnow, wsnowh,             &
2966    wcanwat, wsst, wsstsk, wlai, wh_diabatic, wsina, wtsk, wtiso,                    &
2967    wmax_msftx, wmax_msfty, wrainc, wraincv, wrainsh, wrainshv, wrainnc, wrainncv,   &
2968    wsnownc, wsnowncv, wgraupelnc, wgraupelncv, whailnc, whailncv, wstepave_count,   &
2969    wcldfra, wswdown, wglw, wswnorm, wolr, wemiss, wtmn, wust, wpblh, whfx,          &
2970    wqfx, wlh, wachfx, waclhf, wsnowc, wsr, wpotevp, wsnopcx, wsoiltb)
2971!    wsave_topo_from_real, wavg_fuel_frac, wuah, wvah, wseed1, wseed2)
2972! Subroutine to get LMDZ output and get into WRF standard output variables
2973!   (as they appear in the wrfout)
2974
2975! LMDZ modules
2976    USE indice_sol_mod
2977    USE phys_state_var_mod
2978    USE phys_output_var_mod
2979!    USE pbl_surface_mod
2980    USE phys_local_var_mod
2981    USE comgeomphy
2982
2983! WRF_LMDZ modules
2984    USE output_lmdz_NOmodule
2985    USE lmdz_wrf_variables_mod
2986
2987    IMPLICIT NONE
2988
2989    INTEGER, INTENT(IN)                                  :: is, ie, js, je, ks, ke,  &
2990      ddx, ddy, ddz, bdy, dlmdz, xp, yp, lmdp, Nks, Nz0, nqtot, ivap, iliq,          &
2991      iflag_pbl, iflag_con, iflag_thermals, iflag_wake, itap, nslayers
2992    REAL, INTENT(IN)                                     :: pdtphys
2993    REAL, DIMENSION(dlmdz), INTENT(IN)                   :: pphis, d_ps
2994    REAL, DIMENSION(ddz), INTENT(IN)                     :: presnivs
2995    REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN)             :: paprs
2996    REAL, DIMENSION(dlmdz,ddz), INTENT(IN)               :: pplay, pphi, d_u, d_v, d_t
2997    REAL, DIMENSION(dlmdz,ddz,nqtot), INTENT(IN)         :: qx, d_qx
2998    REAL, DIMENSION(Nz0), INTENT(IN)                     :: orchidlev
2999    REAL, DIMENSION(dlmdz,Nz0,Nks), INTENT(IN)           :: tsoilorchid
3000    INTEGER, INTENT(INOUT)                               :: wstepave_count
3001    REAL, INTENT(INOUT)                                  :: wresm, wzetatop, wtiso,  &
3002      wmax_msftx, wmax_msfty
3003    REAL, DIMENSION(1:nslayers), INTENT(INOUT)           :: wzs, wdzs
3004    REAL, DIMENSION(1:nslayers), INTENT(IN)              :: wslay
3005    REAL, DIMENSION(is:ie,js:je), INTENT(INOUT)          :: wnest_pos, wq2, wt2,     &
3006      wth2, wpsfc, wu10, wv10, wsfroff, wudroff, wgrdflx, wacgrdflx, wsnow, wsnowh,  &
3007      wcanwat, wsst, wsstsk, wlai, wsina, wtsk, wrainc, wraincv, wrainsh, wrainshv,  &
3008      wrainnc, wrainncv, wsnownc, wsnowncv, wgraupelnc, wgraupelncv, whailnc,        &
3009      whailncv, wswdown, wglw, wswnorm, wolr, wemiss, wtmn, wust, wpblh, whfx,       &
3010      wqfx, wlh, wachfx, waclhf, wsnowc, wsr, wpotevp, wsnopcx, wsoiltb
3011    REAL, DIMENSION(is:ie,1:nslayers,js:je), INTENT(INOUT) :: wtslb, wsmois, wsh2o,  &
3012      wsmcrel
3013    REAL, DIMENSION(is:ie,ks:ke,js:je), INTENT(INOUT)    :: wh_diabatic, wcldfra
3014!    EXTERNAL suphel
3015
3016! Local
3017    INTEGER                                              :: dx, dy, ix, iy, ixy, iz
3018    INTEGER                                              :: nsrf
3019    REAL                                                 :: varixy
3020    REAL, DIMENSION(Nz0)                                 :: yvals1
3021    REAL                                                 :: RDAY, RG, RD
3022    CHARACTER(LEN=50)                                    :: fname
3023
3024!!!!!!! Variables
3025! pwrf_grid: WRF grid structure
3026
3027    fname = "'put_lmdz_in_WRFout'"
3028
3029! L. Fita, LMD January 2014. It should work using suphel, but not?
3030!    CALL suphel
3031
3032    RDAY = 86400.
3033    RG = 9.81
3034    RD = 287.
3035
3036    wdzs(1)=wslay(nslayers)
3037    DO iz=1,nslayers-1
3038      wdzs(iz+1) = wslay(nslayers-iz+1) - wslay(nslayers-iz)
3039    END DO
3040!    wmax_msftx =
3041!    wmax_msfty =
3042!    wresm =
3043!    wstepave_count =
3044!    wtiso =
3045!    wzetatop =
3046    wzs = wslay
3047
3048    DO ix=1,ddx
3049      DO iy=1,ddy
3050        ixy=ddx*(iy-1)+ix
3051        wacgrdflx(ix,iy) = wacgrdflx(ix,iy) + 0.
3052        wachfx(ix,iy) = wachfx(ix,iy) + sens(ixy)
3053        waclhf(ix,iy) = waclhf(ix,iy) + zxfluxlat(ixy)
3054        DO iz=1,ddz
3055          wcldfra(ix,iz,iy) = cldfra(ixy,iz)
3056        END DO
3057        wcanwat(ix,iy) = 0.
3058        wglw(ix,iy) = sollwdown(ixy)
3059! These are precipitation fluxes !!
3060!        wgraupelnc(ix,iy) = wgraupelnc(ix,iy) + prmflxs(ixy) + prfs(ixy,iz)
3061!        wgraupelncv(ix,iy) = prmflxs(ixy) + prfs(ixy,iz)
3062        wgrdflx(ix,iy) = 0.
3063        whailnc(ix,iy) = whailnc(ix,iy) + 0.
3064        whailncv(ix,iy) = 0.
3065        wh_diabatic(ix,ks:ke,iy) = 0.
3066        whfx(ix,iy) = sens(ixy)
3067        wlai(ix,iy) = 0.
3068        wlh(ix,iy) = zxfluxlat(ixy)
3069        wnest_pos(ix,iy) = 0.
3070!!        wnoahres(ix,iy) = 0.
3071        wolr(ix,iy) = toplw(ixy)
3072        wpblh(ix,iy) = s_pblh(ixy)
3073! Should be something from evappot_srf ?
3074        wpotevp(ix,iy) = wpotevp(ix,iy) + 0.
3075        wpsfc(ix,iy) = paprs(ixy,1)
3076        wq2(ix,iy) = zq2m(ixy)
3077        wqfx(ix,iy) = 0.
3078        wrainc(ix,iy) = wrainc(ix,iy) + (rain_con(ixy) + snow_con(ixy))*pdtphys
3079        wraincv(ix,iy) = (rain_con(ixy) + snow_con(ixy))*pdtphys
3080        wrainnc(ix,iy) = wrainnc(ix,iy) + (rain_lsc(ixy)+snow_lsc(ixy))*pdtphys
3081        wrainncv(ix,iy) = (rain_lsc(ixy)+snow_lsc(ixy))*pdtphys
3082        wrainsh(ix,iy) = wrainsh(ix,iy) + 0.
3083        wrainshv(ix,iy) = 0.
3084! L. Fita, LMD. February 2014
3085! LMDZ run_off_lic is for continental ice. We could try to get run_off_ter, but....
3086!        wsfroff(ix,iy) = 0.
3087        DO iz=1,Nks
3088          yvals1(iz) = tsoilorchid(ixy,iz,1)*pctsrf(ixy,is_ter) +                    &
3089            tsoilorchid(ixy,iz,2)*pctsrf(ixy,is_lic) +                               &
3090            tsoilorchid(ixy,iz,3)*pctsrf(ixy,is_oce) +                               &
3091            tsoilorchid(ixy,iz,4)*pctsrf(ixy,is_sic)
3092        END DO
3093        CALL interpolate_layers(Nz0,orchidlev,yvals1,nslayers,wslay,wtslb(ix,:,iy))
3094! No way ???
3095!          wsh2o(is:ie,iz,js:je) = wsh2o(is:ie,iz,js:je) +
3096!          wsmcrel(is:ie,iz,js:je) = 0.
3097
3098!        wsinalpha(ix,iy) =
3099        wsina(ix,iy) = 0.
3100        wsnopcx(ix,iy) = wsnopcx(ix,iy) + 0.
3101!!        wsnowh(ix,iy) = snowght(ixy)
3102        wsnownc(ix,iy) = wsnownc(ix,iy) + snow_fall(ixy)*pdtphys
3103        wsnowncv(ix,iy) = snow_fall(ixy)*pdtphys
3104        wsoiltb(ix,iy) = wtslb(ix,nslayers,iy)
3105        IF (rain_fall(ixy) /= 0.) THEN
3106          wsr(ix,iy) = ( snow_fall(ixy) ) / ( rain_fall(ixy))
3107        ELSE
3108          wsr(ix,iy) = 0.
3109        END IF
3110        wsstsk(ix,iy) = ftsol(ixy,is_oce)
3111        wsst(ix,iy) = sst(ixy)
3112        wswdown(ix,iy) = swdn(ixy,1)*pctsrf(ixy,is_ter) +                            &
3113          swdn(ixy,2)*pctsrf(ixy,is_lic) + swdn(ixy,3)*pctsrf(ixy,is_oce) +          &
3114          swdn(ixy,4)*pctsrf(ixy,is_sic)
3115        wswnorm(ix,iy) = 0.
3116        wt2(ix,iy) = zt2m(ixy)
3117! This is not exact!!
3118        CALL temp_to_theta1(zt2m(ixy), paprs(ixy,1), wth2(ix,iy))
3119        wtsk(ix,iy) = ftsol(ixy, is_ter)*pctsrf(ixy,is_ter) + ftsol(ixy, is_lic)*    &
3120          pctsrf(ixy,is_lic) + ftsol(ixy, is_oce)*pctsrf(ixy,is_oce) +               &
3121          ftsol(ixy, is_sic)*pctsrf(ixy,is_sic)
3122        wu10(ix,iy) = zu10m(ixy)
3123        wudroff(ix,iy) = wudroff(ix,iy) + 0.
3124        wv10(ix,iy) = zv10m(ixy)
3125      END DO
3126    END DO
3127
3128  END SUBROUTINE put_lmdz_in_WRFout
3129
3130END MODULE wrf_lmdz_mod
3131
Note: See TracBrowser for help on using the repository browser.