Changeset 3163 for LMDZ6/trunk/libf
- Timestamp:
- Jan 26, 2018, 7:11:23 PM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dynphy_lonlat/phylmd/limit_netcdf.F90
r2947 r3163 27 27 USE init_ssrf_m, ONLY: start_init_subsurf 28 28 29 CHARACTER(LEN=20), PARAMETER :: & 29 INTEGER, PARAMETER :: ns=256 30 CHARACTER(LEN=ns), PARAMETER :: & 30 31 fsst(5)=['amipbc_sst_1x1.nc ','amip_sst_1x1.nc ','cpl_atm_sst.nc '& 31 ,'histmth_sst.nc ','sstk.nc '] 32 CHARACTER(LEN=20), PARAMETER :: & 32 ,'histmth_sst.nc ','sstk.nc '], & 33 33 fsic(5)=['amipbc_sic_1x1.nc ','amip_sic_1x1.nc ','cpl_atm_sic.nc '& 34 ,'histmth_sic.nc ','ci.nc '] 35 CHARACTER(LEN=10), PARAMETER :: & 34 ,'histmth_sic.nc ','ci.nc '], & 36 35 vsst(5)=['tosbcs ','tos ','SISUTESW ','tsol_oce ','sstk '], & 37 vsic(5)=['sicbcs ','sic ','SIICECOV ','pourc_sic ','ci '] 38 CHARACTER(LEN=10), PARAMETER :: & 39 frugo='Rugos.nc ', falbe='Albedo.nc ', frelf='Relief.nc ', & 40 vrug='RUGOS ', valb='ALBEDO ', vrel='RELIEF ' 36 vsic(5)=['sicbcs ','sic ','SIICECOV ','pourc_sic ','ci '], & 37 frugo='Rugos.nc ', falbe='Albedo.nc ', frelf='Relief.nc ', & 38 vrug='RUGOS ', valb='ALBEDO ', vrel='RELIEF ', & 39 DegK(11)=['degK ','degree_K ','degreeK ','deg_K '& 40 ,'degsK ','degrees_K ','degreesK ','degs_K '& 41 ,'degree_kelvin ','degrees_kelvin','K '], & 42 DegC(10)=['degC ','degree_C ','degreeC ','deg_C '& 43 ,'degsC ','degrees_C ','degreesC ','degs_C '& 44 ,'degree_Celsius','celsius '], & 45 Perc(2)=['% ','percent '], & 46 Frac(2)=['1.0 ','1 '] 41 47 42 48 CONTAINS … … 362 368 CHARACTER(LEN=20) :: fnam_m, fnam_p ! previous/next files names 363 369 CHARACTER(LEN=20) :: cal_in ! calendar 364 CHARACTER(LEN=20) :: unit_sic ! attribute "units" in sea-ice file 365 CHARACTER(LEN=20) :: unit_sst ! attribute "units" in sst file 370 CHARACTER(LEN=20) :: units ! attribute "units" in sic/sst file 366 371 INTEGER :: ndays_in ! number of days 372 REAL :: value ! mean/max value near equator 367 373 !--- misc 368 374 INTEGER :: i, j, k, l, ll ! loop counters … … 395 401 CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam) 396 402 CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam) 397 398 !--- Read unit for sea-ice and sst only399 IF (mode=='SIC') THEN400 ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic); CALL strclean(unit_sic)401 IF(ierr/=NF90_NOERR) THEN; unit_sic='%'402 CALL msg(0,'No unit in sea-ice file. Take percentage as default value')403 ELSE IF(TRIM(unit_sic)=="%") THEN404 CALL msg(0,'Sea-ice cover is a PERCENTAGE.')405 ELSE IF(ANY(unit_sic==["1.0","1 "])) THEN; unit_sic="1"406 CALL msg(0,'Sea-ice cover is a FRACTION.')407 ELSE408 CALL abort_physic('SIC','Unrecognized sea-ice unit: '//TRIM(unit_sic),1)409 END IF410 END IF411 IF (mode=='SST') THEN412 ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sst); CALL strclean(unit_sst)413 IF(ierr/=NF90_NOERR) THEN414 CALL msg(0,'No unit in sst file. Take default: kelvins.')415 unit_sst='X'416 ELSE IF(ANY(unit_sst==["degC ","DegC "])) THEN; unit_sst='C'417 CALL msg(0,'Sea-surface temperature is in CELCIUS DEGREES.')418 ELSE IF(ANY(unit_sst==["K ","Kelvin"])) THEN; unit_sst='K'419 CALL msg(0,'Sea-surface temperature is in KELVINS.')420 ELSE421 CALL abort_physic('SST','Unrecognized sst unit: '//TRIM(unit_sst),1)422 END IF423 END IF424 403 425 404 !--- Longitude … … 477 456 DO l=1, lmdep 478 457 CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam) 479 !--- Check whether values are acceptable for SIC, depending on unit. 480 !--- Dropped for mid-month boundary conditions datasets (BCS, ix_sic==1) 481 IF(mode=='SIC'.AND.ix_sic/=1) THEN 482 IF(TRIM(unit_sic)=="1".OR.TRIM(unit_sic)=="1.0") THEN 483 IF(ANY(champ>1.0+EPSFRA)) & 484 CALL abort_physic('SIC','Found sea-ice fractions greater than 1.') 485 ELSE IF(TRIM(unit_sic)=="%") THEN 486 IF(ANY(champ>100.0+EPSFRA)) & 487 CALL abort_physic('SIC','Found sea-ice percentages greater than 100.') 488 ! IF(MAXVAL(champ)< 1.01) & 489 ! CALL abort_physic('SIC','All sea-ice percentages lower than 1.') 458 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.) 459 460 !--- FOR SIC/SST FIELDS ONLY 461 IF(l==1.AND.is_in(mode,['SIC','SST'])) THEN 462 463 !--- DETERMINE THE UNIT: READ FROM FILE OR ASSUMED USING FIELD VALUES 464 ierr=NF90_GET_ATT(ncid, varid, 'units', units) 465 IF(ierr==NF90_NOERR) THEN !--- ATTRIBUTE "units" FOUND IN THE FILE 466 CALL strclean(units) 467 IF(mode=='SIC'.AND.is_in(units,Perc)) units="%" 468 IF(mode=='SIC'.AND.is_in(units,Frac)) units="1" 469 IF(mode=='SST'.AND.is_in(units,DegC)) units="C" 470 IF(mode=='SST'.AND.is_in(units,DegK)) units="K" 471 ELSE !--- CHECK THE FIELD VALUES 472 IF(mode=='SIC') value=MAXVAL(champ(:,:)) 473 IF(mode=='SST') value= SUM(champ(:,jmdep/2),DIM=1)/REAL(imdep) 474 IF(mode=='SIC') THEN; units="1"; IF(value>= 10.) units="%"; END IF 475 IF(mode=='SST') THEN; units="C"; IF(value>=100.) units="K"; END IF 490 476 END IF 477 CALL msg(0,'INPUT FILE '//TRIM(title)//' UNIT IS: "'//TRIM(units)//'".') 478 IF(ierr/=NF90_NOERR) CALL msg(0,'WARNING ! UNIT TO BE CHECKED ! ' & 479 //'No "units" attribute, so only based on the fields values.') 480 481 !--- CHECK VALUES ARE IN THE EXPECTED RANGE 482 SELECT CASE(units) 483 CASE('%'); ll=ANY(champ>100.0+EPSFRA); str='percentages > 100.' 484 CASE('1'); ll=ANY(champ> 1.0+EPSFRA); str='fractions > 1.' 485 CASE('C'); ll=ANY(champ<-100.).OR.ANY(champ> 60.); str='<-100 or >60 DegC' 486 CASE('K'); ll=ANY(champ< 180.).OR.ANY(champ>330.); str='<180 or >330 DegK' 487 CASE DEFAULT; CALL abort_physic(mode, 'Unrecognized '//TRIM(title) & 488 //' unit: '//TRIM(units),1) 489 END SELECT 490 491 !--- DROPPED FOR BCS DATA (FRACTIONS CAN BE HIGHER THAN 1) 492 IF(ll.AND.ix_sic/=1.AND.mode=='SIC') & 493 CALL abort_physic(mode,'unrealistic '//TRIM(mode)//' found: '//TRIM(str)) 494 491 495 END IF 492 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.) 496 493 497 IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) 494 498 IF(l==1) THEN 495 CALL msg(5,"-------------------------------------------------------- --")496 CALL msg(5,"$$$ Barycentri queinterpolation for "//TRIM(title)//" $$$")497 CALL msg(5,"-------------------------------------------------------- --")499 CALL msg(5,"--------------------------------------------------------") 500 CALL msg(5,"$$$ Barycentric interpolation for "//TRIM(title)//" $$$") 501 CALL msg(5,"--------------------------------------------------------") 498 502 END IF 499 503 IF(mode=='RUG') champ=LOG(champ) … … 567 571 IF(.NOT.is_bcs) WRITE(lunout, *)'SPLINES TIME INTERPOLATION.' 568 572 WRITE(lunout, *)' Input time vector: ', timeyear 569 WRITE(lunout, *)' Output time vector from 0 to ', ndays-1573 WRITE(lunout, *)' Output time vector: from 0.5 to ', ndays-0.5 570 574 END IF 571 575 END IF … … 614 618 !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- 615 619 IF(mode=='SST') THEN 616 IF(TRIM(unit_sst)=="K") THEN 617 ! Nothing to be done if the sst field is already in kelvins 618 CALL msg(0,'SST field is already in kelvins.') 619 ELSE 620 ! Convert sst field from celcius degrees to kelvins 621 CALL msg(0,'SST field converted from celcius degrees to kelvins.') 622 champan=champan+273.15 623 END IF 624 CALL msg(0,'Filtering SST: SST >= 271.38') 620 SELECT CASE(units) 621 CASE("K"); CALL msg(0,'SST field is already in kelvins.') 622 CASE("C"); CALL msg(0,'SST field converted from celcius degrees to kelvins.') 623 champan(:, :, :)=champan(:, :, :)+273.15 624 END SELECT 625 CALL msg(0,'Filtering SST: Sea Surface Temperature >= 271.38') 625 626 WHERE(champan<271.38) champan=271.38 626 627 END IF … … 628 629 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ------------------------------------- 629 630 IF(mode=='SIC') THEN 630 CALL msg(0,'Filtering SIC: 0.0 < Sea-ice < 1.0') 631 IF(TRIM(unit_sic)=="1") THEN 632 ! Nothing to be done if the sea-ice field is already in fraction of 1 633 ! This is the case for sea-ice in file cpl_atm_sic.nc 634 CALL msg(0,'Sea-ice field already in fraction of 1') 635 ELSE 636 ! Convert sea ice from percentage to fraction of 1 637 CALL msg(0,'Sea-ice field converted from percentage to fraction of 1.') 631 SELECT CASE(units) 632 CASE("1"); CALL msg(0,'SIC field already in fraction of 1') 633 CASE("%"); CALL msg(0,'SIC field converted from percentage to fraction of 1.') 638 634 champan(:, :, :)=champan(:, :, :)/100. 639 END IF640 champan(iip1, :, :)=champan(1, :, :)635 END SELECT 636 CALL msg(0,'Filtering SIC: 0.0 <= Sea-ice <=1.0') 641 637 WHERE(champan>1.0) champan=1.0 642 638 WHERE(champan<0.0) champan=0.0 … … 789 785 !------------------------------------------------------------------------------- 790 786 787 788 !------------------------------------------------------------------------------- 789 ! 790 FUNCTION is_in(s1,s2) RESULT(res) 791 ! 792 !------------------------------------------------------------------------------- 793 IMPLICIT NONE 794 !------------------------------------------------------------------------------- 795 ! Purpose: Check wether s1 is present in the s2(:) list (case insensitive). 796 !------------------------------------------------------------------------------- 797 ! Arguments: 798 CHARACTER(LEN=*), INTENT(IN) :: s1, s2(:) 799 LOGICAL :: res 800 !------------------------------------------------------------------------------- 801 res=.FALSE.; DO k=1,SIZE(s2); res=res.OR.strLow(s1)==strLow(s2(k)); END DO 802 803 END FUNCTION is_in 804 ! 805 !------------------------------------------------------------------------------- 806 807 808 !------------------------------------------------------------------------------- 809 ! 810 ELEMENTAL FUNCTION strLow(s) RESULT(res) 811 ! 812 !------------------------------------------------------------------------------- 813 IMPLICIT NONE 814 !------------------------------------------------------------------------------- 815 ! Purpose: Lower case conversion. 816 !------------------------------------------------------------------------------- 817 ! Arguments: 818 CHARACTER(LEN=*), INTENT(IN) :: s 819 CHARACTER(LEN=ns) :: res 820 !------------------------------------------------------------------------------- 821 ! Local variable: 822 INTEGER :: k, ix 823 !------------------------------------------------------------------------------- 824 res=s 825 DO k=1,LEN(s); ix=IACHAR(s(k:k)) 826 IF(64<ix.AND.ix<91) res(k:k)=ACHAR(ix+32) 827 END DO 828 829 END FUNCTION strLow 830 ! 831 !------------------------------------------------------------------------------- 832 791 833 #endif 792 834 ! of #ifndef CPP_1D
Note: See TracChangeset
for help on using the changeset viewer.
