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

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

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