Changeset 2766

Timestamp:

Jan 19, 2017, 4:56:16 PM (7 years ago)

Author:

dcugnet

Message:

• Fix limit_netcdf to handle correctly files with sea ice concentration in fractions and not percents. CMIP6 official files have "units" attribute equal to "1.0" with a null end character.
• Extraction of few calendar-related functions in cal_tools_m.F90.

Location:

LMDZ5/trunk/libf

Files:

• dynphy_lonlat/phylmd/limit_netcdf.F90 (modified) (10 diffs)
• phylmd/cal_tools_m.F90 (added)

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

@@ -17 +17 @@
 !-------------------------------------------------------------------------------
 
-  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo,          &
-   ioconf_calendar, ioget_calendar, lock_calendar, ioget_mon_len, ioget_year_len
+  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
   USE assert_eq_m,        ONLY: assert_eq
+  USE cal_tools_m,        ONLY: year_len, mid_month
   USE conf_dat_m,         ONLY: conf_dat2d, conf_dat3d
   USE dimphy,             ONLY: klon, zmasq
-  USE geometry_mod, ONLY: longitude_deg, latitude_deg
+  USE geometry_mod,       ONLY: longitude_deg, latitude_deg
   USE phys_state_var_mod, ONLY: pctsrf
-  USE control_mod, ONLY: anneeref
+  USE control_mod,        ONLY: anneeref
   USE init_ssrf_m,        ONLY: start_init_subsurf
 
@@ -123 +123 @@
 
 !--- Beware: anneeref (from gcm.def) is used to determine output time sampling
-  ndays=ioget_year_len(anneeref)
+  ndays=year_len(anneeref)
 
 !--- RUGOSITY TREATMENT --------------------------------------------------------
@@ -375 +375 @@
 !--- Read unit for sea-ice variable only
   IF (mode=='SIC') THEN
-    IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN
+    ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic); CALL strclean(unit_sic)
+    IF(ierr/=NF90_NOERR) THEN
       CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
       unit_sic='X'
+    ELSE IF(ANY(["%  ","1.0","1  "]==unit_sic)) THEN
+      CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
     ELSE
-      CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
+      CALL abort_physic('SIC','Unrecognized unit for sea-ice file: '&
+        //TRIM(unit_sic),1)
     END IF
   END IF
@@ -418 +422 @@
 !--- Rebuilding input time vector (field from input file might be unreliable)
   IF(lmdep==12) THEN
-    timeyear=mid_month(anneeref, cal_in, ndays_in)
+    timeyear=mid_month(anneeref, cal_in)
     CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
   ELSE IF(lmdep==ndays_in) THEN
@@ -426 +430 @@
     WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
       ' records, 12/',ndays_in,' (monthly/daily needed).'
-    CALL abort_physic('mid_months',TRIM(mess),1)
+    CALL abort_physic('mid_month',TRIM(mess),1)
   END IF
 
@@ -433 +437 @@
   IF(extrp) ALLOCATE(work(imdep, jmdep))
   CALL msg(5,'')
-  CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.')
+  CALL msg(5,'READ AND INTERPOLATE HORIZONTALLY ', lmdep, ' FIELDS.')
   CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
   DO l=1, lmdep
     CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
+    !--- Check whether values are acceptable for SIC, depending on unit.
+    IF(mode=='SIC') THEN
+      IF(TRIM(unit_sic)=="1".OR.TRIM(unit_sic)=="1.0") THEN
+        IF(ANY(champ>1.0+EPSFRA)) &
+          CALL abort_physic('SIC','Found sea-ice fractions greater than 1.')
+      ELSE IF(TRIM(unit_sic)=="X".OR.TRIM(unit_sic)=="%") THEN
+        IF(ANY(champ>100.0+EPSFRA)) &
+          CALL abort_physic('SIC','Found sea-ice percentages greater than 100.')
+      END IF
+    END IF
     CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
     IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     IF(l==1) THEN
       CALL msg(5,"----------------------------------------------------------")
-      CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
+      CALL msg(5,"$$$ Barycentrique interpolation for "//TRIM(title)//" $$$")
       CALL msg(5,"----------------------------------------------------------")
     END IF
@@ -504 +518 @@
      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
+        WRITE(lunout,*)'DIFFERENT YEAR LENGTHS:'
+        WRITE(lunout,*)' In the input file: ',ndays_in
+        WRITE(lunout,*)' In the output file: ',ndays
      END IF
      IF(lmdep==ndays_in) THEN
         WRITE(lunout, *)
-        WRITE(lunout, *)'PAS D INTERPOLATION TEMPORELLE.'
-        WRITE(lunout, *)' Fichier journalier en entree.'
+        WRITE(lunout, *)'NO TIME INTERPOLATION.'
+        WRITE(lunout, *)' Daily input file.'
      ELSE
-        WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
-        WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
-        WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays-1
+        WRITE(lunout, *)'TIME INTERPOLATION.'
+        WRITE(lunout, *)' Input time vector: ', timeyear
+        WRITE(lunout, *)' Output time vector from 0 to ', ndays-1
      END IF
   END IF
@@ -542 +556 @@
   DO j=1, jjp1
     CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
-    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
+    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' at time 10 ', chmin, chmax, j
   END DO
 
 !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    CALL msg(5,'Filtrage de la SST: SST >= 271.38')
+    CALL msg(0,'Filtering SST: SST >= 271.38')
     WHERE(champan<271.38) champan=271.38
   END IF
@@ -553 +567 @@
 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0')
-
-    IF (unit_sic=='1') THEN
+    CALL msg(0,'Filtering SIC: 0.0 < Sea-ice < 1.0')
+    IF(TRIM(unit_sic)=="1".OR.TRIM(unit_sic)=="1.0") THEN
        ! Nothing to be done for sea-ice field is already in fraction of 1
        ! This is the case for sea-ice in file cpl_atm_sic.nc
-       CALL msg(5,'Sea-ice field already in fraction of 1')
+       CALL msg(0,'Sea-ice field already in fraction of 1')
     ELSE
        ! Convert sea ice from percentage to fraction of 1
-       CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.')
+       CALL msg(0,'Sea-ice field converted from percentage to fraction of 1.')
        champan(:, :, :)=champan(:, :, :)/100.
     END IF
-
     champan(iip1, :, :)=champan(1, :, :)
     WHERE(champan>1.0) champan=1.0
@@ -647 +659 @@
 !-------------------------------------------------------------------------------
 !
-FUNCTION year_len(y,cal_in)
-!
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER                       :: year_len
-  INTEGER,           INTENT(IN) :: y
-  CHARACTER(LEN=*),  INTENT(IN) :: cal_in
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
-!-------------------------------------------------------------------------------
-!--- 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 number of days in this year
-  year_len=ioget_year_len(y)
-
-!--- Back to original calendar
-  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
-
-END FUNCTION year_len
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-FUNCTION mid_month(y,cal_in,ndays_out)
-!
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  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           ! 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))
-
-END FUNCTION mid_month
-!
-!-------------------------------------------------------------------------------
-
-
-
-!-------------------------------------------------------------------------------
-!
 SUBROUTINE msg(lev,str1,i,str2)
 !