Changeset 3163

Timestamp:

Jan 26, 2018, 7:11:23 PM (6 years ago)

Author:

dcugnet

Message:

In limit_netcdf.f90:

• add a richer dictionary of recognized SST "units" attributes (taken from udunits).
• in case the unit is not specified, try to determine it from the values in the file: mean zonal temperature at the equator for SST and maximum value for SIC.

File:

• LMDZ6/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (modified) (8 diffs)

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

@@ -27 +27 @@
   USE init_ssrf_m,        ONLY: start_init_subsurf
 
-  CHARACTER(LEN=20), PARAMETER :: &
+  INTEGER,           PARAMETER :: ns=256
+  CHARACTER(LEN=ns), PARAMETER :: &
   fsst(5)=['amipbc_sst_1x1.nc   ','amip_sst_1x1.nc     ','cpl_atm_sst.nc      '&
-          ,'histmth_sst.nc      ','sstk.nc             ']
-  CHARACTER(LEN=20), PARAMETER :: &
+          ,'histmth_sst.nc      ','sstk.nc             '],                     &
   fsic(5)=['amipbc_sic_1x1.nc   ','amip_sic_1x1.nc     ','cpl_atm_sic.nc      '&
-          ,'histmth_sic.nc      ','ci.nc               ']
-  CHARACTER(LEN=10), PARAMETER :: &
+          ,'histmth_sic.nc      ','ci.nc               '],                     &
   vsst(5)=['tosbcs    ','tos       ','SISUTESW  ','tsol_oce  ','sstk      '],  &
-  vsic(5)=['sicbcs    ','sic       ','SIICECOV  ','pourc_sic ','ci        ']
-  CHARACTER(LEN=10), PARAMETER :: &
-  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',    &
-   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    '
+  vsic(5)=['sicbcs    ','sic       ','SIICECOV  ','pourc_sic ','ci        '],  &
+  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',                  &
+   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    ',                  &
+  DegK(11)=['degK          ','degree_K      ','degreeK       ','deg_K         '&
+           ,'degsK         ','degrees_K     ','degreesK      ','degs_K        '&
+           ,'degree_kelvin ','degrees_kelvin','K             '],               &
+  DegC(10)=['degC          ','degree_C      ','degreeC       ','deg_C         '&
+           ,'degsC         ','degrees_C     ','degreesC      ','degs_C        '&
+           ,'degree_Celsius','celsius       '], &
+  Perc(2)=['%              ','percent       '], &
+  Frac(2)=['1.0            ','1             ']
 
 CONTAINS
@@ -362 +368 @@
   CHARACTER(LEN=20) :: fnam_m, fnam_p     ! previous/next files names
   CHARACTER(LEN=20) :: cal_in             ! calendar
-  CHARACTER(LEN=20) :: unit_sic           ! attribute "units" in sea-ice file
-  CHARACTER(LEN=20) :: unit_sst           ! attribute "units" in sst     file
+  CHARACTER(LEN=20) :: units              ! attribute "units" in sic/sst file
   INTEGER           :: ndays_in           ! number of days
+  REAL              :: value              ! mean/max value near equator
 !--- misc
   INTEGER           :: i, j, k, l, ll     ! loop counters
@@ -395 +401 @@
   CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
   CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
-
-!--- Read unit for sea-ice and sst only
-  IF (mode=='SIC') THEN
-    ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic); CALL strclean(unit_sic)
-    IF(ierr/=NF90_NOERR) THEN; unit_sic='%'
-      CALL msg(0,'No unit in sea-ice file. Take percentage as default value')
-    ELSE IF(TRIM(unit_sic)=="%") THEN
-      CALL msg(0,'Sea-ice cover is a PERCENTAGE.')
-    ELSE IF(ANY(unit_sic==["1.0","1  "])) THEN; unit_sic="1"
-      CALL msg(0,'Sea-ice cover is a FRACTION.')
-    ELSE
-      CALL abort_physic('SIC','Unrecognized sea-ice unit: '//TRIM(unit_sic),1)
-    END IF
-  END IF
-  IF (mode=='SST') THEN
-    ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sst); CALL strclean(unit_sst)
-    IF(ierr/=NF90_NOERR) THEN
-      CALL msg(0,'No unit in sst file. Take default: kelvins.')
-      unit_sst='X'
-    ELSE IF(ANY(unit_sst==["degC  ","DegC  "])) THEN; unit_sst='C'
-      CALL msg(0,'Sea-surface temperature is in CELCIUS DEGREES.')
-    ELSE IF(ANY(unit_sst==["K     ","Kelvin"])) THEN; unit_sst='K'
-      CALL msg(0,'Sea-surface temperature is in KELVINS.')
-    ELSE
-      CALL abort_physic('SST','Unrecognized sst unit: '//TRIM(unit_sst),1)
-    END IF
-  END IF
 
 !--- Longitude
@@ -477 +456 @@
   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.
-    !--- Dropped for mid-month boundary conditions datasets (BCS, ix_sic==1)
-    IF(mode=='SIC'.AND.ix_sic/=1) 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)=="%") THEN
-        IF(ANY(champ>100.0+EPSFRA)) &
-          CALL abort_physic('SIC','Found sea-ice percentages greater than 100.')
-!        IF(MAXVAL(champ)< 1.01) &
-!          CALL abort_physic('SIC','All sea-ice percentages lower than 1.')
+    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
+
+    !--- FOR SIC/SST FIELDS ONLY
+    IF(l==1.AND.is_in(mode,['SIC','SST'])) THEN
+
+      !--- DETERMINE THE UNIT: READ FROM FILE OR ASSUMED USING FIELD VALUES
+      ierr=NF90_GET_ATT(ncid, varid, 'units', units)
+      IF(ierr==NF90_NOERR) THEN !--- ATTRIBUTE "units" FOUND IN THE FILE
+        CALL strclean(units)
+        IF(mode=='SIC'.AND.is_in(units,Perc)) units="%"
+        IF(mode=='SIC'.AND.is_in(units,Frac)) units="1"
+        IF(mode=='SST'.AND.is_in(units,DegC)) units="C"
+        IF(mode=='SST'.AND.is_in(units,DegK)) units="K"
+      ELSE                      !--- CHECK THE FIELD VALUES
+        IF(mode=='SIC') value=MAXVAL(champ(:,:))
+        IF(mode=='SST') value=   SUM(champ(:,jmdep/2),DIM=1)/REAL(imdep)
+        IF(mode=='SIC') THEN; units="1"; IF(value>= 10.) units="%"; END IF
+        IF(mode=='SST') THEN; units="C"; IF(value>=100.) units="K"; END IF
       END IF
+      CALL msg(0,'INPUT FILE '//TRIM(title)//' UNIT IS: "'//TRIM(units)//'".')
+      IF(ierr/=NF90_NOERR) CALL msg(0,'WARNING ! UNIT TO BE CHECKED ! '      &
+        //'No "units" attribute, so only based on the fields values.')
+
+      !--- CHECK VALUES ARE IN THE EXPECTED RANGE
+      SELECT CASE(units)
+        CASE('%'); ll=ANY(champ>100.0+EPSFRA); str='percentages > 100.'
+        CASE('1'); ll=ANY(champ>  1.0+EPSFRA); str='fractions > 1.'
+        CASE('C'); ll=ANY(champ<-100.).OR.ANY(champ> 60.); str='<-100 or >60 DegC'
+        CASE('K'); ll=ANY(champ< 180.).OR.ANY(champ>330.); str='<180 or >330 DegK'
+        CASE DEFAULT; CALL abort_physic(mode, 'Unrecognized '//TRIM(title)   &
+                                                  //' unit: '//TRIM(units),1)
+      END SELECT
+
+      !--- DROPPED FOR BCS DATA (FRACTIONS CAN BE HIGHER THAN 1)
+      IF(ll.AND.ix_sic/=1.AND.mode=='SIC') &
+        CALL abort_physic(mode,'unrealistic '//TRIM(mode)//' found: '//TRIM(str))
+
     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,"$$$ Barycentrique interpolation for "//TRIM(title)//" $$$")
-      CALL msg(5,"----------------------------------------------------------")
+      CALL msg(5,"--------------------------------------------------------")
+      CALL msg(5,"$$$ Barycentric interpolation for "//TRIM(title)//" $$$")
+      CALL msg(5,"--------------------------------------------------------")
     END IF
     IF(mode=='RUG') champ=LOG(champ)
@@ -567 +571 @@
         IF(.NOT.is_bcs) WRITE(lunout, *)'SPLINES TIME INTERPOLATION.'
         WRITE(lunout, *)' Input time vector: ', timeyear
-        WRITE(lunout, *)' Output time vector from 0 to ', ndays-1
+        WRITE(lunout, *)' Output time vector: from 0.5 to ', ndays-0.5
      END IF
   END IF
@@ -614 +618 @@
 !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    IF(TRIM(unit_sst)=="K") THEN
-       ! Nothing to be done if the sst field is already in kelvins
-      CALL msg(0,'SST field is already in kelvins.')
-    ELSE
-       ! Convert sst field from celcius degrees to kelvins
-      CALL msg(0,'SST field converted from celcius degrees to kelvins.')
-      champan=champan+273.15
-    END IF
-    CALL msg(0,'Filtering SST: SST >= 271.38')
+    SELECT CASE(units)
+      CASE("K"); CALL msg(0,'SST field is already in kelvins.')
+      CASE("C"); CALL msg(0,'SST field converted from celcius degrees to kelvins.')
+      champan(:, :, :)=champan(:, :, :)+273.15
+    END SELECT
+    CALL msg(0,'Filtering SST: Sea Surface Temperature >= 271.38')
     WHERE(champan<271.38) champan=271.38
   END IF
@@ -628 +629 @@
 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    CALL msg(0,'Filtering SIC: 0.0 < Sea-ice < 1.0')
-    IF(TRIM(unit_sic)=="1") THEN
-       ! Nothing to be done if the sea-ice field is already in fraction of 1
-       ! This is the case for sea-ice in file cpl_atm_sic.nc
-       CALL msg(0,'Sea-ice field already in fraction of 1')
-    ELSE
-       ! Convert sea ice from percentage to fraction of 1
-       CALL msg(0,'Sea-ice field converted from percentage to fraction of 1.')
+    SELECT CASE(units)
+      CASE("1"); CALL msg(0,'SIC field already in fraction of 1')
+      CASE("%"); CALL msg(0,'SIC field converted from percentage to fraction of 1.')
        champan(:, :, :)=champan(:, :, :)/100.
-    END IF
-    champan(iip1, :, :)=champan(1, :, :)
+    END SELECT
+    CALL msg(0,'Filtering SIC: 0.0 <= Sea-ice <=1.0')
     WHERE(champan>1.0) champan=1.0
     WHERE(champan<0.0) champan=0.0
@@ -789 +785 @@
 !-------------------------------------------------------------------------------
 
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION is_in(s1,s2) RESULT(res)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Check wether s1 is present in the s2(:) list (case insensitive).
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: s1, s2(:)
+  LOGICAL                      :: res
+!-------------------------------------------------------------------------------
+  res=.FALSE.; DO k=1,SIZE(s2); res=res.OR.strLow(s1)==strLow(s2(k)); END DO
+
+END FUNCTION is_in
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+ELEMENTAL FUNCTION strLow(s) RESULT(res)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Lower case conversion.
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: s
+  CHARACTER(LEN=ns)            :: res
+!-------------------------------------------------------------------------------
+! Local variable:
+  INTEGER :: k, ix
+!-------------------------------------------------------------------------------
+  res=s
+  DO k=1,LEN(s); ix=IACHAR(s(k:k))
+    IF(64<ix.AND.ix<91) res(k:k)=ACHAR(ix+32)
+  END DO
+
+END FUNCTION strLow
+!
+!-------------------------------------------------------------------------------
+
 #endif
 ! of #ifndef CPP_1D