Changeset 2540


Ignore:
Timestamp:
Jun 3, 2016, 8:02:41 PM (8 years ago)
Author:
dcugnet
Message:

Improved time interpolation for "limit.nc" files, depending on the number
of records found in the input SST/SIC files:

  • 14 records => for interannual runs. Record nr. 1 contains the december field of the previous year. Record nr. 14 contains the january field of the next year.
  • 12 records => for climatic runs. Records 1 and 14 are duplicated from records 13 and 2.
  • as much records as need by the LMDZ calendar for the considered year => no time interpolation at all (typically: for guided runs with SST and SIC known every day).

In the first case: typical SST daily increment is about 0.04K at 3 sigmas.
The gap between 2 years (last day of previous - first day of current year)
decreases from 1K at 3 sigmas (former method) to 0.25K (improved method).
Could be improved a bit by taking even more records (16 ?) and a lot by
interpolating the whole time serie together.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90

    r2399 r2540  
    6060!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
    6161!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
     62!  *    04/2016: D. Cugnet   (12/14 recs SST/SIC files: cyclic/interannual runs)
    6263!-------------------------------------------------------------------------------
    6364#ifndef CPP_1D
    … …  
    339340  INTEGER           :: ndays_in           ! number of days
    340341!--- misc
    341   INTEGER           :: i, j, k, l         ! loop counters
     342  INTEGER           :: i, j, k, l, ll     ! loop counters
    342343  REAL, ALLOCATABLE :: work(:,:)          ! used for extrapolation
    343   CHARACTER(LEN=25) :: title              ! for messages
     344  CHARACTER(LEN=128):: title, mess        ! for messages
    344345  LOGICAL           :: extrp              ! flag for extrapolation
    345346  REAL              :: chmin, chmax
    … …  
    392393!--- Time (variable is not needed - it is rebuilt - but calendar is)
    393394  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
    394   ALLOCATE(timeyear(lmdep))
     395  ALLOCATE(timeyear(MAX(lmdep,14)))
    395396  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
    396397  cal_in=' '
    … …  
    405406  CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
    406407 
    407 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
    408   !--- Determining input file number of days, depending on calendar
     408!--- Determining input file number of days, depending on calendar
    409409  ndays_in=year_len(anneeref, cal_in)
    410410
    411 !--- Time vector reconstruction (time vector from file is not trusted)
    412 !--- If input records are not monthly, time sampling has to be constant !
    413   timeyear=mid_months(anneeref, cal_in, lmdep)
    414   IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
    415        ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
     411!--- Rebuilding input time vector (field from input file might be unreliable)
     412  IF(lmdep==12) THEN
     413    timeyear=mid_month(anneeref, cal_in, ndays_in)
     414    CALL msg(0,'Monthly periodic input file (perpetual run).')
     415  ELSE IF(lmdep==14) THEN
     416    timeyear=mid_month(anneeref, cal_in, ndays_in)
     417    CALL msg(0,'Monthly 14-records input file (interannual run).')
     418  ELSE IF(lmdep==ndays_in) THEN
     419    timeyear=[(REAL(k)-0.5,k=1,ndays_in)]
     420    CALL msg(0,'Daily input file (no time interpolation).')
     421  ELSE
     422    WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
     423      ' records, 12/14/',ndays_in,' (periodic/interannual/daily) needed'
     424    CALL abort_physic('mid_months',TRIM(mess),1)
     425  END IF
    416426
    417427!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
    418   ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
     428  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14)))
    419429  IF(extrp) ALLOCATE(work(imdep, jmdep))
    420430  CALL msg(5,'')
    … …  
    436446      WHERE(NINT(mask)/=1) champint=0.001
    437447    END IF
    438     champtime(:, :, l)=champint
     448    !--- Special case for periodic input file: index shifted
     449    ll = l; IF(lmdep==12) ll = l + 1
     450    champtime(:, :, ll)=champint
    439451  END DO
     452  IF(lmdep==12) THEN
     453    champtime(:,:, 1)=champtime(:,:,13)
     454    champtime(:,:,14)=champtime(:,:, 2)
     455  END IF
    440456  CALL ncerr(NF90_CLOSE(ncid), fnam)
    441457
    … …  
    444460
    445461!--- TIME INTERPOLATION ------------------------------------------------------
    446   IF(prt_level>5) THEN
    447      WRITE(lunout, *)
    448      WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
    449      WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
    450      WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
    451   END IF
    452 
    453   ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
    454   skip = .false.
    455   n_extrap = 0
    456   DO j=1, jjp1
    457     DO i=1, iim
    458       yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
    459            vc_beg=0., vc_end=0.)
    460       CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
    461            arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
    462       if (ierr < 0) stop 1
    463       n_extrap = n_extrap + ierr
    464     END DO
    465   END DO
    466   if (n_extrap /= 0) then
    467      WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
    468   end if
     462  IF(prt_level>0) THEN
     463     IF(ndays/=ndays_in) THEN
     464        WRITE(lunout, *)
     465        WRITE(lunout,*)'DIFFERENTES LONGEURS D ANNEES:'
     466        WRITE(lunout,*)' Dans le fichier d entree: ',ndays_in
     467        WRITE(lunout,*)' Dans le fichier de sortie: ',ndays
     468     END IF
     469     IF(lmdep==ndays_in) THEN
     470        WRITE(lunout, *)
     471        WRITE(lunout, *)'PAS D INTERPOLATION TEMPORELLE.'
     472        WRITE(lunout, *)' Fichier journalier en entree.'
     473     ELSE
     474        WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
     475        WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
     476        WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays-1
     477     END IF
     478  END IF
     479  ALLOCATE(yder(MAX(lmdep,14)), champan(iip1, jjp1, ndays))
     480  IF(lmdep==ndays_in) THEN
     481     champan(1:iim,:,:)=champtime
     482  ELSE
     483     skip = .false.
     484     n_extrap = 0
     485     DO j=1, jjp1
     486       DO i=1, iim
     487         yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
     488              vc_beg=0., vc_end=0.)
     489         CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
     490              arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
     491         if (ierr < 0) stop 1
     492         n_extrap = n_extrap + ierr
     493       END DO
     494     END DO
     495     IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
     496  END IF
    469497  champan(iip1, :, :)=champan(1, :, :)
    470498  DEALLOCATE(yder, champtime, timeyear)
    … …  
    752780!-------------------------------------------------------------------------------
    753781!
    754 FUNCTION mid_months(y,cal_in,nm)
     782FUNCTION mid_month(y,cal_in,ndays_out)
    755783!
    756784!-------------------------------------------------------------------------------
    … …  
    758786!-------------------------------------------------------------------------------
    759787! Arguments:
    760   INTEGER,                INTENT(IN) :: y               ! year
    761   CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
    762   INTEGER,                INTENT(IN) :: nm              ! months/year number
    763   REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
     788  INTEGER,          INTENT(IN) :: y           ! year
     789  CHARACTER(LEN=*), INTENT(IN) :: cal_in      ! calendar
     790  INTEGER,          INTENT(IN) :: ndays_out   ! days  number
     791  REAL,          DIMENSION(14) :: mid_month   ! mid-bins times
    764792!-------------------------------------------------------------------------------
    765793! Local variables:
    766   CHARACTER(LEN=99)                  :: mess            ! error message
    767   CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
    768   INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
    769   INTEGER                            :: m               ! months counter
    770   INTEGER                            :: nd              ! number of days
    771 !-------------------------------------------------------------------------------
    772   nd=year_len(y,cal_in)
    773 
    774   IF(nm==12) THEN
    775 
    776   !--- Getting the input calendar to reset at the end of the function
    777     CALL ioget_calendar(cal_out)
    778 
    779   !--- Unlocking calendar and setting it to wanted one
    780     CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
    781 
    782   !--- Getting the length of each month
    783     DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
     794  CHARACTER(LEN=99)      :: mess              ! error message
     795  CHARACTER(LEN=20)      :: cal_out           ! output calendar
     796  INTEGER, DIMENSION(14) :: tlen              ! months lengths (days)
     797  INTEGER                :: m                 ! months counter
     798  INTEGER                :: nd                ! number of days
     799!-------------------------------------------------------------------------------
     800  !--- Save the input calendar to reset it at the end of the function
     801  CALL ioget_calendar(cal_out)
     802
     803  !--- Unlock calendar and set it to "cal_in"
     804  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
     805
     806  !--- Get the length of each month
     807  tlen(1 )=ioget_mon_len(y-1,12)
     808  DO m=1,12; tlen(m+1)=ioget_mon_len(y,m); END DO
     809  tlen(14)=ioget_mon_len(y+1, 1)
     810
     811  !--- Mid-bins times
     812  mid_month(1)=-0.5*REAL(tlen(1))
     813  DO m=2,14; mid_month(m)=mid_month(m-1)+0.5*REAL(tlen(m-1)+tlen(m)); END DO
    784814
    785815  !--- Back to original calendar
    786     CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
    787 
    788   ELSE IF(MODULO(nd,nm)/=0) THEN
    789     WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
    790       nm,' months/year. Months number should divide days number.'
    791     CALL abort_physic('mid_months',TRIM(mess),1)
    792 
    793   ELSE
    794     mnth=[(m,m=1,nm,nd/nm)]
    795   END IF
    796 
    797 !--- Mid-months times
    798   mid_months(1)=0.5*REAL(mnth(1))
    799   DO k=2,nm
    800     mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
    801   END DO
    802 
    803 END FUNCTION mid_months
     816  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
     817
     818END FUNCTION mid_month
    804819!
    805820!-------------------------------------------------------------------------------
    … …  
    863878!
    864879!*******************************************************************************
     880
Note: See TracChangeset for help on using the changeset viewer.