source: LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90 @ 5159

Last change on this file since 5159 was 5159, checked in by abarral, 7 weeks ago

Put dimensions.h and paramet.h into modules
File size: 23.3 KB
Line 
1! $Id$
2
3MODULE etat0phys
4
5  !*******************************************************************************
6  ! Purpose: Create physical initial state using atmospheric fields from a
7  !          database of atmospheric to initialize the model.
8  !-------------------------------------------------------------------------------
9  ! Comments:
10
11  !    *  This module is designed to work for Earth (and with ioipsl)
12
13  !    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
14  !         "start_init_phys" for variables contained in file "ECPHY.nc":
15  !            'ST'     : Surface temperature
16  !            'CDSW'   : Soil moisture
17  !         "start_init_orog" for variables contained in file "Relief.nc":
18  !            'RELIEF' : High resolution orography
19
20  !    * The land mask and corresponding weights can be:
21  !      1) computed using the ocean mask from the ocean model (to ensure ocean
22  !         fractions are the same for atmosphere and ocean) for coupled runs.
23  !         File name: "o2a.nc"  ;  variable name: "OceMask"
24  !      2) computed from topography file "Relief.nc" for forced runs.
25
26  !    * Allowed values for read_climoz flag are 0, 1 and 2:
27  !      0: do not read an ozone climatology
28  !      1: read a single ozone climatology that will be used day and night
29  !      2: read two ozone climatologies, the average day and night climatology
30  !         and the daylight climatology
31  !-------------------------------------------------------------------------------
32  !    * There is a big mess with the longitude size. Should it be iml or iml+1 ?
33  !  I have chosen to use the iml+1 as an argument to this routine and we declare
34  !  internaly smaller fields when needed. This needs to be cleared once and for
35  !  all in LMDZ. A convention is required.
36  !-------------------------------------------------------------------------------
37
38  USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo
39  USE lmdz_assert_eq, ONLY: assert_eq
40  USE dimphy, ONLY: klon
41  USE conf_dat_m, ONLY: conf_dat2d
42  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
43          solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, &
44          sollw, sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs, w01, &
45          sig1, ftsol, clwcon, fm_therm, wake_Cstar, pctsrf, entr_therm, radpas, f0, &
46          zmax0, fevap, rnebcon, falb_dir, falb_dif, wake_fip, agesno, detr_therm, pbl_tke, &
47          phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
48          prw_ancien, u10m, v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
49          ale_bl, ale_bl_trig, alp_bl, &
50          ale_wake, ale_bl_stat, AWAKE_S
51
52  USE comconst_mod, ONLY: pi, dtvr
53  USE lmdz_iniprint, ONLY: lunout, prt_level
54  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
55  USE lmdz_comgeom2
56  USE lmdz_clesphys
57  USE lmdz_paramet
58
59  PRIVATE
60  PUBLIC :: etat0phys_netcdf
61
62
63  INCLUDE "dimsoil.h"
64  REAL, SAVE :: deg2rad
65  REAL, SAVE, ALLOCATABLE :: tsol(:)
66  INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
67  REAL, ALLOCATABLE, SAVE :: lon_phys(:, :), lat_phys(:, :), levphys_ini(:)
68  CHARACTER(LEN = 256), PARAMETER :: oroparam = "oro_params.nc"
69  CHARACTER(LEN = 256), PARAMETER :: orofname = "Relief.nc", orogvar = "RELIEF"
70  CHARACTER(LEN = 256), PARAMETER :: phyfname = "ECPHY.nc", psrfvar = "SP"
71  CHARACTER(LEN = 256), PARAMETER :: qsolvar = "CDSW", tsrfvar = "ST"
72
73
74CONTAINS
75
76
77  !-------------------------------------------------------------------------------
78
79  SUBROUTINE etat0phys_netcdf(masque, phis)
80
81    !-------------------------------------------------------------------------------
82    ! Purpose: Creates initial states
83    !-------------------------------------------------------------------------------
84    ! Notes:  1) This routine is designed to work for Earth
85    !         2) If masque(:,:)/=-99999., masque and phis are already known.
86    !         Otherwise: compute it.
87    !-------------------------------------------------------------------------------
88    USE control_mod
89    USE fonte_neige_mod
90    USE pbl_surface_mod
91    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
92    USE indice_sol_mod
93    USE conf_phys_m, ONLY: conf_phys
94    USE init_ssrf_m, ONLY: start_init_subsurf
95    USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
96            ratqs_inter_, rneb_ancien
97    USE lmdz_alpale
98    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
99
100    IMPLICIT NONE
101    !-------------------------------------------------------------------------------
102    ! Arguments:
103    REAL, INTENT(INOUT) :: masque(:, :) !--- Land mask           dim(iip1,jjp1)
104    REAL, INTENT(INOUT) :: phis  (:, :) !--- Ground geopotential dim(iip1,jjp1)
105    !-------------------------------------------------------------------------------
106    ! Local variables:
107    CHARACTER(LEN = 256) :: modname = "etat0phys_netcdf", fmt
108    INTEGER :: i, j, l, ji, iml, jml
109    LOGICAL :: read_mask
110    REAL :: phystep, dummy
111    REAL, DIMENSION(SIZE(masque, 1), SIZE(masque, 2)) :: masque_tmp, phiso
112    REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder
113    REAL, DIMENSION(klon, nbsrf) :: qsurf, snsrf
114    REAL, DIMENSION(klon, nsoilmx, nbsrf) :: tsoil
115
116    !--- Arguments for conf_phys
117    LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
118    REAL :: solarlong0, seuil_inversion, fact_cldcon, facttemps
119    LOGICAL :: ok_newmicro
120    INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
121    REAL :: ratqsbas, ratqshaut, tau_ratqs
122    LOGICAL :: ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
123    INTEGER :: flag_aerosol
124    INTEGER :: flag_aerosol_strat
125    INTEGER :: flag_volc_surfstrat
126    LOGICAL :: flag_aer_feedback
127    LOGICAL :: flag_bc_internal_mixture
128    REAL :: bl95_b0, bl95_b1
129    INTEGER :: read_climoz                        !--- Read ozone climatology
130    REAL :: alp_offset
131    LOGICAL :: filtre_oro = .FALSE.
132
133    deg2rad = pi / 180.0
134    iml = assert_eq(SIZE(masque, 1), SIZE(phis, 1), TRIM(modname) // " iml")
135    jml = assert_eq(SIZE(masque, 2), SIZE(phis, 2), TRIM(modname) // " jml")
136
137    ! Physics configuration
138    !*******************************************************************************
139    CALL conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
140            callstats, &
141            solarlong0, seuil_inversion, &
142            fact_cldcon, facttemps, ok_newmicro, iflag_radia, &
143            iflag_cldcon, &
144            ratqsbas, ratqshaut, tau_ratqs, &
145            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, &
146            aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat, &
147            flag_aer_feedback, flag_bc_internal_mixture, bl95_b0, bl95_b1, &
148            read_climoz, alp_offset)
149    CALL phys_state_var_init(read_climoz)
150
151    !--- Initial atmospheric CO2 conc. from .def file
152    co2_ppm0 = co2_ppm
153
154    ! Compute ground geopotential, sub-cells quantities and possibly the mask.
155    !*******************************************************************************
156    read_mask = ANY(masque/=-99999.); masque_tmp = masque
157    CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
158
159    CALL getin('filtre_oro', filtre_oro)
160    IF (filtre_oro) CALL filtreoro(size(phis, 1), size(phis, 2), phis, masque_tmp, rlatu)
161
162    WRITE(fmt, "(i4,'i1)')")iml ; fmt = '(' // ADJUSTL(fmt)
163    IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
164      masque = masque_tmp
165      IF(prt_level>=1) THEN
166        WRITE(lunout, *)'BUILT MASK :'
167        WRITE(lunout, fmt) NINT(masque)
168      END IF
169      WHERE(masque(:, :)<EPSFRA) masque(:, :) = 0.
170      WHERE(1. - masque(:, :)<EPSFRA) masque(:, :) = 1.
171    END IF
172    CALL gr_dyn_fi(1, iml, jml, klon, masque, zmasq) !--- Land mask to physical grid
173
174    ! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
175    !*******************************************************************************
176    CALL start_init_phys(rlonu, rlatv, phis)
177
178    ! Some initializations.
179    !*******************************************************************************
180    sn    (:) = 0.0                               !--- Snow
181    radsol(:) = 0.0                               !--- Net radiation at ground
182    rugmer(:) = 0.001                             !--- Ocean rugosity
183    !--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
184    IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz, ok_daily_climoz)
185
186    ! Sub-surfaces initialization.
187    !*******************************************************************************
188    CALL start_init_subsurf(read_mask)
189
190    ! Write physical initial state
191    !*******************************************************************************
192    WRITE(lunout, *)'phystep ', dtvr, iphysiq, nbapp_rad
193    phystep = dtvr * FLOAT(iphysiq)
194    radpas = NINT (86400. / phystep / FLOAT(nbapp_rad))
195    WRITE(lunout, *)'phystep =', phystep, radpas
196
197    ! Init: ftsol, snsrf, qsurf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
198    !*******************************************************************************
199    DO i = 1, nbsrf; ftsol(:, i) = tsol;
200    END DO
201    DO i = 1, nbsrf; snsrf(:, i) = sn;
202    END DO
203    falb_dir(:, :, is_ter) = 0.08
204    falb_dir(:, :, is_lic) = 0.6
205    falb_dir(:, :, is_oce) = 0.5
206    falb_dir(:, :, is_sic) = 0.6
207
208    !ym warning missing init for falb_dif => set to 0
209    falb_dif(:, :, :) = 0
210
211    u10m(:, :) = 0
212    v10m(:, :) = 0
213    treedrg(:, :, :) = 0
214
215    fevap(:, :) = 0.
216    qsurf = 0.
217    DO i = 1, nbsrf; DO j = 1, nsoilmx; tsoil(:, j, i) = tsol;
218    END DO;
219    END DO
220    rain_fall = 0.
221    snow_fall = 0.
222    solsw = 165.
223    solswfdiff = 1.
224    sollw = -53.
225    !ym warning missing init for sollwdown => set to 0
226    sollwdown = 0.
227    t_ancien = 273.15
228    q_ancien = 0.
229    ql_ancien = 0.
230    qs_ancien = 0.
231    prlw_ancien = 0.
232    prsw_ancien = 0.
233    prw_ancien = 0.
234    agesno = 0.
235
236    u_ancien = 0.
237    v_ancien = 0.
238    wake_delta_pbl_TKE(:, :, :) = 0
239    wake_dens(:) = 0
240    awake_dens = 0.
241    cv_gen = 0.
242    ale_bl = 0.
243    ale_bl_trig = 0.
244    alp_bl = 0.
245    ale_wake = 0.
246    ale_bl_stat = 0.
247
248    z0m(:, :) = 0 ! ym missing 5th subsurface initialization
249
250    z0m(:, is_oce) = rugmer(:)
251    z0m(:, is_ter) = 0.01 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
252    z0m(:, is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
253    z0m(:, is_sic) = 0.001
254    z0h(:, :) = z0m(:, :)
255
256    fder = 0.0
257    clwcon = 0.0
258    rnebcon = 0.0
259    ratqs = 0.0
260    run_off_lic_0 = 0.0
261    rugoro = 0.0
262
263    ! Before phyredem calling, surface modules and values to be saved in startphy.nc
264    ! are initialized
265    !*******************************************************************************
266    dummy = 1.0
267    pbl_tke(:, :, :) = 1.e-8
268    zmax0(:) = 40.
269    f0(:) = 1.e-5
270    sig1(:, :) = 0.
271    w01(:, :) = 0.
272    wake_deltat(:, :) = 0.
273    wake_deltaq(:, :) = 0.
274    wake_s(:) = 0.
275    wake_cstar(:) = 0.
276    wake_fip(:) = 0.
277    wake_pe = 0.
278    fm_therm = 0.
279    entr_therm = 0.
280    detr_therm = 0.
281    awake_s = 0.
282
283    CALL fonte_neige_init(run_off_lic_0)
284    CALL pbl_surface_init(fder, snsrf, qsurf, tsoil)
285
286    IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) THEN
287      delta_tsurf = 0.
288      beta_aridity = 0.
289    end IF
290
291    ratqs_inter_ = 0.002
292    rneb_ancien = 0.
293    CALL phyredem("startphy.nc")
294
295    !  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
296    !  WRITE(lunout,*)'entree histclo'
297    CALL histclo()
298
299  END SUBROUTINE etat0phys_netcdf
300
301  !-------------------------------------------------------------------------------
302
303
304  !-------------------------------------------------------------------------------
305
306  SUBROUTINE start_init_orog(lon_in, lat_in, phis, masque)
307
308    !===============================================================================
309    ! Comment:
310    !   This routine launch grid_noro, which computes parameters for SSO scheme as
311    !   described in LOTT & MILLER (1997) and LOTT(1999).
312    !   In case the file oroparam is present and the key read_orop is activated,
313    !   grid_noro is bypassed and sub-cell parameters are read from the file.
314    !===============================================================================
315    USE grid_noro_m, ONLY: grid_noro, read_noro
316    USE logic_mod, ONLY: read_orop
317    IMPLICIT NONE
318    !-------------------------------------------------------------------------------
319    ! Arguments:
320    REAL, INTENT(IN) :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
321    REAL, INTENT(INOUT) :: phis(:, :), masque(:, :) ! dim (iml,jml)
322    !-------------------------------------------------------------------------------
323    ! Local variables:
324    CHARACTER(LEN = 256) :: modname
325    INTEGER :: fid, llm_tmp, ttm_tmp, iml, jml, iml_rel, jml_rel, itau(1)
326    INTEGER :: ierr
327    REAL :: lev(1), date, dt
328    REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:, :), relief_hi(:, :)
329    REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:, :), tmp_var  (:, :)
330    REAL, ALLOCATABLE :: zmea0(:, :), zstd0(:, :), zsig0(:, :)
331    REAL, ALLOCATABLE :: zgam0(:, :), zthe0(:, :), zpic0(:, :), zval0(:, :)
332    !-------------------------------------------------------------------------------
333    modname = "start_init_orog"
334    iml = assert_eq(SIZE(lon_in), SIZE(phis, 1), SIZE(masque, 1), TRIM(modname) // " iml")
335    jml = assert_eq(SIZE(lat_in), SIZE(phis, 2), SIZE(masque, 2), TRIM(modname) // " jml")
336
337    !--- HIGH RESOLUTION OROGRAPHY
338    CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
339
340    ALLOCATE(lat_rel(iml_rel, jml_rel), lon_rel(iml_rel, jml_rel))
341    CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel, &
342            lev, ttm_tmp, itau, date, dt, fid)
343    ALLOCATE(relief_hi(iml_rel, jml_rel))
344    CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
345    CALL flinclo(fid)
346
347    !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
348    ALLOCATE(lon_ini(iml_rel), lat_ini(jml_rel))
349    lon_ini(:) = lon_rel(:, 1); IF(MAXVAL(lon_rel)>pi) lon_ini = lon_ini * deg2rad
350    lat_ini(:) = lat_rel(1, :); IF(MAXVAL(lat_rel)>pi) lat_ini = lat_ini * deg2rad
351
352    !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
353    ALLOCATE(lon_rad(iml_rel), lat_rad(jml_rel))
354    CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
355    DEALLOCATE(lon_ini, lat_ini)
356
357    !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
358    WRITE(lunout, *)
359    WRITE(lunout, *)'*** Compute parameters needed for gravity wave drag code ***'
360
361    !--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
362    ALLOCATE(zmea0(iml, jml), zstd0(iml, jml)) !--- Mean orography and std deviation
363    ALLOCATE(zsig0(iml, jml), zgam0(iml, jml)) !--- Slope and nisotropy
364    zsig0(:, :) = 0   !ym uninitialized variable
365    zgam0(:, :) = 0   !ym uninitialized variable
366    ALLOCATE(zthe0(iml, jml))                !--- Highest slope orientation
367    zthe0(:, :) = 0   !ym uninitialized variable
368    ALLOCATE(zpic0(iml, jml), zval0(iml, jml)) !--- Peaks and valley heights
369
370    !--- READ SUB-CELL SCALES PARAMETERS FROM A FILE (AT RIGHT RESOLUTION)
371    OPEN(UNIT = 66, FILE = oroparam, STATUS = 'OLD', IOSTAT = ierr)
372    IF(ierr==0.AND.read_orop) THEN
373      CLOSE(UNIT = 66)
374      CALL read_noro(lon_in, lat_in, oroparam, &
375              phis, zmea0, zstd0, zsig0, zgam0, zthe0, zpic0, zval0, masque)
376    ELSE
377      !--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
378      CALL grid_noro(lon_rad, lat_rad, relief_hi, lon_in, lat_in, &
379              phis, zmea0, zstd0, zsig0, zgam0, zthe0, zpic0, zval0, masque)
380    END IF
381    phis = phis * 9.81
382    phis(iml, :) = phis(1, :)
383    DEALLOCATE(relief_hi, lon_rad, lat_rad)
384
385    !--- PUT QUANTITIES TO PHYSICAL GRID
386    CALL gr_dyn_fi(1, iml, jml, klon, zmea0, zmea); DEALLOCATE(zmea0)
387    CALL gr_dyn_fi(1, iml, jml, klon, zstd0, zstd); DEALLOCATE(zstd0)
388    CALL gr_dyn_fi(1, iml, jml, klon, zsig0, zsig); DEALLOCATE(zsig0)
389    CALL gr_dyn_fi(1, iml, jml, klon, zgam0, zgam); DEALLOCATE(zgam0)
390    CALL gr_dyn_fi(1, iml, jml, klon, zthe0, zthe); DEALLOCATE(zthe0)
391    CALL gr_dyn_fi(1, iml, jml, klon, zpic0, zpic); DEALLOCATE(zpic0)
392    CALL gr_dyn_fi(1, iml, jml, klon, zval0, zval); DEALLOCATE(zval0)
393
394  END SUBROUTINE start_init_orog
395
396  !-------------------------------------------------------------------------------
397
398
399  !-------------------------------------------------------------------------------
400
401  SUBROUTINE start_init_phys(lon_in, lat_in, phis)
402
403    !===============================================================================
404    ! Purpose:   Compute tsol and qsol, knowing phis.
405    !===============================================================================
406    IMPLICIT NONE
407    !-------------------------------------------------------------------------------
408    ! Arguments:
409    REAL, INTENT(IN) :: lon_in(:), lat_in(:)       ! dim (iml) (jml2)
410    REAL, INTENT(IN) :: phis(:, :)                   ! dim (iml,jml)
411    !-------------------------------------------------------------------------------
412    ! Local variables:
413    CHARACTER(LEN = 256) :: modname
414    REAL :: date, dt
415    INTEGER :: iml, jml, jml2, itau(1)
416    REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:, :)
417    REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
418    REAL, ALLOCATABLE :: ts(:, :), qs(:, :)
419    !-------------------------------------------------------------------------------
420    modname = "start_init_phys"
421    iml = assert_eq(SIZE(lon_in), SIZE(phis, 1), TRIM(modname) // " iml")
422    jml = SIZE(phis, 2); jml2 = SIZE(lat_in)
423
424    WRITE(lunout, *)'Opening the surface analysis'
425    CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
426    WRITE(lunout, *) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys
427
428    ALLOCATE(lat_phys(iml_phys, jml_phys), lon_phys(iml_phys, jml_phys))
429    ALLOCATE(levphys_ini(llm_phys))
430    CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys, &
431            lon_phys, lat_phys, levphys_ini, ttm_phys, itau, date, dt, fid_phys)
432
433    !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
434    ALLOCATE(lon_ini(iml_phys), lat_ini(jml_phys))
435    lon_ini(:) = lon_phys(:, 1); IF(MAXVAL(lon_phys)>pi) lon_ini = lon_ini * deg2rad
436    lat_ini(:) = lat_phys(1, :); IF(MAXVAL(lat_phys)>pi) lat_ini = lat_ini * deg2rad
437
438    ALLOCATE(var_ana(iml_phys, jml_phys), lon_rad(iml_phys), lat_rad(jml_phys))
439    CALL get_var_phys(tsrfvar, ts)                   !--- SURFACE TEMPERATURE
440    CALL get_var_phys(qsolvar, qs)                   !--- SOIL MOISTURE
441    CALL flinclo(fid_phys)
442    DEALLOCATE(var_ana, lon_rad, lat_rad, lon_ini, lat_ini)
443
444    !--- TSOL AND QSOL ON PHYSICAL GRID
445    ALLOCATE(tsol(klon))
446    CALL gr_dyn_fi(1, iml, jml, klon, ts, tsol)
447    CALL gr_dyn_fi(1, iml, jml, klon, qs, qsol)
448    DEALLOCATE(ts, qs)
449
450  CONTAINS
451
452    !-------------------------------------------------------------------------------
453
454    SUBROUTINE get_var_phys(title, field)
455
456      !-------------------------------------------------------------------------------
457      IMPLICIT NONE
458      !-------------------------------------------------------------------------------
459      ! Arguments:
460      CHARACTER(LEN = *), INTENT(IN) :: title
461      REAL, ALLOCATABLE, INTENT(INOUT) :: field(:, :)
462      !-------------------------------------------------------------------------------
463      ! Local variables:
464      INTEGER :: tllm
465      !-------------------------------------------------------------------------------
466      SELECT CASE(title)
467      CASE(psrfvar);         tllm = 0
468      CASE(tsrfvar, qsolvar); tllm = llm_phys
469      END SELECT
470      IF(ALLOCATED(field)) RETURN
471      ALLOCATE(field(iml, jml)); field(:, :) = 0.
472      CALL flinget(fid_phys, title, iml_phys, jml_phys, tllm, ttm_phys, 1, 1, var_ana)
473      CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
474      CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, &
475              lon_in, lat_in, field)
476
477    END SUBROUTINE get_var_phys
478
479    !-------------------------------------------------------------------------------
480
481  END SUBROUTINE start_init_phys
482
483  !-------------------------------------------------------------------------------
484
485
486  !-------------------------------------------------------------------------------
487
488  SUBROUTINE interp_startvar(nam, ibeg, lon, lat, vari, lon2, lat2, varo)
489
490    !-------------------------------------------------------------------------------
491    USE inter_barxy_m, ONLY: inter_barxy
492    IMPLICIT NONE
493    !-------------------------------------------------------------------------------
494    ! Arguments:
495    CHARACTER(LEN = *), INTENT(IN) :: nam
496    LOGICAL, INTENT(IN) :: ibeg
497    REAL, INTENT(IN) :: lon(:), lat(:)   ! dim (ii) (jj)
498    REAL, INTENT(IN) :: vari(:, :)        ! dim (ii,jj)
499    REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
500    REAL, INTENT(OUT) :: varo(:, :)        ! dim (i1) (j1)
501    !-------------------------------------------------------------------------------
502    ! Local variables:
503    CHARACTER(LEN = 256) :: modname
504    INTEGER :: ii, jj, i1, j1, j2
505    REAL, ALLOCATABLE :: vtmp(:, :)
506    !-------------------------------------------------------------------------------
507    modname = "interp_startvar"
508    ii = assert_eq(SIZE(lon), SIZE(vari, 1), TRIM(modname) // " ii")
509    jj = assert_eq(SIZE(lat), SIZE(vari, 2), TRIM(modname) // " jj")
510    i1 = assert_eq(SIZE(lon2), SIZE(varo, 1), TRIM(modname) // " i1")
511    j1 = SIZE(varo, 2); j2 = SIZE(lat2)
512    ALLOCATE(vtmp(i1 - 1, j1))
513    IF(ibeg.AND.prt_level>1) THEN
514      WRITE(lunout, *)"--------------------------------------------------------"
515      WRITE(lunout, *)"$$$ Interpolation barycentrique pour " // TRIM(nam) // " $$$"
516      WRITE(lunout, *)"--------------------------------------------------------"
517    END IF
518    CALL inter_barxy(lon, lat(:jj - 1), vari, lon2(:i1 - 1), lat2, vtmp)
519    CALL gr_int_dyn(vtmp, varo, i1 - 1, j1)
520
521  END SUBROUTINE interp_startvar
522
523  !-------------------------------------------------------------------------------
524
525  !*******************************************************************************
526
527  SUBROUTINE filtreoro(imp1, jmp1, phis, masque, rlatu)
528
529    IMPLICIT NONE
530
531    INTEGER imp1, jmp1
532    REAL, DIMENSION(imp1, jmp1) :: phis, masque
533    REAL, DIMENSION(jmp1) :: rlatu
534    REAL, DIMENSION(imp1) :: wwf
535    REAL, DIMENSION(imp1, jmp1) :: phiso
536    INTEGER :: ifiltre, ifi, ii, i, j
537    REAL :: coslat0, ssz
538
539    coslat0 = 0.5
540    phiso = phis
541    DO j = 2, jmp1 - 1
542      PRINT*, 'avant if ', cos(rlatu(j)), coslat0
543      IF (cos(rlatu(j))<coslat0) THEN
544        ! nb de pts affectes par le filtrage de part et d'autre du pt
545        ifiltre = (coslat0 / cos(rlatu(j)) - 1.) / 2.
546        wwf = 0.
547        DO i = 1, ifiltre
548          wwf(i) = 1.
549        enddo
550        wwf(ifiltre + 1) = (coslat0 / cos(rlatu(j)) - 1.) / 2. - ifiltre
551        DO i = 1, imp1 - 1
552          IF (masque(i, j)>0.9) THEN
553            ssz = phis(i, j)
554            DO ifi = 1, ifiltre + 1
555              ii = i + ifi
556              IF (ii>imp1 - 1) ii = ii - imp1 + 1
557              ssz = ssz + wwf(ifi) * phis(ii, j)
558              ii = i - ifi
559              IF (ii<1) ii = ii + imp1 - 1
560              ssz = ssz + wwf(ifi) * phis(ii, j)
561            enddo
562            phis(i, j) = ssz * cos(rlatu(j)) / coslat0
563          endif
564        enddo
565        PRINT*, 'j=', j, coslat0 / cos(rlatu(j)), (1. + 2. * sum(wwf)) * cos(rlatu(j)) / coslat0
566      endif
567    enddo
568    CALL dump2d(imp1, jmp1, phis, 'phis ')
569    CALL dump2d(imp1, jmp1, masque, 'masque ')
570    CALL dump2d(imp1, jmp1, phis - phiso, 'dphis ')
571
572  END SUBROUTINE filtreoro
573
574
575END MODULE etat0phys
Note: See TracBrowser for help on using the repository browser.