source: LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 @ 2942

Last change on this file since 2942 was 2942, checked in by dcugnet, 7 years ago

Remove test on sea ice concentration.

The aim of the test was to detect if a sea ice concentration field supposedly in
percents (holding the attribute "units" as correct) was in fact in area fraction.
But in some cases, all the fractions can be lower than one on the globe for
particular months, even if the field is really in percents.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 33.1 KB
Line 
1MODULE limit
2!
3!*******************************************************************************
4! Author : L. Fairhead, 27/01/94
5!-------------------------------------------------------------------------------
6! Purpose: Boundary conditions files building for new model using climatologies.
7!          Both grids have to be regular.
8!-------------------------------------------------------------------------------
9! Note: This routine is designed to work for Earth
10!-------------------------------------------------------------------------------
11! Modification history:
12!  * 23/03/1994: Z. X. Li
13!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
14!  *    07/2001: P. Le Van
15!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
16!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
17!-------------------------------------------------------------------------------
18
19  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
20  USE assert_eq_m,        ONLY: assert_eq
21  USE cal_tools_m,        ONLY: year_len, mid_month
22  USE conf_dat_m,         ONLY: conf_dat2d, conf_dat3d
23  USE dimphy,             ONLY: klon, zmasq
24  USE geometry_mod,       ONLY: longitude_deg, latitude_deg
25  USE phys_state_var_mod, ONLY: pctsrf
26  USE control_mod,        ONLY: anneeref
27  USE init_ssrf_m,        ONLY: start_init_subsurf
28
29  CHARACTER(LEN=20), PARAMETER :: &
30  fsst(5)=['amipbc_sst_1x1.nc   ','amip_sst_1x1.nc     ','cpl_atm_sst.nc      '&
31          ,'histmth_sst.nc      ','sstk.nc             ']
32  CHARACTER(LEN=20), PARAMETER :: &
33  fsic(5)=['amipbc_sic_1x1.nc   ','amip_sic_1x1.nc     ','cpl_atm_sic.nc      '&
34          ,'histmth_sic.nc      ','ci.nc               ']
35  CHARACTER(LEN=10), PARAMETER :: &
36  vsst(5)=['tosbcs    ','tos       ','SISUTESW  ','tsol_oce  ','sstk      '],  &
37  vsic(5)=['sicbcs    ','sic       ','SIICECOV  ','pourc_sic ','ci        ']
38  CHARACTER(LEN=10), PARAMETER :: &
39  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',    &
40   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    '
41
42CONTAINS
43
44!-------------------------------------------------------------------------------
45!
46SUBROUTINE limit_netcdf(masque, phis, extrap)
47!
48!-------------------------------------------------------------------------------
49! Author : L. Fairhead, 27/01/94
50!-------------------------------------------------------------------------------
51! Purpose: Boundary conditions files building for new model using climatologies.
52!          Both grids have to be regular.
53!-------------------------------------------------------------------------------
54! Note: This routine is designed to work for Earth
55!-------------------------------------------------------------------------------
56! Modification history:
57!  * 23/03/1994: Z. X. Li
58!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
59!  *    07/2001: P. Le Van
60!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
61!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
62!  *    04/2016: D. Cugnet   (12/14 recs SST/SIC files: cyclic/interannual runs)
63!  *    05/2017: D. Cugnet   (linear time interpolation for BCS files)
64!-------------------------------------------------------------------------------
65#ifndef CPP_1D
66  USE indice_sol_mod
67  USE netcdf,             ONLY: NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,        &
68                  NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,      &
69                  NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,       &
70                  NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
71  USE inter_barxy_m,      ONLY: inter_barxy
72  USE netcdf95,           ONLY: nf95_def_var, nf95_put_att, nf95_put_var
73  USE comconst_mod, ONLY: pi
74  USE phys_cal_mod, ONLY: calend
75  IMPLICIT NONE
76!-------------------------------------------------------------------------------
77! Arguments:
78  include "iniprint.h"
79  include "dimensions.h"
80  include "paramet.h"
81  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
82  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis   ! ground geopotential
83  LOGICAL,                    INTENT(IN)    :: extrap ! SST extrapolation flag
84!-------------------------------------------------------------------------------
85! Local variables:
86  include "comgeom2.h"
87
88!--- INPUT NETCDF FILES NAMES --------------------------------------------------
89  CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
90  CHARACTER(LEN=10) :: varname
91
92!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
93  REAL               :: fi_ice(klon), verif(klon)
94  REAL, POINTER      :: phy_rug(:,:)=>NULL(), phy_ice(:,:)=>NULL()
95  REAL, POINTER      :: phy_sst(:,:)=>NULL(), phy_alb(:,:)=>NULL()
96  REAL, ALLOCATABLE  :: phy_bil(:,:), pctsrf_t(:,:,:)
97  INTEGER            :: nbad
98
99!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
100  INTEGER :: ierr, nid, ndim, ntim, k, dims(2), ix_sic, ix_sst
101  INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
102  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC, varid_longitude, varid_latitude
103  INTEGER :: NF90_FORMAT
104  INTEGER :: ndays                   !--- Depending on the output calendar
105  CHARACTER(LEN=256) :: str
106
107!--- INITIALIZATIONS -----------------------------------------------------------
108#ifdef NC_DOUBLE
109  NF90_FORMAT=NF90_DOUBLE
110#else
111  NF90_FORMAT=NF90_FLOAT
112#endif
113  CALL inigeom
114
115!--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
116   IF(ALL(masque==-99999.)) THEN
117    CALL start_init_orog0(rlonv,rlatu,phis,masque)
118    CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq)          !--- To physical grid
119    ALLOCATE(pctsrf(klon,nbsrf))
120    CALL start_init_subsurf(.FALSE.)
121  !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
122    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
123    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
124  END IF
125
126!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
127  ndays=year_len(anneeref)
128
129!--- RUGOSITY TREATMENT --------------------------------------------------------
130  CALL msg(1,'Traitement de la rugosite')
131  CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
132
133!--- OCEAN TREATMENT -----------------------------------------------------------
134  CALL msg(1,'Traitement de la glace oceanique')
135
136! Input SIC file selection
137! Open file only to test if available
138  DO ix_sic=1,SIZE(fsic)
139     IF ( NF90_OPEN(TRIM(fsic(ix_sic)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
140        icefile=fsic(ix_sic); varname=vsic(ix_sic); EXIT
141     END IF
142  END DO
143  IF(ix_sic==SIZE(fsic)+1) THEN
144     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
145     WRITE(lunout,*) 'One of following files must be available : '
146     DO k=1,SIZE(fsic); WRITE(lunout,*) TRIM(fsic(k)); END DO
147     CALL abort_physic('limit_netcdf','No sea-ice file was found',1)
148  END IF
149  CALL ncerr(NF90_CLOSE(nid),icefile)
150  CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
151
152  CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
153
154  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
155  DO k=1,ndays
156     fi_ice=phy_ice(:,k)
157     WHERE(fi_ice>=1.0  ) fi_ice=1.0
158     WHERE(fi_ice<EPSFRA) fi_ice=0.0
159     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
160     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
161     SELECT CASE(ix_sic)
162        CASE(3)                                   ! SIC=pICE*(1-LIC-TER) (CPL)
163        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
164        CASE(4)                                   ! SIC=pICE            (HIST)
165        pctsrf_t(:,is_sic,k)=fi_ice(:)
166        CASE DEFAULT                              ! SIC=pICE-LIC   (AMIP,ERAI)
167        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
168     END SELECT
169     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
170     WHERE(1.0-zmasq<EPSFRA)
171        pctsrf_t(:,is_sic,k)=0.0
172        pctsrf_t(:,is_oce,k)=0.0
173     ELSEWHERE
174        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
175           pctsrf_t(:,is_sic,k)=1.0-zmasq
176           pctsrf_t(:,is_oce,k)=0.0
177        ELSEWHERE
178           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
179           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
180              pctsrf_t(:,is_oce,k)=0.0
181              pctsrf_t(:,is_sic,k)=1.0-zmasq
182           END WHERE
183        END WHERE
184     END WHERE
185     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
186     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb points = ',nbad
187     nbad=COUNT(ABS(SUM(pctsrf_t(:,:,k),DIM=2)-1.0)>EPSFRA)
188     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
189  END DO
190  DEALLOCATE(phy_ice)
191
192!--- SST TREATMENT -------------------------------------------------------------
193  CALL msg(1,'Traitement de la sst')
194
195! Input SST file selection
196! Open file only to test if available
197  DO ix_sst=1,SIZE(fsst)
198     IF ( NF90_OPEN(TRIM(fsst(ix_sst)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
199       sstfile=fsst(ix_sst); varname=vsst(ix_sst); EXIT
200     END IF
201  END DO
202  IF(ix_sst==SIZE(fsst)+1) THEN
203     WRITE(lunout,*) 'ERROR! No sst input file was found.'
204     WRITE(lunout,*) 'One of following files must be available : '
205     DO k=1,SIZE(fsst); WRITE(lunout,*) TRIM(fsst(k)); END DO
206     CALL abort_physic('limit_netcdf','No sst file was found',1)
207  END IF
208  CALL ncerr(NF90_CLOSE(nid),sstfile)
209  CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
210
211  CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
212
213!--- ALBEDO TREATMENT ----------------------------------------------------------
214  CALL msg(1,'Traitement de l albedo')
215  CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
216
217!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
218  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
219
220!--- OUTPUT FILE WRITING -------------------------------------------------------
221  CALL msg(5,'Ecriture du fichier limit : debut')
222  fnam="limit.nc"
223
224  !--- File creation
225  CALL ncerr(NF90_CREATE(fnam,NF90_CLOBBER,nid),fnam)
226  CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites"),fnam)
227  str='File produced using ce0l executable.'
228  str=TRIM(str)//NEW_LINE(' ')//'Sea Ice Concentration built from'
229  SELECT CASE(ix_sic)
230    CASE(1); str=TRIM(str)//' Amip mid-month boundary condition (BCS).'
231    CASE(2); str=TRIM(str)//' Amip monthly mean observations.'
232    CASE(3); str=TRIM(str)//' IPSL coupled model outputs.'
233    CASE(4); str=TRIM(str)//' LMDZ model outputs.'
234    CASE(5); str=TRIM(str)//' ci.nc file.'
235  END SELECT
236  str=TRIM(str)//NEW_LINE(' ')//'Sea Surface Temperature built from'
237  SELECT CASE(ix_sst)
238    CASE(1); str=TRIM(str)//' Amip mid-month boundary condition (BCS).'
239    CASE(2); str=TRIM(str)//' Amip monthly mean observations.'
240    CASE(3); str=TRIM(str)//' IPSL coupled model outputs.'
241    CASE(4); str=TRIM(str)//' LMDZ model outputs.'
242    CASE(5); str=TRIM(str)//' sstk.nc file.'
243  END SELECT
244  CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"history",TRIM(str)),fnam)
245
246  !--- Dimensions creation
247  CALL ncerr(NF90_DEF_DIM(nid,"points_physiques",klon,ndim),fnam)
248  CALL ncerr(NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim),fnam)
249
250  dims=[ndim,ntim]
251
252  !--- Variables creation
253  CALL ncerr(NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,[ntim],id_tim),fnam)
254  CALL ncerr(NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE),fnam)
255  CALL ncerr(NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC),fnam)
256  CALL ncerr(NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER),fnam)
257  CALL ncerr(NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC),fnam)
258  CALL ncerr(NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST),fnam)
259  CALL ncerr(NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS),fnam)
260  CALL ncerr(NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB),fnam)
261  CALL ncerr(NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG),fnam)
262  call nf95_def_var(nid, "longitude", NF90_FLOAT, ndim, varid_longitude)
263  call nf95_def_var(nid, "latitude",  NF90_FLOAT, ndim, varid_latitude)
264
265  !--- Attributes creation
266  CALL ncerr(NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee"),fnam)
267  CALL ncerr(NF90_PUT_ATT(nid,id_tim, "calendar",calend),fnam)
268  CALL ncerr(NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean"),fnam)
269  CALL ncerr(NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer"),fnam)
270  CALL ncerr(NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre"),fnam)
271  CALL ncerr(NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice"),fnam)
272  CALL ncerr(NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer"),fnam)
273  CALL ncerr(NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol"),fnam)
274  CALL ncerr(NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface"),fnam)
275  CALL ncerr(NF90_PUT_ATT(nid,id_RUG, "title","Rugosite"),fnam)
276
277  call nf95_put_att(nid, varid_longitude, "standard_name", "longitude")
278  call nf95_put_att(nid, varid_longitude, "units", "degrees_east")
279
280  call nf95_put_att(nid, varid_latitude, "standard_name", "latitude")
281  call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
282
283  CALL ncerr(NF90_ENDDEF(nid),fnam)
284
285  !--- Variables saving
286  CALL ncerr(NF90_PUT_VAR(nid,id_tim,[(REAL(k),k=1,ndays)]),fnam)
287  CALL ncerr(NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),[1,1],[klon,ndays]),fnam)
288  CALL ncerr(NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),[1,1],[klon,ndays]),fnam)
289  CALL ncerr(NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),[1,1],[klon,ndays]),fnam)
290  CALL ncerr(NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),[1,1],[klon,ndays]),fnam)
291  CALL ncerr(NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),[1,1],[klon,ndays]),fnam)
292  CALL ncerr(NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),[1,1],[klon,ndays]),fnam)
293  CALL ncerr(NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),[1,1],[klon,ndays]),fnam)
294  CALL ncerr(NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),[1,1],[klon,ndays]),fnam)
295  call nf95_put_var(nid, varid_longitude, longitude_deg)
296  call nf95_put_var(nid, varid_latitude, latitude_deg)
297
298  CALL ncerr(NF90_CLOSE(nid),fnam)
299
300  CALL msg(6,'Ecriture du fichier limit : fin')
301
302  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
303
304
305!===============================================================================
306!
307  CONTAINS
308!
309!===============================================================================
310
311
312!-------------------------------------------------------------------------------
313!
314SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask)
315!
316!-----------------------------------------------------------------------------
317! Comments:
318!   There are two assumptions concerning the NetCDF files, that are satisfied
319!   with files that are conforming NC convention:
320!     1) The last dimension of the variables used is the time record.
321!     2) Dimensional variables have the same names as corresponding dimensions.
322!-----------------------------------------------------------------------------
323  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
324       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
325       NF90_GET_ATT
326  USE pchsp_95_m, only: pchsp_95
327  USE pchfe_95_m, only: pchfe_95
328  USE arth_m, only: arth
329  USE indice_sol_mod
330
331  IMPLICIT NONE
332  include "dimensions.h"
333  include "paramet.h"
334  include "comgeom2.h"
335!-----------------------------------------------------------------------------
336! Arguments:
337  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
338  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
339  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
340  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
341  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
342  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
343  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
344!------------------------------------------------------------------------------
345! Local variables:
346!--- NetCDF
347  INTEGER           :: ncid, varid        ! NetCDF identifiers
348  CHARACTER(LEN=30) :: dnam               ! dimension name
349!--- dimensions
350  INTEGER           :: dids(4)            ! NetCDF dimensions identifiers
351  REAL, ALLOCATABLE :: dlon_ini(:)        ! initial longitudes vector
352  REAL, ALLOCATABLE :: dlat_ini(:)        ! initial latitudes  vector
353  REAL, POINTER     :: dlon(:), dlat(:)   ! reordered lon/lat  vectors
354!--- fields
355  INTEGER :: imdep, jmdep, lmdep          ! dimensions of 'champ'
356  REAL, ALLOCATABLE :: champ(:,:)         ! wanted field on initial grid
357  REAL, ALLOCATABLE :: yder(:), timeyear(:)
358  REAL              :: champint(iim,jjp1) ! interpolated field
359  REAL, ALLOCATABLE :: champtime(:,:,:)
360  REAL, ALLOCATABLE :: champan(:,:,:)
361!--- input files
362  CHARACTER(LEN=20) :: fnam_m, fnam_p     ! previous/next files names
363  CHARACTER(LEN=20) :: cal_in             ! calendar
364  CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
365  INTEGER           :: ndays_in           ! number of days
366!--- misc
367  INTEGER           :: i, j, k, l, ll     ! loop counters
368  REAL, ALLOCATABLE :: work(:,:)          ! used for extrapolation
369  CHARACTER(LEN=128):: title, mess        ! for messages
370  LOGICAL           :: is_bcs             ! flag for BCS data
371  LOGICAL           :: extrp              ! flag for extrapolation
372  REAL              :: chmin, chmax, timeday, al
373  INTEGER ierr, idx
374  integer n_extrap ! number of extrapolated points
375  logical skip
376
377!------------------------------------------------------------------------------
378!---Variables depending on keyword 'mode' -------------------------------------
379  NULLIFY(champo)
380
381  SELECT CASE(mode)
382  CASE('RUG'); title='Rugosite'
383  CASE('SIC'); title='Sea-ice'
384  CASE('SST'); title='SST'
385  CASE('ALB'); title='Albedo'
386  END SELECT
387  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
388  is_bcs=(mode=='SIC'.AND.ix_sic==1).OR.(mode=='SST'.AND.ix_sst==1)
389  idx=INDEX(fnam,'.nc')-1
390
391!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
392  CALL msg(5,' Now reading file : '//TRIM(fnam))
393  CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam)
394  CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
395  CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
396
397!--- Read unit for sea-ice variable only
398  IF (mode=='SIC') THEN
399    ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic); CALL strclean(unit_sic)
400    IF(ierr/=NF90_NOERR) THEN
401      CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
402      unit_sic='X'
403    ELSE IF(ANY(["%  ","1.0","1  "]==unit_sic)) THEN
404      CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
405    ELSE
406      CALL abort_physic('SIC','Unrecognized unit for sea-ice file: '&
407        //TRIM(unit_sic),1)
408    END IF
409  END IF
410
411!--- Longitude
412  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam)
413  ALLOCATE(dlon_ini(imdep), dlon(imdep))
414  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
415  CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam)
416  CALL msg(5,'variable '//TRIM(dnam)//' dimension ', imdep)
417
418!--- Latitude
419  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam)
420  ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
421  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
422  CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam)
423  CALL msg(5,'variable '//TRIM(dnam)//' dimension ', jmdep)
424
425!--- Time (variable is not needed - it is rebuilt - but calendar is)
426  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
427  ALLOCATE(timeyear(lmdep+2))
428  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
429  cal_in=' '
430  IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN
431    SELECT CASE(mode)
432      CASE('RUG', 'ALB'); cal_in='360d'
433      CASE('SIC', 'SST'); cal_in='gregorian'
434    END SELECT
435    CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
436     &//TRIM(fnam)//'. Choosing default value.')
437  END IF
438  CALL strclean(cal_in)                     !--- REMOVE (WEIRD) NULL CHARACTERS
439  CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
440 
441!--- Determining input file number of days, depending on calendar
442  ndays_in=year_len(anneeref, cal_in)
443
444!--- Rebuilding input time vector (field from input file might be unreliable)
445  IF(lmdep==12) THEN
446    timeyear=mid_month(anneeref, cal_in)
447    CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
448  ELSE IF(lmdep==ndays_in) THEN
449    timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
450    CALL msg(0,'Daily input file (no time interpolation).')
451  ELSE
452    WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
453      ' records, 12/',ndays_in,' (monthly/daily needed).'
454    CALL abort_physic('mid_month',TRIM(mess),1)
455  END IF
456
457!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
458  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
459  IF(extrp) ALLOCATE(work(imdep, jmdep))
460  CALL msg(5,'')
461  CALL msg(5,'READ AND INTERPOLATE HORIZONTALLY ', lmdep, ' FIELDS.')
462  CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
463  DO l=1, lmdep
464    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
465    !--- Check whether values are acceptable for SIC, depending on unit.
466    !--- Dropped for mid-month boundary conditions datasets (BCS, ix_sic==1)
467    IF(mode=='SIC'.AND.ix_sic/=1) THEN
468      IF(TRIM(unit_sic)=="1".OR.TRIM(unit_sic)=="1.0") THEN
469        IF(ANY(champ>1.0+EPSFRA)) &
470          CALL abort_physic('SIC','Found sea-ice fractions greater than 1.')
471      ELSE IF(TRIM(unit_sic)=="X".OR.TRIM(unit_sic)=="%") THEN
472        IF(ANY(champ>100.0+EPSFRA)) &
473          CALL abort_physic('SIC','Found sea-ice percentages greater than 100.')
474!        IF(MAXVAL(champ)< 1.01) &
475!          CALL abort_physic('SIC','All sea-ice percentages lower than 1.')
476      END IF
477    END IF
478    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
479    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
480    IF(l==1) THEN
481      CALL msg(5,"----------------------------------------------------------")
482      CALL msg(5,"$$$ Barycentrique interpolation for "//TRIM(title)//" $$$")
483      CALL msg(5,"----------------------------------------------------------")
484    END IF
485    IF(mode=='RUG') champ=LOG(champ)
486    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
487    IF(mode=='RUG') THEN
488      champint=EXP(champint)
489      WHERE(NINT(mask)/=1) champint=0.001
490    END IF
491    champtime(:, :, l+1)=champint
492  END DO
493  CALL ncerr(NF90_CLOSE(ncid), fnam)
494
495!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
496  fnam_m=fnam(1:idx)//'_m.nc'
497  IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
498    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
499    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
500    CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
501    CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m)
502    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
503    CALL ncerr(NF90_CLOSE(ncid), fnam_m)
504    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
505    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
506    IF(mode=='RUG') champ=LOG(champ)
507    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
508    IF(mode=='RUG') THEN
509      champint=EXP(champint)
510      WHERE(NINT(mask)/=1) champint=0.001
511    END IF
512    champtime(:, :, 1)=champint
513  ELSE
514    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
515    champtime(:, :, 1)=champtime(:, :, lmdep+1)
516  END IF
517
518!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
519  fnam_p=fnam(1:idx)//'_p.nc'
520  IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
521    CALL msg(0,'Reading next year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
522    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
523    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
524    CALL ncerr(NF90_CLOSE(ncid), fnam_p)
525    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
526    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
527    IF(mode=='RUG') champ=LOG(champ)
528    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
529    IF(mode=='RUG') THEN
530      champint=EXP(champint)
531      WHERE(NINT(mask)/=1) champint=0.001
532    END IF
533    champtime(:, :, lmdep+2)=champint
534  ELSE
535    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
536    champtime(:, :, lmdep+2)=champtime(:, :, 2)
537  END IF
538  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
539  IF(extrp) DEALLOCATE(work)
540
541!--- TIME INTERPOLATION ------------------------------------------------------
542  IF(prt_level>0) THEN
543     IF(ndays/=ndays_in) THEN
544        WRITE(lunout, *)
545        WRITE(lunout,*)'DIFFERENT YEAR LENGTHS:'
546        WRITE(lunout,*)' In the input file: ',ndays_in
547        WRITE(lunout,*)' In the output file: ',ndays
548     END IF
549     IF(lmdep==ndays_in) THEN
550        WRITE(lunout, *)
551        WRITE(lunout, *)'NO TIME INTERPOLATION.'
552        WRITE(lunout, *)' Daily input file.'
553     ELSE
554        IF(     is_bcs) WRITE(lunout, *)'LINEAR TIME INTERPOLATION.'
555        IF(.NOT.is_bcs) WRITE(lunout, *)'SPLINES TIME INTERPOLATION.'
556        WRITE(lunout, *)' Input time vector: ', timeyear
557        WRITE(lunout, *)' Output time vector from 0 to ', ndays-1
558     END IF
559  END IF
560  ALLOCATE(champan(iip1, jjp1, ndays))
561
562  IF(lmdep==ndays_in) THEN  !--- DAILY DATA: NO     TIME INTERPOLATION
563     champan(1:iim,:,:)=champtime
564  ELSE IF(is_bcs) THEN      !--- BCS   DATA: LINEAR TIME INTERPOLATION
565    l=1
566    DO k=1, ndays
567      timeday = (REAL(k)-0.5)*REAL(ndays_in)/ndays
568      IF(timeyear(l+1)<timeday) l=l+1
569      al=(timeday-timeyear(l))/(timeyear(l+1)-timeyear(l))
570      DO j=1, jjp1
571        DO i=1, iim
572          champan(i,j,k) = champtime(i,j,l)+al*(champtime(i,j,l+1)-champtime(i,j,l))
573        END DO
574      END DO
575    END DO
576  ELSE                      !--- AVE   DATA: SPLINE TIME INTERPOLATION
577     skip = .false.
578     n_extrap = 0
579     ALLOCATE(yder(lmdep+2))
580     DO j=1, jjp1
581       DO i=1, iim
582         yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
583              vc_beg=0., vc_end=0.)
584         CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
585              arth(0.5, real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
586         if (ierr < 0) stop 1
587         n_extrap = n_extrap + ierr
588       END DO
589     END DO
590     IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
591     DEALLOCATE(yder)
592  END IF
593  champan(iip1, :, :)=champan(1, :, :)
594  DEALLOCATE(champtime, timeyear)
595
596!--- Checking the result
597  DO j=1, jjp1
598    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
599    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' at time 10 ', chmin, chmax, j
600  END DO
601
602!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
603  IF(mode=='SST') THEN
604    CALL msg(0,'Filtering SST: SST >= 271.38')
605    WHERE(champan<271.38) champan=271.38
606  END IF
607
608!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
609  IF(mode=='SIC') THEN
610    CALL msg(0,'Filtering SIC: 0.0 < Sea-ice < 1.0')
611    IF(TRIM(unit_sic)=="1".OR.TRIM(unit_sic)=="1.0") THEN
612       ! Nothing to be done for sea-ice field is already in fraction of 1
613       ! This is the case for sea-ice in file cpl_atm_sic.nc
614       CALL msg(0,'Sea-ice field already in fraction of 1')
615    ELSE
616       ! Convert sea ice from percentage to fraction of 1
617       CALL msg(0,'Sea-ice field converted from percentage to fraction of 1.')
618       champan(:, :, :)=champan(:, :, :)/100.
619    END IF
620    champan(iip1, :, :)=champan(1, :, :)
621    WHERE(champan>1.0) champan=1.0
622    WHERE(champan<0.0) champan=0.0
623 END IF
624
625!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
626  ALLOCATE(champo(klon, ndays))
627  DO k=1, ndays
628    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
629  END DO
630  DEALLOCATE(champan)
631
632END SUBROUTINE get_2Dfield
633!
634!-------------------------------------------------------------------------------
635
636
637!-------------------------------------------------------------------------------
638!
639SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
640!
641!-------------------------------------------------------------------------------
642  USE grid_noro_m, ONLY: grid_noro0
643  IMPLICIT NONE
644!===============================================================================
645! Purpose:  Compute "phis" just like it would be in start_init_orog.
646!===============================================================================
647! Arguments:
648  REAL,             INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
649  REAL,             INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
650!-------------------------------------------------------------------------------
651! Local variables:
652  CHARACTER(LEN=256) :: modname="start_init_orog0"
653  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
654  REAL               :: lev(1), date, dt, deg2rad
655  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
656  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
657!-------------------------------------------------------------------------------
658  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
659  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
660  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
661  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
662  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
663  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
664
665!--- HIGH RESOLUTION OROGRAPHY
666  CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
667
668  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
669  CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,   &
670                lev, ttm_tmp, itau, date, dt, fid)
671  ALLOCATE(relief_hi(iml_rel,jml_rel))
672  CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
673  CALL flinclo(fid)
674
675!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
676  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
677  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
678  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
679
680!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
681  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
682  CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
683  DEALLOCATE(lon_ini,lat_ini)
684
685!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
686  WRITE(lunout,*)
687  WRITE(lunout,*)'*** Compute surface geopotential ***'
688
689!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
690  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
691  phis = phis * 9.81
692  phis(iml,:) = phis(1,:)
693  DEALLOCATE(relief_hi,lon_rad,lat_rad)
694
695END SUBROUTINE start_init_orog0
696!
697!-------------------------------------------------------------------------------
698
699
700!-------------------------------------------------------------------------------
701!
702SUBROUTINE msg(lev,str1,i,str2)
703!
704!-------------------------------------------------------------------------------
705! Arguments:
706  INTEGER,                    INTENT(IN) :: lev
707  CHARACTER(LEN=*),           INTENT(IN) :: str1
708  INTEGER,          OPTIONAL, INTENT(IN) :: i
709  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
710!-------------------------------------------------------------------------------
711  IF(prt_level>lev) THEN
712    IF(PRESENT(str2)) THEN
713      WRITE(lunout,*) TRIM(str1), i, TRIM(str2)
714    ELSE IF(PRESENT(i)) THEN
715      WRITE(lunout,*) TRIM(str1), i
716    ELSE
717      WRITE(lunout,*) TRIM(str1)
718    END IF
719  END IF
720
721END SUBROUTINE msg
722!
723!-------------------------------------------------------------------------------
724
725
726!-------------------------------------------------------------------------------
727!
728SUBROUTINE ncerr(ncres,fnam)
729!
730!-------------------------------------------------------------------------------
731! Purpose: NetCDF errors handling.
732!-------------------------------------------------------------------------------
733  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
734  IMPLICIT NONE
735!-------------------------------------------------------------------------------
736! Arguments:
737  INTEGER,          INTENT(IN) :: ncres
738  CHARACTER(LEN=*), INTENT(IN) :: fnam
739!-------------------------------------------------------------------------------
740  IF(ncres/=NF90_NOERR) THEN
741    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
742    CALL abort_physic('limit_netcdf',NF90_STRERROR(ncres),1)
743  END IF
744
745END SUBROUTINE ncerr
746!
747!-------------------------------------------------------------------------------
748
749
750!-------------------------------------------------------------------------------
751!
752SUBROUTINE strclean(s)
753!
754!-------------------------------------------------------------------------------
755  IMPLICIT NONE
756!-------------------------------------------------------------------------------
757! Purpose: Remove tail null characters from the input string.
758!-------------------------------------------------------------------------------
759! Parameters:
760  CHARACTER(LEN=*), INTENT(INOUT) :: s
761!-------------------------------------------------------------------------------
762! Local variable:
763  INTEGER :: k
764!-------------------------------------------------------------------------------
765  k=LEN_TRIM(s); DO WHILE(ICHAR(s(k:k))==0); s(k:k)=' '; k=LEN_TRIM(s); END DO
766
767END SUBROUTINE strclean
768!
769!-------------------------------------------------------------------------------
770
771#endif
772! of #ifndef CPP_1D
773END SUBROUTINE limit_netcdf
774
775END MODULE limit
776!
777!*******************************************************************************
778
Note: See TracBrowser for help on using the repository browser.