Changeset 2947 for LMDZ5/trunk/libf/dynphy_lonlat

Timestamp:

Jul 13, 2017, 6:43:58 PM (7 years ago)

Author:

dcugnet

Message:

• The "units" attribute of SST variable is now checked to allow celcius degrees to kelvins conversion if required.
• Few cosmetic changes for output text file to make it more readable (and help to detect if ce0l has really done what was expected).

File:

• LMDZ5/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (modified) (14 diffs)

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

@@ -128 +128 @@
 
 !--- RUGOSITY TREATMENT --------------------------------------------------------
-  CALL msg(1,'Traitement de la rugosite')
+  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE LA RUGOSITE ***")
   CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
 
 !--- OCEAN TREATMENT -----------------------------------------------------------
-  CALL msg(1,'Traitement de la glace oceanique')
+  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE LA GLACE OCEANIQUE ***")
 
 ! Input SIC file selection
@@ -148 +148 @@
   END IF
   CALL ncerr(NF90_CLOSE(nid),icefile)
-  CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
+  CALL msg(0,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
 
   CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
@@ -191 +191 @@
 
 !--- SST TREATMENT -------------------------------------------------------------
-  CALL msg(1,'Traitement de la sst')
+  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE LA SST ***")
 
 ! Input SST file selection
@@ -207 +207 @@
   END IF
   CALL ncerr(NF90_CLOSE(nid),sstfile)
-  CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
+  CALL msg(0,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
 
   CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
 
 !--- ALBEDO TREATMENT ----------------------------------------------------------
-  CALL msg(1,'Traitement de l albedo')
+  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE L'ALBEDO ***")
   CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
 
@@ -219 +219 @@
 
 !--- OUTPUT FILE WRITING -------------------------------------------------------
-  CALL msg(5,'Ecriture du fichier limit : debut')
+  CALL msg(0,""); CALL msg(0,' *** Ecriture du fichier limit : debut ***')
   fnam="limit.nc"
 
@@ -298 +298 @@
   CALL ncerr(NF90_CLOSE(nid),fnam)
 
-  CALL msg(6,'Ecriture du fichier limit : fin')
+  CALL msg(0,""); CALL msg(0,' *** Ecriture du fichier limit : fin ***')
 
   DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
@@ -362 +362 @@
   CHARACTER(LEN=20) :: fnam_m, fnam_p     ! previous/next files names
   CHARACTER(LEN=20) :: cal_in             ! calendar
-  CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
+  CHARACTER(LEN=20) :: unit_sic           ! attribute "units" in sea-ice file
+  CHARACTER(LEN=20) :: unit_sst           ! attribute "units" in sst     file
   INTEGER           :: ndays_in           ! number of days
 !--- misc
@@ -395 +396 @@
   CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
 
-!--- Read unit for sea-ice variable only
+!--- 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(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))
+      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('SIC','Unrecognized unit for sea-ice file: '&
-        //TRIM(unit_sic),1)
+      CALL abort_physic('SST','Unrecognized sst unit: '//TRIM(unit_sst),1)
     END IF
   END IF
@@ -433 +447 @@
       CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
+    CALL msg(0,'WARNING: missing "calendar" attribute for "time" in '&
      &//TRIM(fnam)//'. Choosing default value.')
   END IF
   CALL strclean(cal_in)                     !--- REMOVE (WEIRD) NULL CHARACTERS
-  CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
+  CALL msg(0,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
   
 !--- Determining input file number of days, depending on calendar
@@ -469 +483 @@
         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
+      ELSE IF(TRIM(unit_sic)=="%") THEN
         IF(ANY(champ>100.0+EPSFRA)) &
           CALL abort_physic('SIC','Found sea-ice percentages greater than 100.')
@@ -542 +556 @@
   IF(prt_level>0) THEN
      IF(ndays/=ndays_in) THEN
-        WRITE(lunout, *)
         WRITE(lunout,*)'DIFFERENT YEAR LENGTHS:'
-        WRITE(lunout,*)' In the input file: ',ndays_in
+        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, *)'NO TIME INTERPOLATION.'
         WRITE(lunout, *)' Daily input file.'
@@ -602 +614 @@
 !--- 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')
     WHERE(champan<271.38) champan=271.38
@@ -609 +629 @@
   IF(mode=='SIC') 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
+    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')
@@ -709 +729 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
 !-------------------------------------------------------------------------------
-  IF(prt_level>lev) THEN
+  IF(prt_level>=lev) THEN
     IF(PRESENT(str2)) THEN
       WRITE(lunout,*) TRIM(str1), i, TRIM(str2)