Index: LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90 (revision 1507) +++ LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90 (revision 1508) @@ -42,5 +42,5 @@ REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask #ifndef CPP_EARTH - WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' + CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1) #else !------------------------------------------------------------------------------- @@ -52,18 +52,15 @@ #include "indicesol.h" -!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) -------- - LOGICAL, PARAMETER :: fracterre=.TRUE. - !--- INPUT NETCDF FILES NAMES -------------------------------------------------- CHARACTER(LEN=25) :: icefile, sstfile, dumstr CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', & famipsic='amipbc_sic_1x1.nc ', & - fclimsst='amipbc_sst_1x1_clim.nc ', & - fclimsic='amipbc_sic_1x1_clim.nc ', & fcpldsst='cpl_atm_sst.nc ', & fcpldsic='cpl_atm_sic.nc ', & + fhistsst='histmth_sst.nc ', & + fhistsic='histmth_sic.nc ', & frugo ='Rugos.nc ', & falbe ='Albedo.nc ' - + CHARACTER(LEN=10) :: varname !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ REAL, DIMENSION(klon) :: fi_ice, verif @@ -80,5 +77,4 @@ INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC INTEGER :: NF90_FORMAT - LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files INTEGER :: ndays !--- Depending on the output calendar @@ -104,110 +100,97 @@ !--- RUGOSITY TREATMENT -------------------------------------------------------- - WRITE(lunout,*) 'Traitement de la rugosite' - CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) + IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite' + varname='RUGOS' + CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) !--- OCEAN TREATMENT ----------------------------------------------------------- - PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE. + IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique' ! Input SIC file selection - icefile='fake' - IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic) - IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic) - IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic) - SELECT CASE(icefile) - CASE(famipsic); dumstr='Amip.' - CASE(fclimsic); dumstr='Amip climatologique.' - CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. - CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1) - END SELECT +! Open file only to test if available + IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + icefile=TRIM(famipsic) + varname='sicbcs' + ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + icefile=TRIM(fcpldsic) + varname='SIICECOV' + ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + icefile=TRIM(fhistsic) + varname='pourc_sic' + ELSE + WRITE(lunout,*) 'ERROR! No sea-ice input file was found.' + WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic) + CALL abort_gcm('limit_netcdf','No sea-ice file was found',1) + END IF ierr=NF90_CLOSE(nid) - WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr) - - CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL) + IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile) + + CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice) ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) DO k=1,ndays - fi_ice=phy_ice(:,k) - WHERE(fi_ice>=1.0 ) fi_ice=1.0 - WHERE(fi_ice=1.0 ) fi_ice=1.0 + WHERE(fi_ice=1.0-zmasq) - pctsrf_t(:,is_sic,k)=1.0-zmasq - pctsrf_t(:,is_oce,k)=0.0 + pctsrf_t(:,is_sic,k)=1.0-zmasq + pctsrf_t(:,is_oce,k)=0.0 ELSEWHERE - pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) - WHERE(pctsrf_t(:,is_oce,k)0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad - nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) - IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad - ELSE - pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) - WHERE(NINT(pctsrf(:,is_ter))==1) - pctsrf_t(:,is_sic,k)=0. - pctsrf_t(:,is_oce,k)=0. - WHERE(fi_ice>=1.e-5) - pctsrf_t(:,is_lic,k)=fi_ice - pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k) - ELSEWHERE - pctsrf_t(:,is_lic,k)=0.0 - END WHERE - ELSEWHERE - pctsrf_t(:,is_lic,k) = 0.0 - WHERE(fi_ice>=1.e-5) - pctsrf_t(:,is_ter,k)=0.0 - pctsrf_t(:,is_sic,k)=fi_ice - pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k) - ELSEWHERE - pctsrf_t(:,is_sic,k)=0.0 - pctsrf_t(:,is_oce,k)=1.0 - END WHERE - END WHERE - verif=sum(pctsrf_t(:,:,k),dim=2) - nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5) - IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad - END IF + END WHERE + nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) + IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad + nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) + IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad END DO DEALLOCATE(phy_ice) !--- SST TREATMENT ------------------------------------------------------------- - WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE. + IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst' ! Input SST file selection - sstfile='fake' - IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst) - IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst) - IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst) - SELECT CASE(icefile) - CASE(famipsic); dumstr='Amip.' - CASE(fclimsic); dumstr='Amip climatologique.' - CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. - CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1) - END SELECT +! Open file only to test if available + IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + sstfile=TRIM(famipsst) + varname='tosbcs' + ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + sstfile=TRIM(fcpldsst) + varname='SISUTESW' + ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + sstfile=TRIM(fhistsst) + varname='tsol_oce' + ELSE + WRITE(lunout,*) 'ERROR! No sst input file was found.' + WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst) + CALL abort_gcm('limit_netcdf','No sst file was found',1) + END IF ierr=NF90_CLOSE(nid) - WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr) - - CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap) + IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile) + + CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap) !--- ALBEDO TREATMENT ---------------------------------------------------------- - WRITE(lunout,*) 'Traitement de l albedo' - CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb) + IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo' + varname='ALBEDO' + CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb) !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- @@ -215,5 +198,5 @@ !--- OUTPUT FILE WRITING ------------------------------------------------------- - WRITE(lunout,*) 'Ecriture du fichier limit' + IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut' !--- File creation @@ -264,4 +247,6 @@ ierr=NF90_CLOSE(nid) + IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin' + DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) @@ -276,5 +261,5 @@ !------------------------------------------------------------------------------- ! -SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL) +SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask) ! !----------------------------------------------------------------------------- @@ -304,11 +289,11 @@ ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name + CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels INTEGER, INTENT(IN) :: ndays ! current year number of days - REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) + REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask - LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE) !------------------------------------------------------------------------------ ! Local variables: @@ -316,5 +301,4 @@ INTEGER :: ncid, varid ! NetCDF identifiers CHARACTER(LEN=30) :: dnam ! dimension name - CHARACTER(LEN=80) :: varname ! NetCDF variable name !--- dimensions INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers @@ -331,4 +315,5 @@ !--- input files CHARACTER(LEN=20) :: cal_in ! calendar + CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file INTEGER :: ndays_in ! number of days !--- misc @@ -337,26 +322,45 @@ CHARACTER(LEN=25) :: title ! for messages LOGICAL :: extrp ! flag for extrapolation + LOGICAL :: oldice ! flag for old way ice computation REAL :: chmin, chmax INTEGER ierr integer n_extrap ! number of extrapolated points logical skip + !------------------------------------------------------------------------------ !---Variables depending on keyword 'mode' ------------------------------------- NULLIFY(champo) + SELECT CASE(mode) - CASE('RUG'); varname='RUGOS'; title='Rugosite' - CASE('SIC'); varname='sicbcs'; title='Sea-ice' - CASE('SST'); varname='tosbcs'; title='SST' - CASE('ALB'); varname='ALBEDO'; title='Albedo' + CASE('RUG'); title='Rugosite' + CASE('SIC'); title='Sea-ice' + CASE('SST'); title='SST' + CASE('ALB'); title='Albedo' END SELECT + + extrp=.FALSE. + oldice=.FALSE. IF ( PRESENT(flag) ) THEN IF ( flag .AND. mode=='SST' ) extrp=.TRUE. + IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. END IF !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- + IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam) - ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) + ierr=NF90_INQ_VARID(ncid, trim(varname), varid); CALL ncerr(ierr, fnam) ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam) + +!--- Read unit for sea-ice variable only + IF (mode=='SIC') THEN + ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic) + IF(ierr/=NF90_NOERR) THEN + IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value' + unit_sic='X' + ELSE + IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic + END IF + END IF !--- Longitude @@ -365,5 +369,5 @@ ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam) - WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep + IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep !--- Latitude @@ -372,5 +376,5 @@ ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam) - WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep + IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep !--- Time (variable is not needed - it is rebuilt - but calendar is) @@ -385,10 +389,11 @@ CASE('SIC', 'SST'); cal_in='gregorian' END SELECT - WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & + IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' END IF - WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & + IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & cal_in + !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- !--- Determining input file number of days, depending on calendar @@ -398,7 +403,6 @@ !--- If input records are not monthly, time sampling has to be constant ! timeyear=mid_months(anneeref, cal_in, lmdep) - IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' & - // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, & - ' enregistrements.' + IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), & + ' ne comportent pas 12, mais ', lmdep, ' enregistrements.' !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- @@ -406,7 +410,6 @@ IF(extrp) ALLOCATE(work(imdep, jmdep)) - WRITE(lunout, *) - WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, & - ' CHAMPS.' + IF (prt_level>5) WRITE(lunout, *) + IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.' ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) DO l=1, lmdep @@ -419,12 +422,9 @@ work) - IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN - IF(l==1) THEN - WRITE(lunout, *) & - '-------------------------------------------------------------------------' - WRITE(lunout, *) & - 'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' - WRITE(lunout, *) & - '-------------------------------------------------------------------------' + IF(ibar .AND. .NOT.oldice) THEN + IF(l==1 .AND. prt_level>5) THEN + WRITE(lunout, *) '-------------------------------------------------------------------------' + WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$' + WRITE(lunout, *) '-------------------------------------------------------------------------' END IF IF(mode=='RUG') champ=LOG(champ) @@ -453,8 +453,11 @@ !--- TIME INTERPOLATION ------------------------------------------------------ - WRITE(lunout, *) - WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' - WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear - WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays + IF (prt_level>5) THEN + WRITE(lunout, *) + WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' + WRITE(lunout, *)' Vecteur temps en entree: ', timeyear + WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays + END IF + ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) skip = .false. @@ -471,5 +474,5 @@ END DO if (n_extrap /= 0) then - print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap + WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap end if champan(iip1, :, :)=champan(1, :, :) @@ -479,10 +482,10 @@ DO j=1, jjp1 CALL minmax(iip1, champan(1, j, 10), chmin, chmax) - WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j + IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j END DO !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- IF(mode=='SST') THEN - WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' + IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' WHERE(champan<271.38) champan=271.38 END IF @@ -490,10 +493,20 @@ !--- SPECIAL FILTER FOR SIC: 0.05) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' + + IF (unit_sic=='1') 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 + IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1' + ELSE + ! Convert sea ice from percentage to fraction of 1 + IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' + champan(:, :, :)=champan(:, :, :)/100. + END IF + champan(iip1, :, :)=champan(1, :, :) WHERE(champan>1.0) champan=1.0 WHERE(champan<0.0) champan=0.0 - END IF + END IF !--- DYNAMICAL TO PHYSICAL GRID ---------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90 (revision 1507) +++ LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90 (revision 1508) @@ -42,5 +42,5 @@ REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask #ifndef CPP_EARTH - WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' + CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1) #else !------------------------------------------------------------------------------- @@ -52,18 +52,15 @@ #include "indicesol.h" -!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) -------- - LOGICAL, PARAMETER :: fracterre=.TRUE. - !--- INPUT NETCDF FILES NAMES -------------------------------------------------- CHARACTER(LEN=25) :: icefile, sstfile, dumstr CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', & famipsic='amipbc_sic_1x1.nc ', & - fclimsst='amipbc_sst_1x1_clim.nc ', & - fclimsic='amipbc_sic_1x1_clim.nc ', & fcpldsst='cpl_atm_sst.nc ', & fcpldsic='cpl_atm_sic.nc ', & + fhistsst='histmth_sst.nc ', & + fhistsic='histmth_sic.nc ', & frugo ='Rugos.nc ', & falbe ='Albedo.nc ' - + CHARACTER(LEN=10) :: varname !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ REAL, DIMENSION(klon) :: fi_ice, verif @@ -80,5 +77,4 @@ INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC INTEGER :: NF90_FORMAT - LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files INTEGER :: ndays !--- Depending on the output calendar @@ -104,110 +100,97 @@ !--- RUGOSITY TREATMENT -------------------------------------------------------- - WRITE(lunout,*) 'Traitement de la rugosite' - CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) + IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite' + varname='RUGOS' + CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) !--- OCEAN TREATMENT ----------------------------------------------------------- - PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE. + IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique' ! Input SIC file selection - icefile='fake' - IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic) - IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic) - IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic) - SELECT CASE(icefile) - CASE(famipsic); dumstr='Amip.' - CASE(fclimsic); dumstr='Amip climatologique.' - CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. - CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1) - END SELECT +! Open file only to test if available + IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + icefile=TRIM(famipsic) + varname='sicbcs' + ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + icefile=TRIM(fcpldsic) + varname='SIICECOV' + ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + icefile=TRIM(fhistsic) + varname='pourc_sic' + ELSE + WRITE(lunout,*) 'ERROR! No sea-ice input file was found.' + WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic) + CALL abort_gcm('limit_netcdf','No sea-ice file was found',1) + END IF ierr=NF90_CLOSE(nid) - WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr) - - CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL) + IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile) + + CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice) ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) DO k=1,ndays - fi_ice=phy_ice(:,k) - WHERE(fi_ice>=1.0 ) fi_ice=1.0 - WHERE(fi_ice=1.0 ) fi_ice=1.0 + WHERE(fi_ice=1.0-zmasq) - pctsrf_t(:,is_sic,k)=1.0-zmasq - pctsrf_t(:,is_oce,k)=0.0 + pctsrf_t(:,is_sic,k)=1.0-zmasq + pctsrf_t(:,is_oce,k)=0.0 ELSEWHERE - pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) - WHERE(pctsrf_t(:,is_oce,k)0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad - nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) - IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad - ELSE - pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) - WHERE(NINT(pctsrf(:,is_ter))==1) - pctsrf_t(:,is_sic,k)=0. - pctsrf_t(:,is_oce,k)=0. - WHERE(fi_ice>=1.e-5) - pctsrf_t(:,is_lic,k)=fi_ice - pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k) - ELSEWHERE - pctsrf_t(:,is_lic,k)=0.0 - END WHERE - ELSEWHERE - pctsrf_t(:,is_lic,k) = 0.0 - WHERE(fi_ice>=1.e-5) - pctsrf_t(:,is_ter,k)=0.0 - pctsrf_t(:,is_sic,k)=fi_ice - pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k) - ELSEWHERE - pctsrf_t(:,is_sic,k)=0.0 - pctsrf_t(:,is_oce,k)=1.0 - END WHERE - END WHERE - verif=sum(pctsrf_t(:,:,k),dim=2) - nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5) - IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad - END IF + END WHERE + nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) + IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad + nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) + IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad END DO DEALLOCATE(phy_ice) !--- SST TREATMENT ------------------------------------------------------------- - WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE. + IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst' ! Input SST file selection - sstfile='fake' - IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst) - IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst) - IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst) - SELECT CASE(icefile) - CASE(famipsic); dumstr='Amip.' - CASE(fclimsic); dumstr='Amip climatologique.' - CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. - CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1) - END SELECT +! Open file only to test if available + IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + sstfile=TRIM(famipsst) + varname='tosbcs' + ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + sstfile=TRIM(fcpldsst) + varname='SISUTESW' + ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN + sstfile=TRIM(fhistsst) + varname='tsol_oce' + ELSE + WRITE(lunout,*) 'ERROR! No sst input file was found.' + WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst) + CALL abort_gcm('limit_netcdf','No sst file was found',1) + END IF ierr=NF90_CLOSE(nid) - WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr) - - CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap) + IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile) + + CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap) !--- ALBEDO TREATMENT ---------------------------------------------------------- - WRITE(lunout,*) 'Traitement de l albedo' - CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb) + IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo' + varname='ALBEDO' + CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb) !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- @@ -215,5 +198,5 @@ !--- OUTPUT FILE WRITING ------------------------------------------------------- - WRITE(lunout,*) 'Ecriture du fichier limit' + IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut' !--- File creation @@ -264,4 +247,6 @@ ierr=NF90_CLOSE(nid) + IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin' + DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) @@ -276,5 +261,5 @@ !------------------------------------------------------------------------------- ! -SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL) +SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask) ! !----------------------------------------------------------------------------- @@ -304,11 +289,11 @@ ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name + CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels INTEGER, INTENT(IN) :: ndays ! current year number of days - REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) + REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask - LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE) !------------------------------------------------------------------------------ ! Local variables: @@ -316,5 +301,4 @@ INTEGER :: ncid, varid ! NetCDF identifiers CHARACTER(LEN=30) :: dnam ! dimension name - CHARACTER(LEN=80) :: varname ! NetCDF variable name !--- dimensions INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers @@ -331,4 +315,5 @@ !--- input files CHARACTER(LEN=20) :: cal_in ! calendar + CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file INTEGER :: ndays_in ! number of days !--- misc @@ -337,26 +322,45 @@ CHARACTER(LEN=25) :: title ! for messages LOGICAL :: extrp ! flag for extrapolation + LOGICAL :: oldice ! flag for old way ice computation REAL :: chmin, chmax INTEGER ierr integer n_extrap ! number of extrapolated points logical skip + !------------------------------------------------------------------------------ !---Variables depending on keyword 'mode' ------------------------------------- NULLIFY(champo) + SELECT CASE(mode) - CASE('RUG'); varname='RUGOS'; title='Rugosite' - CASE('SIC'); varname='sicbcs'; title='Sea-ice' - CASE('SST'); varname='tosbcs'; title='SST' - CASE('ALB'); varname='ALBEDO'; title='Albedo' + CASE('RUG'); title='Rugosite' + CASE('SIC'); title='Sea-ice' + CASE('SST'); title='SST' + CASE('ALB'); title='Albedo' END SELECT + + extrp=.FALSE. + oldice=.FALSE. IF ( PRESENT(flag) ) THEN IF ( flag .AND. mode=='SST' ) extrp=.TRUE. + IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. END IF !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- + IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam) - ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) + ierr=NF90_INQ_VARID(ncid, trim(varname), varid); CALL ncerr(ierr, fnam) ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam) + +!--- Read unit for sea-ice variable only + IF (mode=='SIC') THEN + ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic) + IF(ierr/=NF90_NOERR) THEN + IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value' + unit_sic='X' + ELSE + IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic + END IF + END IF !--- Longitude @@ -365,5 +369,5 @@ ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam) - WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep + IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep !--- Latitude @@ -372,5 +376,5 @@ ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam) - WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep + IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep !--- Time (variable is not needed - it is rebuilt - but calendar is) @@ -385,10 +389,11 @@ CASE('SIC', 'SST'); cal_in='gregorian' END SELECT - WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & + IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' END IF - WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & + IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & cal_in + !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- !--- Determining input file number of days, depending on calendar @@ -398,7 +403,6 @@ !--- If input records are not monthly, time sampling has to be constant ! timeyear=mid_months(anneeref, cal_in, lmdep) - IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' & - // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, & - ' enregistrements.' + IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), & + ' ne comportent pas 12, mais ', lmdep, ' enregistrements.' !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- @@ -406,7 +410,6 @@ IF(extrp) ALLOCATE(work(imdep, jmdep)) - WRITE(lunout, *) - WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, & - ' CHAMPS.' + IF (prt_level>5) WRITE(lunout, *) + IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.' ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) DO l=1, lmdep @@ -419,12 +422,9 @@ work) - IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN - IF(l==1) THEN - WRITE(lunout, *) & - '-------------------------------------------------------------------------' - WRITE(lunout, *) & - 'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' - WRITE(lunout, *) & - '-------------------------------------------------------------------------' + IF(ibar .AND. .NOT.oldice) THEN + IF(l==1 .AND. prt_level>5) THEN + WRITE(lunout, *) '-------------------------------------------------------------------------' + WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$' + WRITE(lunout, *) '-------------------------------------------------------------------------' END IF IF(mode=='RUG') champ=LOG(champ) @@ -453,8 +453,11 @@ !--- TIME INTERPOLATION ------------------------------------------------------ - WRITE(lunout, *) - WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' - WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear - WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays + IF (prt_level>5) THEN + WRITE(lunout, *) + WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' + WRITE(lunout, *)' Vecteur temps en entree: ', timeyear + WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays + END IF + ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) skip = .false. @@ -471,5 +474,5 @@ END DO if (n_extrap /= 0) then - print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap + WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap end if champan(iip1, :, :)=champan(1, :, :) @@ -479,10 +482,10 @@ DO j=1, jjp1 CALL minmax(iip1, champan(1, j, 10), chmin, chmax) - WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j + IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j END DO !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- IF(mode=='SST') THEN - WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' + IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' WHERE(champan<271.38) champan=271.38 END IF @@ -490,10 +493,20 @@ !--- SPECIAL FILTER FOR SIC: 0.05) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' + + IF (unit_sic=='1') 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 + IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1' + ELSE + ! Convert sea ice from percentage to fraction of 1 + IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' + champan(:, :, :)=champan(:, :, :)/100. + END IF + champan(iip1, :, :)=champan(1, :, :) WHERE(champan>1.0) champan=1.0 WHERE(champan<0.0) champan=0.0 - END IF + END IF !--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------