source: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90 @ 3771

Last change on this file since 3771 was 3771, checked in by lguez, 4 years ago

Bug fix: initial value of surface humidity

Rename qsolsrf to qsurf as it should contain surface humidity. Fix the
value, which should be between 0 and 1. A better initial value would
be saturation humidity but we do not have yet the surface pressure in
etat0phys_netcdf, it will be defined later in etat0dyn_netcdf. The
old value, 150, probably came from a confusion with soil water
content. This initial value of surface humidity is used in the first
time step of physics to compute a virtual temperature at the surface,
in cdrag.
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)         :: qsurf, 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, qsurf, 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  qsurf = 0.
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, qsurf, 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.