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

Last change on this file since 5153 was 5139, checked in by abarral, 8 weeks ago

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