Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 =================================================================== --- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (revision 3125) +++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90 (revision 3163) @@ -27,16 +27,22 @@ 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,7 +368,7 @@ 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,31 +401,4 @@ 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,23 +456,48 @@ 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,5 +571,5 @@ 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,13 +618,10 @@ !--- 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,15 +629,10 @@ !--- SPECIAL FILTER FOR SIC: 0.01.0) champan=1.0 WHERE(champan<0.0) champan=0.0 @@ -789,4 +785,50 @@ !------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! +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