Ignore:
Timestamp:
Nov 21, 2019, 4:43:45 PM (4 years ago)
Author:
lguez
Message:

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/regr_horiz_time_climoz_m.F90

    r3278 r3605  
    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
    54   USE control_mod,        ONLY: anneeref
     61!!  USE control_mod,        ONLY: anneeref
     62  USE time_phylmdz_mod,   ONLY: annee_ref
    5563  USE ioipsl,             ONLY: ioget_year_len, ioget_calendar
    5664  USE regr_conserv_m,     ONLY: regr_conserv
     
    5866  USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
    5967  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
    6075!-------------------------------------------------------------------------------
    6176! Arguments:
     
    8398  CHARACTER(LEN=20) :: cal_in              ! Calendar
    8499  REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:)   ! Ozone climatologies
     100  REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:)   ! Ozone climatologies
    85101  REAL, ALLOCATABLE :: o3_in2  (:,:,:,:)   ! Ozone climatologies
     102  REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:)   ! Ozone climatologies
    86103  ! last index: 1 for the day-night average, 2 for the daylight field.
    87104  REAL :: NaN
     
    91108  REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13   ,:)
    92109  REAL, ALLOCATABLE :: o3_out3       (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
     110  REAL, ALLOCATABLE :: o3_out3_glo   (:,:,:,:) !   (nbp_lat,:,ntim_ou,:)
    93111  REAL, ALLOCATABLE :: o3_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
    94112  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,:)
    95115! Dimension number  | Interval                | Contains  | For variables:
    96116!   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
     
    116136  INTEGER, ALLOCATABLE :: sta(:), cnt(:)
    117137  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)
     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)
    160345        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+2:)=v1(nlon_in-dx1+2:)+2.*pi
    193 
    194     !--- Extend input longitudes vector until last interval contains boundslon_reg(nlon_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(anneeref, cal_in)
    210   IF(interpt) THEN
    211     ntim_ou=ioget_year_len(anneeref)
    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)
     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    !=============================================================================
    259595    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 
     596    !=============================================================================
     597
     598    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
     599
     600  ENDIF ! is_omp_master
     601
     602  first=.FALSE.
    396603END SUBROUTINE regr_horiz_time_climoz
    397604!
     
    408615!-------------------------------------------------------------------------------
    409616  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!
    410621!-------------------------------------------------------------------------------
    411622! Arguments:
     
    419630  INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
    420631  INTEGER :: vlonID, vlatID, ncerr,  is
     632  REAL,ALLOCATABLE    :: latitude_glo_(:)
    421633  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)
     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)
    426643
    427644  !--- 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)
     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)
    452669
    453670  !--- 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&
     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&
    459676      &_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
     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
    465682
    466683  !--- Global attributes:
    467684  ! The following commands, copying attributes, may fail. That is OK.
    468685  ! 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 
     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 
    485711CONTAINS
    486712
Note: See TracChangeset for help on using the changeset viewer.