source: LMDZ6/branches/Optimisation_LMDZ/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

Last change on this file was 3630, checked in by Laurent Fairhead, 5 years ago

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