Index: LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 =================================================================== --- LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (revision 2539) +++ LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (revision 2540) @@ -60,4 +60,5 @@ ! * 11/2009: L. Guez (ozone day & night climatos, see etat0_netcdf.F90) ! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs) +! * 04/2016: D. Cugnet (12/14 recs SST/SIC files: cyclic/interannual runs) !------------------------------------------------------------------------------- #ifndef CPP_1D @@ -339,7 +340,7 @@ INTEGER :: ndays_in ! number of days !--- misc - INTEGER :: i, j, k, l ! loop counters + INTEGER :: i, j, k, l, ll ! loop counters REAL, ALLOCATABLE :: work(:,:) ! used for extrapolation - CHARACTER(LEN=25) :: title ! for messages + CHARACTER(LEN=128):: title, mess ! for messages LOGICAL :: extrp ! flag for extrapolation REAL :: chmin, chmax @@ -392,5 +393,5 @@ !--- Time (variable is not needed - it is rebuilt - but calendar is) CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam) - ALLOCATE(timeyear(lmdep)) + ALLOCATE(timeyear(MAX(lmdep,14))) CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam) cal_in=' ' @@ -405,16 +406,25 @@ CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep) -!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- - !--- Determining input file number of days, depending on calendar +!--- Determining input file number of days, depending on calendar ndays_in=year_len(anneeref, cal_in) -!--- Time vector reconstruction (time vector from file is not trusted) -!--- If input records are not monthly, time sampling has to be constant ! - timeyear=mid_months(anneeref, cal_in, lmdep) - IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), & - ' ne comportent pas 12, mais ', lmdep, ' enregistrements.' +!--- Rebuilding input time vector (field from input file might be unreliable) + IF(lmdep==12) THEN + timeyear=mid_month(anneeref, cal_in, ndays_in) + CALL msg(0,'Monthly periodic input file (perpetual run).') + ELSE IF(lmdep==14) THEN + timeyear=mid_month(anneeref, cal_in, ndays_in) + CALL msg(0,'Monthly 14-records input file (interannual run).') + ELSE IF(lmdep==ndays_in) THEN + timeyear=[(REAL(k)-0.5,k=1,ndays_in)] + CALL msg(0,'Daily input file (no time interpolation).') + ELSE + WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep, & + ' records, 12/14/',ndays_in,' (periodic/interannual/daily) needed' + CALL abort_physic('mid_months',TRIM(mess),1) + END IF !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- - ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep)) + ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14))) IF(extrp) ALLOCATE(work(imdep, jmdep)) CALL msg(5,'') @@ -436,6 +446,12 @@ WHERE(NINT(mask)/=1) champint=0.001 END IF - champtime(:, :, l)=champint + !--- Special case for periodic input file: index shifted + ll = l; IF(lmdep==12) ll = l + 1 + champtime(:, :, ll)=champint END DO + IF(lmdep==12) THEN + champtime(:,:, 1)=champtime(:,:,13) + champtime(:,:,14)=champtime(:,:, 2) + END IF CALL ncerr(NF90_CLOSE(ncid), fnam) @@ -444,27 +460,39 @@ !--- TIME INTERPOLATION ------------------------------------------------------ - IF(prt_level>5) THEN - WRITE(lunout, *) - WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' - WRITE(lunout, *)' Vecteur temps en entree: ', timeyear - WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays - END IF - - ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) - skip = .false. - n_extrap = 0 - DO j=1, jjp1 - DO i=1, iim - yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, & - vc_beg=0., vc_end=0.) - CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, & - arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr) - if (ierr < 0) stop 1 - n_extrap = n_extrap + ierr - END DO - END DO - if (n_extrap /= 0) then - WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap - end if + IF(prt_level>0) THEN + IF(ndays/=ndays_in) THEN + WRITE(lunout, *) + WRITE(lunout,*)'DIFFERENTES LONGEURS D ANNEES:' + WRITE(lunout,*)' Dans le fichier d entree: ',ndays_in + WRITE(lunout,*)' Dans le fichier de sortie: ',ndays + END IF + IF(lmdep==ndays_in) THEN + WRITE(lunout, *) + WRITE(lunout, *)'PAS D INTERPOLATION TEMPORELLE.' + WRITE(lunout, *)' Fichier journalier en entree.' + ELSE + WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' + WRITE(lunout, *)' Vecteur temps en entree: ', timeyear + WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays-1 + END IF + END IF + ALLOCATE(yder(MAX(lmdep,14)), champan(iip1, jjp1, ndays)) + IF(lmdep==ndays_in) THEN + champan(1:iim,:,:)=champtime + ELSE + skip = .false. + n_extrap = 0 + DO j=1, jjp1 + DO i=1, iim + yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, & + vc_beg=0., vc_end=0.) + CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, & + arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr) + if (ierr < 0) stop 1 + n_extrap = n_extrap + ierr + END DO + END DO + IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap + END IF champan(iip1, :, :)=champan(1, :, :) DEALLOCATE(yder, champtime, timeyear) @@ -752,5 +780,5 @@ !------------------------------------------------------------------------------- ! -FUNCTION mid_months(y,cal_in,nm) +FUNCTION mid_month(y,cal_in,ndays_out) ! !------------------------------------------------------------------------------- @@ -758,48 +786,35 @@ !------------------------------------------------------------------------------- ! Arguments: - INTEGER, INTENT(IN) :: y ! year - CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar - INTEGER, INTENT(IN) :: nm ! months/year number - REAL, DIMENSION(nm) :: mid_months ! mid-month times + INTEGER, INTENT(IN) :: y ! year + CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar + INTEGER, INTENT(IN) :: ndays_out ! days number + REAL, DIMENSION(14) :: mid_month ! mid-bins times !------------------------------------------------------------------------------- ! Local variables: - CHARACTER(LEN=99) :: mess ! error message - CHARACTER(LEN=20) :: cal_out ! calendar (for outputs) - INTEGER, DIMENSION(nm) :: mnth ! months lengths (days) - INTEGER :: m ! months counter - INTEGER :: nd ! number of days -!------------------------------------------------------------------------------- - nd=year_len(y,cal_in) - - IF(nm==12) THEN - - !--- Getting the input calendar to reset at the end of the function - CALL ioget_calendar(cal_out) - - !--- Unlocking calendar and setting it to wanted one - CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in)) - - !--- Getting the length of each month - DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO + CHARACTER(LEN=99) :: mess ! error message + CHARACTER(LEN=20) :: cal_out ! output calendar + INTEGER, DIMENSION(14) :: tlen ! months lengths (days) + INTEGER :: m ! months counter + INTEGER :: nd ! number of days +!------------------------------------------------------------------------------- + !--- Save the input calendar to reset it at the end of the function + CALL ioget_calendar(cal_out) + + !--- Unlock calendar and set it to "cal_in" + CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in)) + + !--- Get the length of each month + tlen(1 )=ioget_mon_len(y-1,12) + DO m=1,12; tlen(m+1)=ioget_mon_len(y,m); END DO + tlen(14)=ioget_mon_len(y+1, 1) + + !--- Mid-bins times + mid_month(1)=-0.5*REAL(tlen(1)) + DO m=2,14; mid_month(m)=mid_month(m-1)+0.5*REAL(tlen(m-1)+tlen(m)); END DO !--- Back to original calendar - CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out)) - - ELSE IF(MODULO(nd,nm)/=0) THEN - WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',& - nm,' months/year. Months number should divide days number.' - CALL abort_physic('mid_months',TRIM(mess),1) - - ELSE - mnth=[(m,m=1,nm,nd/nm)] - END IF - -!--- Mid-months times - mid_months(1)=0.5*REAL(mnth(1)) - DO k=2,nm - mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k)) - END DO - -END FUNCTION mid_months + CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out)) + +END FUNCTION mid_month ! !------------------------------------------------------------------------------- @@ -863,2 +878,3 @@ ! !******************************************************************************* +