source: LMDZ6/trunk/libf/phylmd/regr_horiz_time_climoz_m.F90 @ 3741

Last change on this file since 3741 was 3465, checked in by Laurent Fairhead, 6 years ago

Further modifications for DYNAMICO/LMDZ convergence. These are based
on Yann's LMDZ6_V2 sources. Compiles on irene and converges with revision 3459
in a bucket configuration
YM/LF

File size: 32.6 KB
Line 
1MODULE regr_horiz_time_climoz_m
2
3  USE interpolation,     ONLY: locate
4  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
5  USE nrtype,            ONLY: pi
6  USE netcdf,   ONLY: NF90_CLOBBER, NF90_FLOAT,     NF90_GET_VAR, NF90_OPEN,   &
7                      NF90_NOWRITE, NF90_NOERR,     NF90_GET_ATT, NF90_GLOBAL
8  USE netcdf95, ONLY: NF95_DEF_DIM, NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION,    &
9                      NF95_DEF_VAR, NF95_INQ_VARID, NF95_INQUIRE_VARIABLE,     &
10          NF95_OPEN,  NF95_CREATE,  NF95_GET_ATT,   NF95_GW_VAR,  HANDLE_ERR,  &
11          NF95_CLOSE, NF95_ENDDEF,  NF95_PUT_ATT,   NF95_PUT_VAR, NF95_COPY_ATT
12  USE print_control_mod, ONLY: lunout
13  USE dimphy
14  IMPLICIT NONE
15  PRIVATE
16  PUBLIC :: regr_horiz_time_climoz
17  REAL, PARAMETER :: deg2rad=pi/180.
18  CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3         ','tro3_daylight']
19
20  INTEGER :: nlat_ou, nlon_ou
21  REAL, ALLOCATABLE :: latitude_glo(:)
22!$OMP THREADPRIVATE(latitude_glo)
23  INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
24!$OMP THREADPRIVATE(ind_cell_glo_glo)
25
26CONTAINS
27
28!-------------------------------------------------------------------------------
29!
30SUBROUTINE regr_horiz_time_climoz(read_climoz,interpt)
31!
32!-------------------------------------------------------------------------------
33! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
34!   * Read ozone climatology from netcdf file
35!   * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
36!   * If interpt=T, interpolate linearly in time (one record each day)
37!     If interpt=F, keep original time sampling  (14 months).
38!   * Save it to a new netcdf file.
39!-------------------------------------------------------------------------------
40! Remarks:
41!   * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
42!   * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
43!   * Output grid cells centers coordinates given by [rlonv,] rlatu.
44!   * Output grid cells edges   coordinates given by [rlonu,] rlatv.
45!   * Input file [longitudes and] latitudes given in degrees.
46!   * Input file pressure levels are given in Pa or hPa.
47!   * All coordinates variables are stricly monotonic.
48!   * Monthly fields are interpolated linearly in time to get daily values.
49!   * Fields are known at the middle of the months, so interpolation requires an
50!     additional record both for 1st half of january and 2nd half of december:
51!     - For a 14-records "climoz.nc": records 1 and 14.
52!     - For 12-records files:
53!       record 12 of "climoz_m.nc" if available, or record 1  of "climoz.nc".
54!       record 1  of "climoz_p.nc" if available, or record 12 of "climoz.nc".
55!   * Calendar is taken into account to get one record each day (not 360 always).
56!   * Missing values are filled in from sky to ground by copying lowest valid one.
57!     Attribute "missing_value" or "_FillValue" must be present in input file.
58!-------------------------------------------------------------------------------
59  USE assert_m,           ONLY: assert
60  USE cal_tools_m,        ONLY: year_len, mid_month
61!!  USE control_mod,        ONLY: anneeref
62  USE time_phylmdz_mod,   ONLY: annee_ref
63  USE ioipsl,             ONLY: ioget_year_len, ioget_calendar
64  USE regr_conserv_m,     ONLY: regr_conserv
65  USE regr_lint_m,        ONLY: regr_lint
66  USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
67  USE slopes_m,           ONLY: slopes
68#ifdef CPP_XIOS
69  USE xios
70#endif
71  USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
72  USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
73  USE mod_grid_phy_lmdz, ONLY: klon_glo
74
75!-------------------------------------------------------------------------------
76! Arguments:
77  INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
78!                         1: read a single ozone climatology used day and night
79!                         2: same + read also a daylight climatology
80  LOGICAL, INTENT(IN) :: interpt     ! TRUE  => daily interpolation
81                                     ! FALSE => no interpolation (14 months)
82!-------------------------------------------------------------------------------
83! Local variables:
84
85!--- Input files variables
86  INTEGER :: nlon_in                       ! Number of longitudes
87  INTEGER :: nlat_in                       ! Number of latitudes
88  INTEGER :: nlev_in                       ! Number of pressure levels
89  INTEGER :: nmth_in                       ! Number of months
90  REAL, POINTER     :: lon_in(:)           ! Longitudes   (ascending order, rad)
91  REAL, POINTER     :: lat_in(:)           ! Latitudes    (ascending order, rad)
92  REAL, POINTER     :: lev_in(:)           ! Pressure levels (ascen. order, hPa)
93  REAL, ALLOCATABLE :: lon_in_edge(:)      ! Longitude intervals edges
94                                           !              (ascending order,  / )
95  REAL, ALLOCATABLE :: sinlat_in_edge(:)   ! Sinus of latitude intervals edges
96                                           !              (ascending order,  / )
97  LOGICAL :: ldec_lon, ldec_lat, ldec_lev  ! Decreasing order in input file
98  CHARACTER(LEN=20) :: cal_in              ! Calendar
99  REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:)   ! Ozone climatologies
100  REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:)   ! Ozone climatologies
101  REAL, ALLOCATABLE :: o3_in2  (:,:,:,:)   ! Ozone climatologies
102  REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:)   ! Ozone climatologies
103  ! last index: 1 for the day-night average, 2 for the daylight field.
104  REAL :: NaN
105
106!--- Partially or totally regridded variables      (:,:,nlev_in,:,read_climoz)
107  REAL, ALLOCATABLE :: o3_regr_lon   (:,:,:,:,:) ! (nlon_ou,nlat_in,:,0:13   ,:)
108  REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13   ,:)
109  REAL, ALLOCATABLE :: o3_out3       (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
110  REAL, ALLOCATABLE :: o3_out3_glo   (:,:,:,:) !   (nbp_lat,:,ntim_ou,:)
111  REAL, ALLOCATABLE :: o3_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
112  REAL, ALLOCATABLE :: o3_out2         (:,:,:,:) !         (nlat_ou,:,ntim_ou,:)
113  REAL, ALLOCATABLE :: o3_out2_glo     (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
114  REAL, ALLOCATABLE :: o3_out          (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
115! Dimension number  | Interval                | Contains  | For variables:
116!   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
117!   2 (latitude)    | [rlatv(j), rlatv(j-1)]  | rlatu(j)  | all but o3_regr_lon
118!   3 (press level) |                         |   lev(k)  | all
119! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
120! Dimension 4 is: month number                             (all vars but o3_out)
121!                 days elapsed since Jan. 1st 0h at mid-day (o3_out only)
122  REAL, ALLOCATABLE :: v1(:)
123
124!--- For NetCDF:
125  INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
126  INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
127  INTEGER, POINTER :: dIDs(:)
128  CHARACTER(LEN=20) :: cal_ou     !--- Calendar; no time inter => same as input
129  CHARACTER(LEN=80) :: press_unit !--- Pressure unit
130  REAL    :: tmidmonth(0:13)      !--- Elapsed days since Jan-1 0h at mid-months
131                                  ! Additional records 0, 13 for interpolation
132  REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
133  LOGICAL :: lprev, lnext         !--- Flags: previous/next files are present
134  LOGICAL :: l3D, l2D             !--- Flag:  input fields are 3D or zonal
135  INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
136  INTEGER, ALLOCATABLE :: sta(:), cnt(:)
137  CHARACTER(LEN=80) :: sub, dim_nam, msg
138  REAL :: null_array(0)
139  LOGICAL,SAVE :: first=.TRUE.
140!$OMP THREADPRIVATE(first) 
141  REAL, ALLOCATABLE :: test_o3_in(:,:)
142  REAL, ALLOCATABLE :: test_o3_out(:)
143
144
145  IF (grid_type==unstructured) THEN
146    IF (first) THEN
147      IF (is_master) THEN
148        ALLOCATE(latitude_glo(klon_glo))
149        ALLOCATE(ind_cell_glo_glo(klon_glo))
150      ELSE
151        ALLOCATE(latitude_glo(0))
152        ALLOCATE(ind_cell_glo_glo(0))
153      ENDIF
154      CALL gather(latitude_deg,  latitude_glo)
155      CALL gather(ind_cell_glo,  ind_cell_glo_glo)
156    ENDIF
157  ENDIF
158   
159  IF (is_omp_master) THEN
160    nlat_ou=nbp_lat
161    nlon_ou=nbp_lon
162   
163   !-------------------------------------------------------------------------------
164    IF (is_mpi_root) THEN
165      sub="regr_horiz_time_climoz"
166      WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
167      CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
168
169      CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
170      lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
171      lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
172
173      !--- Get coordinates from the input file. Converts lon/lat in radians.
174      !    Few inversions because "regr_conserv" and gcm need ascending vectors.
175      CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
176      CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
177      l3D=ndims==4; l2D=ndims==3
178      IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
179      IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
180      DO i=1,ndims
181        CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
182        CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
183        ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE
184        SELECT CASE(ii)
185          CASE(1)                                         !--- LONGITUDE
186            CALL NF95_GW_VAR(fID_in, varid, lon_in)
187            ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
188            nlon_in=dln; lon_in=lon_in*deg2rad
189          CASE(2)                                         !--- LATITUDE
190            CALL NF95_GW_VAR(fID_in, varid, lat_in)
191            ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
192            nlat_in=dln; lat_in=lat_in*deg2rad
193          CASE(3)                                         !--- PRESSURE LEVELS
194            CALL NF95_GW_VAR(fID_in, varid, lev_in)
195            ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
196            nlev_in=dln
197            CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
198            k=LEN_TRIM(press_unit)
199            DO WHILE(ICHAR(press_unit(k:k))==0)
200              press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
201            END DO
202            IF(press_unit ==  "Pa") THEN
203              lev_in = lev_in/100.                        !--- CONVERT TO hPa
204            ELSE IF(press_unit /= "hPa") THEN
205              CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
206            END IF
207          CASE(4)                                         !--- TIME
208            CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
209            cal_in='gregorian'
210            IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        & 
211            WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
212            TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
213            k=LEN_TRIM(cal_in)
214            DO WHILE(ICHAR(cal_in(k:k))==0)
215              cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
216            END DO
217        END SELECT
218      END DO
219
220      !--- Prepare quantities for time interpolation
221      tmidmonth=mid_month(annee_ref, cal_in)
222      IF(interpt) THEN
223        ntim_ou=ioget_year_len(annee_ref)
224        ALLOCATE(tmidday(ntim_ou))
225        tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
226        CALL ioget_calendar(cal_ou)
227      ELSE
228        ntim_ou=14
229        cal_ou=cal_in
230      END IF
231    ENDIF
232
233    IF (grid_type==unstructured) THEN
234      CALL bcast_mpi(nlon_in)
235      CALL bcast_mpi(nlat_in)
236      CALL bcast_mpi(nlev_in)
237      CALL bcast_mpi(l3d)
238      CALL bcast_mpi(tmidmonth)
239      CALL bcast_mpi(tmidday)
240      CALL bcast_mpi(ntim_ou)
241
242#ifdef CPP_XIOS   
243      IF (is_mpi_root) THEN
244        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
245        IF (l3D) THEN
246          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
247        ELSE
248          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
249        ENDIF
250      ELSE
251        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
252        IF (l3D) THEN
253          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
254        ELSE
255          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
256        ENDIF
257      ENDIF
258      CALL  xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
259      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
260      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
261      CALL  xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
262      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
263      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
264#endif
265     
266      IF (first) THEN
267        first=.FALSE.
268        RETURN
269      ENDIF
270    ENDIF
271   
272   
273    IF (is_mpi_root) THEN     
274      !--- Longitudes management:
275      !    * Need to shift data if the origin of input file longitudes /= -pi
276      !    * Need to add some margin in longitude to ensure input interval contains
277      !      all the output intervals => at least one longitudes slice has to be
278      !      duplicated, possibly more for undersampling.
279      IF(l3D) THEN
280        IF (grid_type==unstructured) THEN
281          dx2=0
282        ELSE
283          !--- Compute input edges longitudes vector (no end point yet)
284          ALLOCATE(v1(nlon_in+1))
285          v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
286          FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
287          v1(nlon_in+1)=v1(1)+2.*pi
288          DEALLOCATE(lon_in)
289
290          !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
291          v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
292
293          !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
294          dx1=locate(v1,boundslon_reg(1,west))-1
295          v1=CSHIFT(v1,SHIFT=dx1,DIM=1)
296          v1(nlon_in-dx1+2:)=v1(nlon_in-dx1+2:)+2.*pi
297   
298          !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
299          dx2=0; DO WHILE(v1(1+dx2)+2.*pi<=boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
300
301          !--- Final edges longitudes vector (with margin and end point)
302          ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
303          DEALLOCATE(v1)
304        ENDIF
305      END IF
306
307      !--- Compute sinus of intervals edges latitudes:
308      ALLOCATE(sinlat_in_edge(nlat_in+1))
309      sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
310      FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
311      DEALLOCATE(lat_in)
312
313
314
315      !--- Check for contiguous years:
316      ib=0; ie=13
317      IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
318        WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
319      ELSE 
320        IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
321        IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
322        IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
323        IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
324        IF(.NOT.lprev) ib=1
325        IF(.NOT.lnext) ie=12
326      END IF
327      ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1 
328      IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
329      IF(l2D) cnt=[        nlat_in,nlev_in,1] 
330      IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
331      IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
332
333      !--- Read full current file and one record each available contiguous file
334      DO iv=1,read_climoz
335        msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
336        CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
337        IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
338        IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
339        CALL handle_err(TRIM(msg), ncerr, fID_in)
340        IF(lprev) THEN; sta(ndims)=12 
341          CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
342          IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
343          IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
344          CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
345        END IF
346        IF(lnext) THEN; sta(ndims)=1 
347          CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
348          IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
349          IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
350          CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
351        END IF
352      END DO
353      IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
354      IF(lprev) CALL NF95_CLOSE(fID_in_m)
355      IF(lnext) CALL NF95_CLOSE(fID_in_p)
356
357      !--- Revert decreasing coordinates vector
358      IF(l3D) THEN
359        IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
360        IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
361        IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
362       
363        IF (grid_type /= unstructured) THEN
364          !--- Shift values for longitude and duplicate some longitudes slices
365          o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
366          o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
367        ENDIF
368      ELSE
369        IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
370        IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
371      END IF
372
373     !--- Deal with missing values
374      DO m=1, read_climoz
375        WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
376        IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
377          IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
378            WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
379          END IF
380        END IF
381        WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
382        WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
383
384        !--- Check top layer contains no NaNs & search NaNs from top to ground
385        msg=TRIM(sub)//": NaNs in top layer !"
386        IF(l3D) THEN
387          IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
388          DO k = 2,nlev_in
389            WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
390          END DO
391        ELSE
392          IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
393            WRITE(lunout,*)msg
394            !--- Fill in latitudes where all values are missing
395            DO l=1,nmth_in
396              !--- Next to south pole
397              j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
398              IF(j>1) &
399                o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
400              !--- Next to north pole
401              j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
402              IF(j<nlat_in) &
403                o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
404            END DO
405          END IF
406
407          !--- Fill in high latitudes missing values
408          !--- Highest level been filled-in, so has always valid values.
409          DO k = 2,nlev_in
410            WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
411          END DO
412        END IF
413      END DO
414
415    ENDIF
416   
417    !=============================================================================
418    IF(l3D) THEN                                                   !=== 3D FIELDS
419    !=============================================================================
420     IF (grid_type==unstructured) THEN
421#ifdef CPP_XIOS
422       nlat_ou=klon_mpi
423       
424       IF (is_mpi_root) THEN
425         ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz))
426         o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
427       ELSE
428         ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
429       ENDIF
430       ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
431       
432       CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
433       CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
434#endif
435     ELSE
436         
437       !--- Regrid in longitude
438        ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
439        CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
440                            xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
441                            vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
442        DEALLOCATE(o3_in3)
443
444        !--- Regrid in latitude: averaging with respect to SIN(lat) is
445        !                        equivalent to weighting by COS(lat)
446        !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
447        ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
448        CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
449                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
450                        vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
451                   slope = slopes(2,o3_regr_lon, sinlat_in_edge))
452        DEALLOCATE(o3_regr_lon)
453
454     ENDIF
455
456     !--- Duplicate previous/next record(s) if they are not available
457     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
458     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
459     
460     !--- Regrid in time by linear interpolation:
461     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
462     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
463     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
464     DEALLOCATE(o3_regr_lonlat)
465
466     nlat_ou=nbp_lat
467     IF (grid_type==unstructured) THEN
468#ifdef CPP_XIOS
469       CALL xios_send_field('o3_out',o3_out3)
470       ndims=3
471       ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
472       CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
473#endif
474     ENDIF
475
476    !--- Create the output file and get the variable IDs:
477    CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
478                     ndims, cal_ou)
479
480    IF (is_mpi_root) THEN
481      !--- Write remaining coordinate variables:
482      CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
483      IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
484      IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
485
486      !--- Write to file (the order of "rlatu" is inverted in the output file):
487        IF (grid_type==unstructured) THEN
488
489          ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
490          DO i=1,klon_glo
491            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
492          ENDDO
493
494          DO m = 1, read_climoz
495            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
496          END DO
497         
498        ELSE
499          DO m = 1, read_climoz
500            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
501          END DO
502      ENDIF
503      CALL NF95_CLOSE(fID_ou)
504
505
506    ENDIF
507
508
509    !=============================================================================
510    ELSE                                                         !=== ZONAL FIELDS
511    !=============================================================================
512   
513     IF (grid_type==unstructured) THEN
514#ifdef CPP_XIOS
515       nlat_ou=klon_mpi
516
517       IF (is_mpi_root) THEN
518         ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz))
519         o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
520       ELSE
521         ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
522       ENDIF
523       ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
524       CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
525       CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
526       IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
527       IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
528#endif       
529     
530     ELSE
531        !--- Regrid in latitude: averaging with respect to SIN(lat) is
532        !                        equivalent to weighting by COS(lat)
533        !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
534        ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
535        CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
536                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
537                        vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
538                     slope = slopes(1,o3_in2, sinlat_in_edge))
539        DEALLOCATE(o3_in2)
540
541        !--- Duplicate previous/next record(s) if they are not available
542        IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
543        IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
544
545     ENDIF
546     
547      !--- Regrid in time by linear interpolation:
548      ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
549      IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
550      IF(.NOT.interpt) o3_out2=o3_regr_lat
551      DEALLOCATE(o3_regr_lat)
552
553      nlat_ou=nbp_lat
554   
555      IF (grid_type==unstructured) THEN
556        ndims=3
557        ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
558        CALL gather_mpi(o3_out2, o3_out2_glo)
559      ENDIF
560     
561      !--- Create the output file and get the variable IDs:
562      CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
563                         ndims, cal_ou)
564
565      IF (is_mpi_root) THEN
566     
567        !--- Write remaining coordinate variables:
568        CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
569        IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
570        IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
571
572        IF (grid_type==unstructured) THEN
573
574          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
575          DO i=1,klon_glo
576            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
577          ENDDO
578
579
580          DO m = 1, read_climoz
581            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
582          END DO
583        ELSE
584          !--- Write to file (the order of "rlatu" is inverted in the output file):
585          DO m = 1, read_climoz
586            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
587          END DO
588        ENDIF
589       
590        CALL NF95_CLOSE(fID_ou)
591     
592      ENDIF
593
594    !=============================================================================
595    END IF
596    !=============================================================================
597
598    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
599
600  ENDIF ! is_omp_master
601
602  first=.FALSE.
603END SUBROUTINE regr_horiz_time_climoz
604!
605!-------------------------------------------------------------------------------
606
607
608!-------------------------------------------------------------------------------
609!
610SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
611                       vID_ou, ndims, cal_ou)
612!-------------------------------------------------------------------------------
613! Purpose:  This subroutine creates the NetCDF output file, defines
614!     dimensions and variables, and writes some of the coordinate variables.
615!-------------------------------------------------------------------------------
616  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
617  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
618  USE mod_phys_lmdz_para, ONLY: is_mpi_root
619  USE mod_grid_phy_lmdz, ONLY: klon_glo
620!
621!-------------------------------------------------------------------------------
622! Arguments:
623  INTEGER, INTENT(IN)  :: fID_in, nlev_in, ntim_ou
624  INTEGER, INTENT(OUT) :: fID_ou, vlevID,  vtimID
625  INTEGER, INTENT(OUT) :: vID_ou(:)      ! dim(1/2) 1: O3day&night 2: O3daylight
626  INTEGER, INTENT(IN)  :: ndims          ! fields rank (3 or 4)
627  CHARACTER(LEN=*), INTENT(IN) :: cal_ou ! calendar
628!-------------------------------------------------------------------------------
629! Local variables:
630  INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
631  INTEGER :: vlonID, vlatID, ncerr,  is
632  REAL,ALLOCATABLE    :: latitude_glo_(:)
633  CHARACTER(LEN=80) :: sub
634  INTEGER :: i
635
636
637!-------------------------------------------------------------------------------
638 
639  IF (is_mpi_root) THEN 
640    sub="prepare_out"
641    WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
642    CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
643
644  !--- Dimensions:
645    IF(ndims==4) &
646    CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
647    CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
648    CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
649    CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
650
651    !--- Define coordinate variables:
652    IF(ndims==4) &
653    CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
654    CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
655    CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
656    CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
657    IF(ndims==4) &
658    CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
659    CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
660    CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
661    CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
662    IF(ndims==4) &
663    CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
664    CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
665    CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
666    CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
667    CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
668    CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
669
670  !--- Define the main variables:
671    IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
672    IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
673    CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
674    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
675    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
676      &_in_air")
677    IF(SIZE(vID_ou) == 2) THEN
678      CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
679      CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
680        &ylight")
681    END IF
682
683  !--- Global attributes:
684  ! The following commands, copying attributes, may fail. That is OK.
685  ! It should just mean that the attribute is not defined in the input file.
686    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
687    CALL handle_err_copy_att("Conventions")
688    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
689    CALL handle_err_copy_att("title")
690    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
691    CALL handle_err_copy_att("institution")
692    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
693    CALL handle_err_copy_att("source")
694    CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
695    CALL NF95_ENDDEF(fID_ou)
696
697    IF (grid_type==unstructured) THEN
698      ALLOCATE(latitude_glo_(klon_glo))
699      DO i=1,klon_glo
700        latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
701      ENDDO
702      CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
703    ELSE
704      !--- Write one of the coordinate variables:
705      IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
706      CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
707    !    (convert from rad to degrees and sort in ascending order)
708    ENDIF
709  ENDIF
710 
711CONTAINS
712
713!-------------------------------------------------------------------------------
714!
715SUBROUTINE handle_err_copy_att(att_name)
716!
717!-------------------------------------------------------------------------------
718  USE netcdf, ONLY: NF90_NOERR, NF90_strerror
719!-------------------------------------------------------------------------------
720! Arguments:
721  CHARACTER(LEN=*), INTENT(IN) :: att_name
722!-------------------------------------------------------------------------------
723  IF(ncerr /= NF90_NOERR) &
724    WRITE(lunout,*)TRIM(sub)//" prepare_out NF95_COPY_ATT "//TRIM(att_name)//  &
725                      " -- "//TRIM(NF90_strerror(ncerr))
726
727END SUBROUTINE handle_err_copy_att
728!
729!-------------------------------------------------------------------------------
730
731END SUBROUTINE prepare_out
732!
733!-------------------------------------------------------------------------------
734
735END MODULE regr_horiz_time_climoz_m
736!
737!-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.