GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/regr_horiz_time_climoz_m.F90 Lines: 0 288 0.0 %
Date: 2023-06-30 12:51:15 Branches: 0 1017 0.0 %

Line Branch Exec Source
1
MODULE 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_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,  nf95_get_var,  &
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
26
CONTAINS
27
28
!-------------------------------------------------------------------------------
29
!
30
SUBROUTINE 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, ALLOCATABLE:: lon_in(:)           ! Longitudes   (ascending order, rad)
91
  REAL, ALLOCATABLE:: lat_in(:)           ! Latitudes    (ascending order, rad)
92
  REAL, ALLOCATABLE:: 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, ALLOCATABLE :: 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
        CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
336
        IF(l3D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
337
        IF(l2D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
338
        IF(lprev) THEN; sta(ndims)=12
339
          CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
340
          IF(l3D) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
341
          IF(l2d) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
342
        END IF
343
        IF(lnext) THEN; sta(ndims)=1
344
          CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
345
          IF(l3D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
346
          IF(l2D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
347
        END IF
348
      END DO
349
      IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
350
      IF(lprev) CALL NF95_CLOSE(fID_in_m)
351
      IF(lnext) CALL NF95_CLOSE(fID_in_p)
352
353
      !--- Revert decreasing coordinates vector
354
      IF(l3D) THEN
355
        IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
356
        IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
357
        IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
358
359
        IF (grid_type /= unstructured) THEN
360
          !--- Shift values for longitude and duplicate some longitudes slices
361
          o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
362
          o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
363
        ENDIF
364
      ELSE
365
        IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
366
        IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
367
      END IF
368
369
     !--- Deal with missing values
370
      DO m=1, read_climoz
371
        WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
372
        IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
373
          IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
374
            WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
375
          END IF
376
        END IF
377
        WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
378
        WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
379
380
        !--- Check top layer contains no NaNs & search NaNs from top to ground
381
        msg=TRIM(sub)//": NaNs in top layer !"
382
        IF(l3D) THEN
383
          IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
384
          DO k = 2,nlev_in
385
            WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
386
          END DO
387
        ELSE
388
          IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
389
            WRITE(lunout,*)msg
390
            !--- Fill in latitudes where all values are missing
391
            DO l=1,nmth_in
392
              !--- Next to south pole
393
              j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
394
              IF(j>1) &
395
                o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
396
              !--- Next to north pole
397
              j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
398
              IF(j<nlat_in) &
399
                o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
400
            END DO
401
          END IF
402
403
          !--- Fill in high latitudes missing values
404
          !--- Highest level been filled-in, so has always valid values.
405
          DO k = 2,nlev_in
406
            WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
407
          END DO
408
        END IF
409
      END DO
410
411
    ENDIF
412
413
    !=============================================================================
414
    IF(l3D) THEN                                                   !=== 3D FIELDS
415
    !=============================================================================
416
     IF (grid_type==unstructured) THEN
417
#ifdef CPP_XIOS
418
       nlat_ou=klon_mpi
419
420
       IF (is_mpi_root) THEN
421
         ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz))
422
         o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
423
       ELSE
424
         ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
425
       ENDIF
426
       ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
427
428
       CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
429
       CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
430
#endif
431
     ELSE
432
433
       !--- Regrid in longitude
434
        ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
435
        CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
436
                            xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
437
                            vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
438
        DEALLOCATE(o3_in3)
439
440
        !--- Regrid in latitude: averaging with respect to SIN(lat) is
441
        !                        equivalent to weighting by COS(lat)
442
        !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
443
        ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
444
        CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
445
                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
446
                        vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
447
                   slope = slopes(2,o3_regr_lon, sinlat_in_edge))
448
        DEALLOCATE(o3_regr_lon)
449
450
     ENDIF
451
452
     !--- Duplicate previous/next record(s) if they are not available
453
     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
454
     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
455
456
     !--- Regrid in time by linear interpolation:
457
     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
458
     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
459
     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
460
     DEALLOCATE(o3_regr_lonlat)
461
462
     nlat_ou=nbp_lat
463
     IF (grid_type==unstructured) THEN
464
#ifdef CPP_XIOS
465
       CALL xios_send_field('o3_out',o3_out3)
466
       ndims=3
467
       ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
468
       CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
469
#endif
470
     ENDIF
471
472
    !--- Create the output file and get the variable IDs:
473
    CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
474
                     ndims, cal_ou)
475
476
    IF (is_mpi_root) THEN
477
      !--- Write remaining coordinate variables:
478
      CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
479
      IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
480
      IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
481
482
      !--- Write to file (the order of "rlatu" is inverted in the output file):
483
        IF (grid_type==unstructured) THEN
484
485
          ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
486
          DO i=1,klon_glo
487
            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
488
          ENDDO
489
490
          DO m = 1, read_climoz
491
            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
492
          END DO
493
494
        ELSE
495
          DO m = 1, read_climoz
496
            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
497
          END DO
498
      ENDIF
499
      CALL NF95_CLOSE(fID_ou)
500
501
502
    ENDIF
503
504
505
    !=============================================================================
506
    ELSE                                                         !=== ZONAL FIELDS
507
    !=============================================================================
508
509
     IF (grid_type==unstructured) THEN
510
#ifdef CPP_XIOS
511
       nlat_ou=klon_mpi
512
513
       IF (is_mpi_root) THEN
514
         ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz))
515
         o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
516
       ELSE
517
         ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
518
       ENDIF
519
       ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
520
       CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
521
       CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
522
       IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
523
       IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
524
#endif
525
526
     ELSE
527
        !--- Regrid in latitude: averaging with respect to SIN(lat) is
528
        !                        equivalent to weighting by COS(lat)
529
        !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
530
        ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
531
        CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
532
                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
533
                        vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
534
                     slope = slopes(1,o3_in2, sinlat_in_edge))
535
        DEALLOCATE(o3_in2)
536
537
        !--- Duplicate previous/next record(s) if they are not available
538
        IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
539
        IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
540
541
     ENDIF
542
543
      !--- Regrid in time by linear interpolation:
544
      ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
545
      IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
546
      IF(.NOT.interpt) o3_out2=o3_regr_lat
547
      DEALLOCATE(o3_regr_lat)
548
549
      nlat_ou=nbp_lat
550
551
      IF (grid_type==unstructured) THEN
552
        ndims=3
553
        ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
554
        CALL gather_mpi(o3_out2, o3_out2_glo)
555
      ENDIF
556
557
      !--- Create the output file and get the variable IDs:
558
      CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
559
                         ndims, cal_ou)
560
561
      IF (is_mpi_root) THEN
562
563
        !--- Write remaining coordinate variables:
564
        CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
565
        IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
566
        IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
567
568
        IF (grid_type==unstructured) THEN
569
570
          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
571
          DO i=1,klon_glo
572
            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
573
          ENDDO
574
575
576
          DO m = 1, read_climoz
577
            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
578
          END DO
579
        ELSE
580
          !--- Write to file (the order of "rlatu" is inverted in the output file):
581
          DO m = 1, read_climoz
582
            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
583
          END DO
584
        ENDIF
585
586
        CALL NF95_CLOSE(fID_ou)
587
588
      ENDIF
589
590
    !=============================================================================
591
    END IF
592
    !=============================================================================
593
594
    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
595
596
  ENDIF ! is_omp_master
597
598
  first=.FALSE.
599
END SUBROUTINE regr_horiz_time_climoz
600
!
601
!-------------------------------------------------------------------------------
602
603
604
!-------------------------------------------------------------------------------
605
!
606
SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
607
                       vID_ou, ndims, cal_ou)
608
!-------------------------------------------------------------------------------
609
! Purpose:  This subroutine creates the NetCDF output file, defines
610
!     dimensions and variables, and writes some of the coordinate variables.
611
!-------------------------------------------------------------------------------
612
  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
613
  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
614
  USE mod_phys_lmdz_para, ONLY: is_mpi_root
615
  USE mod_grid_phy_lmdz, ONLY: klon_glo
616
!
617
!-------------------------------------------------------------------------------
618
! Arguments:
619
  INTEGER, INTENT(IN)  :: fID_in, nlev_in, ntim_ou
620
  INTEGER, INTENT(OUT) :: fID_ou, vlevID,  vtimID
621
  INTEGER, INTENT(OUT) :: vID_ou(:)      ! dim(1/2) 1: O3day&night 2: O3daylight
622
  INTEGER, INTENT(IN)  :: ndims          ! fields rank (3 or 4)
623
  CHARACTER(LEN=*), INTENT(IN) :: cal_ou ! calendar
624
!-------------------------------------------------------------------------------
625
! Local variables:
626
  INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
627
  INTEGER :: vlonID, vlatID, ncerr,  is
628
  REAL,ALLOCATABLE    :: latitude_glo_(:)
629
  CHARACTER(LEN=80) :: sub
630
  INTEGER :: i
631
632
633
!-------------------------------------------------------------------------------
634
635
  IF (is_mpi_root) THEN
636
    sub="prepare_out"
637
    WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
638
    CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
639
640
  !--- Dimensions:
641
    IF(ndims==4) &
642
    CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
643
    CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
644
    CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
645
    CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
646
647
    !--- Define coordinate variables:
648
    IF(ndims==4) &
649
    CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
650
    CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
651
    CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
652
    CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
653
    IF(ndims==4) &
654
    CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
655
    CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
656
    CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
657
    CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
658
    IF(ndims==4) &
659
    CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
660
    CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
661
    CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
662
    CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
663
    CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
664
    CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
665
666
  !--- Define the main variables:
667
    IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
668
    IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
669
    CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
670
    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
671
    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
672
      &_in_air")
673
    IF(SIZE(vID_ou) == 2) THEN
674
      CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
675
      CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
676
        &ylight")
677
    END IF
678
679
  !--- Global attributes:
680
  ! The following commands, copying attributes, may fail. That is OK.
681
  ! It should just mean that the attribute is not defined in the input file.
682
    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
683
    CALL handle_err_copy_att("Conventions")
684
    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
685
    CALL handle_err_copy_att("title")
686
    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
687
    CALL handle_err_copy_att("institution")
688
    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
689
    CALL handle_err_copy_att("source")
690
    CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
691
    CALL NF95_ENDDEF(fID_ou)
692
693
    IF (grid_type==unstructured) THEN
694
      ALLOCATE(latitude_glo_(klon_glo))
695
      DO i=1,klon_glo
696
        latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
697
      ENDDO
698
      CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
699
    ELSE
700
      !--- Write one of the coordinate variables:
701
      IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
702
      CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
703
    !    (convert from rad to degrees and sort in ascending order)
704
    ENDIF
705
  ENDIF
706
707
CONTAINS
708
709
!-------------------------------------------------------------------------------
710
!
711
SUBROUTINE handle_err_copy_att(att_name)
712
!
713
!-------------------------------------------------------------------------------
714
  USE netcdf, ONLY: NF90_NOERR, NF90_strerror
715
!-------------------------------------------------------------------------------
716
! Arguments:
717
  CHARACTER(LEN=*), INTENT(IN) :: att_name
718
!-------------------------------------------------------------------------------
719
  IF(ncerr /= NF90_NOERR) &
720
    WRITE(lunout,*)TRIM(sub)//" prepare_out NF95_COPY_ATT "//TRIM(att_name)//  &
721
                      " -- "//TRIM(NF90_strerror(ncerr))
722
723
END SUBROUTINE handle_err_copy_att
724
!
725
!-------------------------------------------------------------------------------
726
727
END SUBROUTINE prepare_out
728
!
729
!-------------------------------------------------------------------------------
730
731
END MODULE regr_horiz_time_climoz_m
732
!
733
!-------------------------------------------------------------------------------