Ignore:
Timestamp:
Nov 12, 2018, 1:52:29 PM (6 years ago)
Author:
Laurent Fairhead
Message:

Inclusion of Yann's latest (summer/fall 2018) modifications for
convergence of DYNAMICO/LMDZ physics
YM/LF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_horiz_time_climoz_m.F90

    r3411 r3413  
    22
    33  USE interpolation,     ONLY: locate
    4   USE mod_grid_phy_lmdz, ONLY: nlon_ou => nbp_lon, nlat_ou => nbp_lat
     4  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
    55  USE nrtype,            ONLY: pi
    66  USE netcdf,   ONLY: NF90_CLOBBER, NF90_FLOAT,     NF90_GET_VAR, NF90_OPEN,   &
     
    1111          NF95_CLOSE, NF95_ENDDEF,  NF95_PUT_ATT,   NF95_PUT_VAR, NF95_COPY_ATT
    1212  USE print_control_mod, ONLY: lunout
     13  USE dimphy
    1314  IMPLICIT NONE
    1415  PRIVATE
     
    1617  REAL, PARAMETER :: deg2rad=pi/180.
    1718  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)
    1825
    1926CONTAINS
     
    5259  USE assert_m,           ONLY: assert
    5360  USE cal_tools_m,        ONLY: year_len, mid_month
     61!!  USE control_mod,        ONLY: anneeref
    5462  USE time_phylmdz_mod,   ONLY: annee_ref
    5563  USE ioipsl,             ONLY: ioget_year_len, ioget_calendar
     
    5866  USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
    5967  USE slopes_m,           ONLY: slopes
     68  USE xios
     69  USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
     70  USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
     71  USE mod_grid_phy_lmdz, ONLY: klon_glo
     72
    6073!-------------------------------------------------------------------------------
    6174! Arguments:
     
    8396  CHARACTER(LEN=20) :: cal_in              ! Calendar
    8497  REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:)   ! Ozone climatologies
     98  REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:)   ! Ozone climatologies
    8599  REAL, ALLOCATABLE :: o3_in2  (:,:,:,:)   ! Ozone climatologies
     100  REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:)   ! Ozone climatologies
    86101  ! last index: 1 for the day-night average, 2 for the daylight field.
    87102  REAL :: NaN
     
    91106  REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13   ,:)
    92107  REAL, ALLOCATABLE :: o3_out3       (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
     108  REAL, ALLOCATABLE :: o3_out3_glo   (:,:,:,:) !   (nbp_lat,:,ntim_ou,:)
    93109  REAL, ALLOCATABLE :: o3_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
    94110  REAL, ALLOCATABLE :: o3_out2         (:,:,:,:) !         (nlat_ou,:,ntim_ou,:)
     111  REAL, ALLOCATABLE :: o3_out2_glo     (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
     112  REAL, ALLOCATABLE :: o3_out          (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
    95113! Dimension number  | Interval                | Contains  | For variables:
    96114!   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
     
    116134  INTEGER, ALLOCATABLE :: sta(:), cnt(:)
    117135  CHARACTER(LEN=80) :: sub, dim_nam, msg
    118 !-------------------------------------------------------------------------------
    119   sub="regr_horiz_time_climoz"
    120   WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
    121   CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
    122 
    123   CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
    124   lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
    125   lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
    126 
    127   !--- Get coordinates from the input file. Converts lon/lat in radians.
    128   !    Few inversions because "regr_conserv" and gcm need ascending vectors.
    129   CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
    130   CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
    131   l3D=ndims==4; l2D=ndims==3
    132   IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
    133   IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
    134   DO i=1,ndims
    135     CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
    136     CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
    137     ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE
    138     SELECT CASE(ii)
    139       CASE(1)                                         !--- LONGITUDE
    140         CALL NF95_GW_VAR(fID_in, varid, lon_in)
    141         ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
    142         nlon_in=dln; lon_in=lon_in*deg2rad
    143       CASE(2)                                         !--- LATITUDE
    144         CALL NF95_GW_VAR(fID_in, varid, lat_in)
    145         ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
    146         nlat_in=dln; lat_in=lat_in*deg2rad
    147       CASE(3)                                         !--- PRESSURE LEVELS
    148         CALL NF95_GW_VAR(fID_in, varid, lev_in)
    149         ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
    150         nlev_in=dln
    151         CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
    152         k=LEN_TRIM(press_unit)
    153         DO WHILE(ICHAR(press_unit(k:k))==0)
    154           press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
    155         END DO
    156         IF(press_unit ==  "Pa") THEN
    157           lev_in = lev_in/100.                        !--- CONVERT TO hPa
    158         ELSE IF(press_unit /= "hPa") THEN
    159           CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
     136  REAL :: null_array(0)
     137  LOGICAL,SAVE :: first=.TRUE.
     138!$OMP THREADPRIVATE(first) 
     139  REAL, ALLOCATABLE :: test_o3_in(:,:)
     140  REAL, ALLOCATABLE :: test_o3_out(:)
     141
     142
     143  IF (grid_type==unstructured) THEN
     144    IF (first) THEN
     145      IF (is_master) THEN
     146        ALLOCATE(latitude_glo(klon_glo))
     147        ALLOCATE(ind_cell_glo_glo(klon_glo))
     148      ELSE
     149        ALLOCATE(latitude_glo(0))
     150        ALLOCATE(ind_cell_glo_glo(0))
     151      ENDIF
     152      CALL gather(latitude_deg,  latitude_glo)
     153      CALL gather(ind_cell_glo,  ind_cell_glo_glo)
     154    ENDIF
     155  ENDIF
     156   
     157  IF (is_omp_master) THEN
     158    nlat_ou=nbp_lat
     159    nlon_ou=nbp_lon
     160   
     161   !-------------------------------------------------------------------------------
     162    IF (is_mpi_root) THEN
     163      sub="regr_horiz_time_climoz"
     164      WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
     165      CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
     166
     167      CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
     168      lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
     169      lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
     170
     171      !--- Get coordinates from the input file. Converts lon/lat in radians.
     172      !    Few inversions because "regr_conserv" and gcm need ascending vectors.
     173      CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
     174      CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
     175      l3D=ndims==4; l2D=ndims==3
     176      IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
     177      IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
     178      DO i=1,ndims
     179        CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
     180        CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
     181        ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE
     182        SELECT CASE(ii)
     183          CASE(1)                                         !--- LONGITUDE
     184            CALL NF95_GW_VAR(fID_in, varid, lon_in)
     185            ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
     186            nlon_in=dln; lon_in=lon_in*deg2rad
     187          CASE(2)                                         !--- LATITUDE
     188            CALL NF95_GW_VAR(fID_in, varid, lat_in)
     189            ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
     190            nlat_in=dln; lat_in=lat_in*deg2rad
     191          CASE(3)                                         !--- PRESSURE LEVELS
     192            CALL NF95_GW_VAR(fID_in, varid, lev_in)
     193            ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
     194            nlev_in=dln
     195            CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
     196            k=LEN_TRIM(press_unit)
     197            DO WHILE(ICHAR(press_unit(k:k))==0)
     198              press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
     199            END DO
     200            IF(press_unit ==  "Pa") THEN
     201              lev_in = lev_in/100.                        !--- CONVERT TO hPa
     202            ELSE IF(press_unit /= "hPa") THEN
     203              CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
     204            END IF
     205          CASE(4)                                         !--- TIME
     206            CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
     207            cal_in='gregorian'
     208            IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        & 
     209            WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
     210            TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
     211            k=LEN_TRIM(cal_in)
     212            DO WHILE(ICHAR(cal_in(k:k))==0)
     213              cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
     214            END DO
     215        END SELECT
     216      END DO
     217
     218      !--- Prepare quantities for time interpolation
     219      tmidmonth=mid_month(annee_ref, cal_in)
     220      IF(interpt) THEN
     221        ntim_ou=ioget_year_len(annee_ref)
     222        ALLOCATE(tmidday(ntim_ou))
     223        tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
     224        CALL ioget_calendar(cal_ou)
     225      ELSE
     226        ntim_ou=14
     227        cal_ou=cal_in
     228      END IF
     229    ENDIF
     230
     231    IF (grid_type==unstructured) THEN
     232      CALL bcast_mpi(nlon_in)
     233      CALL bcast_mpi(nlat_in)
     234      CALL bcast_mpi(nlev_in)
     235      CALL bcast_mpi(l3d)
     236      CALL bcast_mpi(tmidmonth)
     237      CALL bcast_mpi(tmidday)
     238      CALL bcast_mpi(ntim_ou)
     239   
     240      IF (is_mpi_root) THEN
     241        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
     242        IF (l3D) THEN
     243          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
     244        ELSE
     245          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
     246        ENDIF
     247      ELSE
     248        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
     249        IF (l3D) THEN
     250          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
     251        ELSE
     252          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
     253        ENDIF
     254      ENDIF
     255      CALL  xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
     256      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
     257      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
     258      CALL  xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
     259      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
     260      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
     261
     262      IF (first) THEN
     263        first=.FALSE.
     264        RETURN
     265      ENDIF
     266    ENDIF
     267   
     268   
     269    IF (is_mpi_root) THEN     
     270      !--- Longitudes management:
     271      !    * Need to shift data if the origin of input file longitudes /= -pi
     272      !    * Need to add some margin in longitude to ensure input interval contains
     273      !      all the output intervals => at least one longitudes slice has to be
     274      !      duplicated, possibly more for undersampling.
     275      IF(l3D) THEN
     276        IF (grid_type==unstructured) THEN
     277          dx2=0
     278        ELSE
     279          !--- Compute input edges longitudes vector (no end point yet)
     280          ALLOCATE(v1(nlon_in+1))
     281          v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
     282          FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
     283          v1(nlon_in+1)=v1(1)+2.*pi
     284          DEALLOCATE(lon_in)
     285
     286          !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
     287          v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
     288
     289          !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
     290          dx1=locate(v1,boundslon_reg(1,west))-1
     291          v1=CSHIFT(v1,SHIFT=dx1,DIM=1)
     292          v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
     293   
     294          !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
     295          dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
     296
     297          !--- Final edges longitudes vector (with margin and end point)
     298          ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
     299          DEALLOCATE(v1)
     300        ENDIF
     301      END IF
     302
     303      !--- Compute sinus of intervals edges latitudes:
     304      ALLOCATE(sinlat_in_edge(nlat_in+1))
     305      sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
     306      FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
     307      DEALLOCATE(lat_in)
     308
     309
     310
     311      !--- Check for contiguous years:
     312      ib=0; ie=13
     313      IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
     314        WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
     315      ELSE 
     316        IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
     317        IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
     318        IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
     319        IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
     320        IF(.NOT.lprev) ib=1
     321        IF(.NOT.lnext) ie=12
     322      END IF
     323      ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1 
     324      IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
     325      IF(l2D) cnt=[        nlat_in,nlev_in,1] 
     326      IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
     327      IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
     328
     329      !--- Read full current file and one record each available contiguous file
     330      DO iv=1,read_climoz
     331        msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
     332        CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
     333        IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
     334        IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
     335        CALL handle_err(TRIM(msg), ncerr, fID_in)
     336        IF(lprev) THEN; sta(ndims)=12 
     337          CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
     338          IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
     339          IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
     340          CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
    160341        END IF
    161       CASE(4)                                         !--- TIME
    162         CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
    163         cal_in='gregorian'
    164         IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        &
    165           WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
    166           TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
    167         k=LEN_TRIM(cal_in)
    168         DO WHILE(ICHAR(cal_in(k:k))==0)
    169           cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
    170         END DO
    171     END SELECT
    172   END DO
    173 
    174   !--- Longitudes management:
    175   !    * Need to shift data if the origin of input file longitudes /= -pi
    176   !    * Need to add some margin in longitude to ensure input interval contains
    177   !      all the output intervals => at least one longitudes slice has to be
    178   !      duplicated, possibly more for undersampling.
    179   IF(l3D) THEN
    180     !--- Compute input edges longitudes vector (no end point yet)
    181     ALLOCATE(v1(nlon_in+1))
    182     v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
    183     FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
    184     v1(nlon_in+1)=v1(1)+2.*pi
    185     DEALLOCATE(lon_in)
    186 
    187     !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
    188     v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
    189 
    190     !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
    191     dx1=locate(v1,boundslon_reg(1,west))-1
    192     v1=CSHIFT(v1,SHIFT=dx1,DIM=1); v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
    193 
    194     !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
    195     dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
    196 
    197     !--- Final edges longitudes vector (with margin and end point)
    198     ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
    199     DEALLOCATE(v1)
    200   END IF
    201 
    202   !--- Compute sinus of intervals edges latitudes:
    203   ALLOCATE(sinlat_in_edge(nlat_in+1))
    204   sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
    205   FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
    206   DEALLOCATE(lat_in)
    207 
    208   !--- Prepare quantities for time interpolation
    209   tmidmonth=mid_month(annee_ref, cal_in)
    210   IF(interpt) THEN
    211     ntim_ou=ioget_year_len(annee_ref)
    212     ALLOCATE(tmidday(ntim_ou))
    213     tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
    214     CALL ioget_calendar(cal_ou)
    215   ELSE
    216     ntim_ou=14
    217     cal_ou=cal_in
    218   END IF
    219 
    220   !--- Create the output file and get the variable IDs:
    221   CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
    222                    ndims, cal_ou)
    223 
    224   !--- Write remaining coordinate variables:
    225   CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
    226   IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
    227   IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
    228 
    229   !--- Check for contiguous years:
    230   ib=0; ie=13
    231   IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
    232     WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
    233   ELSE
    234     IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
    235     IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
    236     IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
    237     IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
    238     IF(.NOT.lprev) ib=1
    239     IF(.NOT.lnext) ie=12
    240   END IF
    241   ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1
    242   IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
    243   IF(l2D) cnt=[        nlat_in,nlev_in,1]
    244   IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
    245   IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
    246 
    247   !--- Read full current file and one record each available contiguous file
    248   DO iv=1,read_climoz
    249     msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
    250     CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
    251     IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
    252     IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
    253     CALL handle_err(TRIM(msg), ncerr, fID_in)
    254     IF(lprev) THEN; sta(ndims)=12
    255       CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
    256       IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
    257       IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
    258       CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
     342        IF(lnext) THEN; sta(ndims)=1 
     343          CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
     344          IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
     345          IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
     346          CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
     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       nlat_ou=klon_mpi
     418       
     419       IF (is_mpi_root) THEN
     420         ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz))
     421         o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
     422       ELSE
     423         ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
     424       ENDIF
     425       ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
     426       
     427       CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
     428       CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
     429
     430     ELSE
     431         
     432       !--- Regrid in longitude
     433        ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
     434        CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
     435                            xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
     436                            vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
     437        DEALLOCATE(o3_in3)
     438
     439        !--- Regrid in latitude: averaging with respect to SIN(lat) is
     440        !                        equivalent to weighting by COS(lat)
     441        !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
     442        ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
     443        CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
     444                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
     445                        vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
     446                   slope = slopes(2,o3_regr_lon, sinlat_in_edge))
     447        DEALLOCATE(o3_regr_lon)
     448
     449     ENDIF
     450
     451     !--- Duplicate previous/next record(s) if they are not available
     452     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
     453     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
     454     
     455     !--- Regrid in time by linear interpolation:
     456     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
     457     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
     458     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
     459     DEALLOCATE(o3_regr_lonlat)
     460
     461     nlat_ou=nbp_lat
     462     IF (grid_type==unstructured) THEN
     463       CALL xios_send_field('o3_out',o3_out3)
     464       ndims=3
     465       ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
     466       CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
     467     ENDIF
     468
     469    !--- Create the output file and get the variable IDs:
     470    CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
     471                     ndims, cal_ou)
     472
     473    IF (is_mpi_root) THEN
     474      !--- Write remaining coordinate variables:
     475      CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
     476      IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
     477      IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
     478
     479      !--- Write to file (the order of "rlatu" is inverted in the output file):
     480        IF (grid_type==unstructured) THEN
     481
     482          ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
     483          DO i=1,klon_glo
     484            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
     485          ENDDO
     486
     487          DO m = 1, read_climoz
     488            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
     489          END DO
     490         
     491        ELSE
     492          DO m = 1, read_climoz
     493            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
     494          END DO
     495      ENDIF
     496      CALL NF95_CLOSE(fID_ou)
     497
     498
     499    ENDIF
     500
     501
     502    !=============================================================================
     503    ELSE                                                         !=== ZONAL FIELDS
     504    !=============================================================================
     505   
     506     IF (grid_type==unstructured) THEN
     507       nlat_ou=klon_mpi
     508
     509       IF (is_mpi_root) THEN
     510         ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz))
     511         o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
     512       ELSE
     513         ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
     514       ENDIF
     515       ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
     516       CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
     517       CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
     518       IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
     519       IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
     520       
     521     
     522     ELSE
     523        !--- Regrid in latitude: averaging with respect to SIN(lat) is
     524        !                        equivalent to weighting by COS(lat)
     525        !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
     526        ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
     527        CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
     528                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
     529                        vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
     530                     slope = slopes(1,o3_in2, sinlat_in_edge))
     531        DEALLOCATE(o3_in2)
     532
     533        !--- Duplicate previous/next record(s) if they are not available
     534        IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
     535        IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
     536
     537     ENDIF
     538     
     539      !--- Regrid in time by linear interpolation:
     540      ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
     541      IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
     542      IF(.NOT.interpt) o3_out2=o3_regr_lat
     543      DEALLOCATE(o3_regr_lat)
     544
     545      nlat_ou=nbp_lat
     546   
     547      IF (grid_type==unstructured) THEN
     548        ndims=3
     549        ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
     550        CALL gather_mpi(o3_out2, o3_out2_glo)
     551      ENDIF
     552     
     553      !--- Create the output file and get the variable IDs:
     554      CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
     555                         ndims, cal_ou)
     556
     557      IF (is_mpi_root) THEN
     558     
     559        !--- Write remaining coordinate variables:
     560        CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
     561        IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
     562        IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
     563
     564        IF (grid_type==unstructured) THEN
     565
     566          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
     567          DO i=1,klon_glo
     568            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
     569          ENDDO
     570
     571
     572          DO m = 1, read_climoz
     573            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
     574          END DO
     575        ELSE
     576          !--- Write to file (the order of "rlatu" is inverted in the output file):
     577          DO m = 1, read_climoz
     578            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
     579          END DO
     580        ENDIF
     581       
     582        CALL NF95_CLOSE(fID_ou)
     583     
     584      ENDIF
     585
     586    !=============================================================================
    259587    END IF
    260     IF(lnext) THEN; sta(ndims)=1
    261       CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
    262       IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
    263       IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
    264       CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
    265     END IF
    266   END DO
    267   IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
    268   IF(lprev) CALL NF95_CLOSE(fID_in_m)
    269   IF(lnext) CALL NF95_CLOSE(fID_in_p)
    270 
    271   !--- Revert decreasing coordinates vector
    272   IF(l3D) THEN
    273     IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
    274     IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
    275     IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
    276     !--- Shift values for longitude and duplicate some longitudes slices
    277     o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
    278     o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
    279   ELSE
    280     IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
    281     IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
    282   END IF
    283 
    284  !--- Deal with missing values
    285   DO m=1, read_climoz
    286     WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
    287     IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
    288       IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
    289         WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
    290       END IF
    291     END IF
    292     WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
    293     WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
    294 
    295     !--- Check top layer contains no NaNs & search NaNs from top to ground
    296     msg=TRIM(sub)//": NaNs in top layer !"
    297     IF(l3D) THEN
    298       IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
    299       DO k = 2,nlev_in
    300         WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
    301       END DO
    302     ELSE
    303       IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
    304         WRITE(lunout,*)msg
    305         !--- Fill in latitudes where all values are missing
    306         DO l=1,nmth_in
    307           !--- Next to south pole
    308           j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
    309           IF(j>1) &
    310             o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
    311           !--- Next to north pole
    312           j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
    313           IF(j<nlat_in) &
    314             o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
    315         END DO
    316       END IF
    317 
    318       !--- Fill in high latitudes missing values
    319       !--- Highest level been filled-in, so has always valid values.
    320       DO k = 2,nlev_in
    321         WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
    322       END DO
    323     END IF
    324   END DO
    325   CALL NF95_CLOSE(fID_in)
    326 
    327   !=============================================================================
    328   IF(l3D) THEN                                                   !=== 3D FIELDS
    329   !=============================================================================
    330     !--- Regrid in longitude
    331     ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
    332     CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
    333                         xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
    334                         vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
    335     DEALLOCATE(o3_in3)
    336 
    337     !--- Regrid in latitude: averaging with respect to SIN(lat) is
    338     !                        equivalent to weighting by COS(lat)
    339     !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
    340     ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
    341     CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
    342                     xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
    343                     vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
    344                  slope = slopes(2,o3_regr_lon, sinlat_in_edge))
    345     DEALLOCATE(o3_regr_lon)
    346 
    347     !--- Duplicate previous/next record(s) if they are not available
    348     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
    349     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
    350 
    351     !--- Regrid in time by linear interpolation:
    352     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
    353     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
    354     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
    355     DEALLOCATE(o3_regr_lonlat)
    356 
    357     !--- Write to file (the order of "rlatu" is inverted in the output file):
    358     DO m = 1, read_climoz
    359       CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
    360     END DO
    361 
    362   !=============================================================================
    363   ELSE                                                         !=== ZONAL FIELDS
    364   !=============================================================================
    365     !--- Regrid in latitude: averaging with respect to SIN(lat) is
    366     !                        equivalent to weighting by COS(lat)
    367     !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
    368     ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
    369     CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
    370                     xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
    371                     vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
    372                  slope = slopes(1,o3_in2, sinlat_in_edge))
    373     DEALLOCATE(o3_in2)
    374 
    375     !--- Duplicate previous/next record(s) if they are not available
    376     IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
    377     IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
    378 
    379     !--- Regrid in time by linear interpolation:
    380     ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
    381     IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
    382     IF(.NOT.interpt) o3_out2=o3_regr_lat
    383     DEALLOCATE(o3_regr_lat)
    384 
    385     !--- Write to file (the order of "rlatu" is inverted in the output file):
    386     DO m = 1, read_climoz
    387       CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
    388     END DO
    389 
    390   !=============================================================================
    391   END IF
    392   !=============================================================================
    393 
    394   CALL NF95_CLOSE(fID_ou)
    395 
     588    !=============================================================================
     589
     590    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
     591
     592  ENDIF ! is_omp_master
     593
     594  first=.FALSE.
    396595END SUBROUTINE regr_horiz_time_climoz
    397596!
     
    408607!-------------------------------------------------------------------------------
    409608  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
     609  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
     610  USE mod_phys_lmdz_para, ONLY: is_mpi_root
     611  USE mod_grid_phy_lmdz, ONLY: klon_glo
     612!
    410613!-------------------------------------------------------------------------------
    411614! Arguments:
     
    419622  INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
    420623  INTEGER :: vlonID, vlatID, ncerr,  is
     624  REAL,ALLOCATABLE    :: latitude_glo_(:)
    421625  CHARACTER(LEN=80) :: sub
    422 !-------------------------------------------------------------------------------
    423   sub="prepare_out"
    424   WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
    425   CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
     626  INTEGER :: i
     627
     628
     629!-------------------------------------------------------------------------------
     630 
     631  IF (is_mpi_root) THEN 
     632    sub="prepare_out"
     633    WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
     634    CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
    426635
    427636  !--- Dimensions:
    428   IF(ndims==4) &
    429   CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
    430   CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
    431   CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
    432   CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
    433 
    434   !--- Define coordinate variables:
    435   IF(ndims==4) &
    436   CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
    437   CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
    438   CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
    439   CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
    440   IF(ndims==4) &
    441   CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
    442   CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
    443   CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
    444   CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
    445   IF(ndims==4) &
    446   CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
    447   CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
    448   CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
    449   CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
    450   CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
    451   CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
     637    IF(ndims==4) &
     638    CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
     639    CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
     640    CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
     641    CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
     642
     643    !--- Define coordinate variables:
     644    IF(ndims==4) &
     645    CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
     646    CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
     647    CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
     648    CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
     649    IF(ndims==4) &
     650    CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
     651    CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
     652    CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
     653    CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
     654    IF(ndims==4) &
     655    CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
     656    CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
     657    CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
     658    CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
     659    CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
     660    CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
    452661
    453662  !--- Define the main variables:
    454   IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
    455   IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
    456   CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
    457   CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
    458   CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
     663    IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
     664    IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
     665    CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
     666    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
     667    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
    459668      &_in_air")
    460   IF(SIZE(vID_ou) == 2) THEN
    461     CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
    462     CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
    463       &ylight")
    464   END IF
     669    IF(SIZE(vID_ou) == 2) THEN
     670      CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
     671      CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
     672        &ylight")
     673    END IF
    465674
    466675  !--- Global attributes:
    467676  ! The following commands, copying attributes, may fail. That is OK.
    468677  ! It should just mean that the attribute is not defined in the input file.
    469   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
    470   CALL handle_err_copy_att("Conventions")
    471   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
    472   CALL handle_err_copy_att("title")
    473   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
    474   CALL handle_err_copy_att("institution")
    475   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
    476   CALL handle_err_copy_att("source")
    477   CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
    478   CALL NF95_ENDDEF(fID_ou)
    479 
    480   !--- Write one of the coordinate variables:
    481   IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
    482   CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
    483   !    (convert from rad to degrees and sort in ascending order)
    484 
     678    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
     679    CALL handle_err_copy_att("Conventions")
     680    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
     681    CALL handle_err_copy_att("title")
     682    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
     683    CALL handle_err_copy_att("institution")
     684    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
     685    CALL handle_err_copy_att("source")
     686    CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
     687    CALL NF95_ENDDEF(fID_ou)
     688
     689    IF (grid_type==unstructured) THEN
     690      ALLOCATE(latitude_glo_(klon_glo))
     691      DO i=1,klon_glo
     692        latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
     693      ENDDO
     694      CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
     695    ELSE
     696      !--- Write one of the coordinate variables:
     697      IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
     698      CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
     699    !    (convert from rad to degrees and sort in ascending order)
     700    ENDIF
     701  ENDIF
     702 
    485703CONTAINS
    486704
Note: See TracChangeset for help on using the changeset viewer.