Ignore:
Timestamp:
Sep 2, 2010, 3:45:23 PM (14 years ago)
Author:
lguez
Message:

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90

    r1404 r1425  
    267267#endif
    268268! of #ifdef CPP_EARTH
    269   STOP
    270269
    271270
    … …  
    281280SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
    282281!
    283 !-------------------------------------------------------------------------------
     282!-----------------------------------------------------------------------------
    284283! Comments:
    285284!   There are two assumptions concerning the NetCDF files, that are satisfied
    … …  
    287286!     1) The last dimension of the variables used is the time record.
    288287!     2) Dimensional variables have the same names as corresponding dimensions.
    289 !-------------------------------------------------------------------------------
    290   USE netcdf, ONLY : NF90_OPEN,    NF90_INQ_VARID, NF90_INQUIRE_VARIABLE,      &
    291                      NF90_CLOSE,   NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION,    &
    292                      NF90_GET_VAR, NF90_GET_ATT
     288!-----------------------------------------------------------------------------
     289  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
     290       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
     291       NF90_GET_ATT
    293292  USE dimphy, ONLY : klon
    294293  USE phys_state_var_mod, ONLY : pctsrf
    295294  USE control_mod
     295  use pchsp_95_m, only: pchsp_95
     296  use pchfe_95_m, only: pchfe_95
     297  use arth_m, only: arth
     298
    296299  IMPLICIT NONE
    297300#include "dimensions.h"
    … …  
    300303#include "indicesol.h"
    301304#include "iniprint.h"
    302 !-------------------------------------------------------------------------------
     305!-----------------------------------------------------------------------------
    303306! Arguments:
    304307  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
    … …  
    306309  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
    307310  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
    308   REAL,    POINTER,  DIMENSION(:,:) :: champo   ! output field = f(t)
    309   LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST)  old ice (SIC)
    310   REAL,    OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask
     311  REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
     312  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
     313  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
    311314  LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
    312 !-------------------------------------------------------------------------------
     315!------------------------------------------------------------------------------
    313316! Local variables:
    314317!--- NetCDF
    315   INTEGER :: ierr, ncid, varid                  ! NetCDF identifiers
     318  INTEGER :: ncid, varid                  ! NetCDF identifiers
    316319  CHARACTER(LEN=30)               :: dnam       ! dimension name
    317320  CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
    … …  
    323326!--- fields
    324327  INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
    325   REAL, ALLOCATABLE, DIMENSION(:,:) :: champ      ! wanted field on initial grid
     328  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ   ! wanted field on initial grid
    326329  REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
    327   REAL,              DIMENSION(iim,jjp1) :: champint   ! interpolated field
    328   REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champtime
    329   REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champan
    330   REAL                              :: time, by
     330  REAL,              DIMENSION(iim, jjp1) :: champint   ! interpolated field
     331  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champtime
     332  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champan
    331333!--- input files
    332334  CHARACTER(LEN=20)                 :: cal_in   ! calendar
    … …  
    334336!--- misc
    335337  INTEGER :: i, j, k, l                         ! loop counters
    336   REAL, ALLOCATABLE, DIMENSION(:,:) :: work     ! used for extrapolation
     338  REAL, ALLOCATABLE, DIMENSION(:, :) :: work     ! used for extrapolation
    337339  CHARACTER(LEN=25)                 :: title    ! for messages
    338340  LOGICAL                           :: extrp    ! flag for extrapolation
    339341  REAL                              :: chmin, chmax
    340 !-------------------------------------------------------------------------------
    341 !---Variables depending on keyword 'mode' --------------------------------------
     342  INTEGER ierr
     343  integer n_extrap ! number of extrapolated points
     344  logical skip
     345!------------------------------------------------------------------------------
     346!---Variables depending on keyword 'mode' -------------------------------------
    342347  NULLIFY(champo)
    343348  SELECT CASE(mode)
    … …  
    352357  END IF
    353358
    354 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
    355   ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
    356   ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
    357   ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
     359!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
     360  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
     361  ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
     362  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
    358363
    359364!--- Longitude
    360   ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
    361   CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
    362   ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
    363   ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
    364   WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
     365  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
     366  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
     367  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
     368  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
     369  WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
    365370
    366371!--- Latitude
    367   ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
    368   CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
    369   ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
    370   ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
    371   WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
     372  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
     373  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
     374  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
     375  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
     376  WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
    372377
    373378!--- Time (variable is not needed - it is rebuilt - but calendar is)
    374   ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
    375   CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
    376   ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
     379  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
     380  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
     381  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
    377382  cal_in=' '
    378   ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
     383  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
    379384  IF(ierr/=NF90_NOERR) THEN
    380385    SELECT CASE(mode)
    381       CASE('RUG','ALB'); cal_in='360d'
    382       CASE('SIC','SST'); cal_in='gregorian'
     386      CASE('RUG', 'ALB'); cal_in='360d'
     387      CASE('SIC', 'SST'); cal_in='gregorian'
    383388    END SELECT
    384     WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
    385      &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
     389    WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
     390         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
    386391  END IF
    387   WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
    388 
    389 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
     392  WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
     393       cal_in
     394
     395!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
    390396  !--- Determining input file number of days, depending on calendar
    391   ndays_in=year_len(anneeref,cal_in)
     397  ndays_in=year_len(anneeref, cal_in)
    392398
    393399!--- Time vector reconstruction (time vector from file is not trusted)
    394400!--- If input records are not monthly, time sampling has to be constant !
    395   timeyear=mid_months(anneeref,cal_in,lmdep)
    396   IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode)   &
    397     //' ne comportent pas 12, mais ',lmdep,' enregistrements.'
    398 
    399 !--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
    400   ALLOCATE(champ(imdep,jmdep),champtime(iim,jjp1,lmdep))
    401   IF(extrp) ALLOCATE(work(imdep,jmdep))
    402 
    403   WRITE(lunout,*)
    404   WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
    405   ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
    406   DO l=1,lmdep
    407     ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
    408     CALL ncerr(ierr,fnam)
    409     CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
    410 
    411     IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
     401  timeyear=mid_months(anneeref, cal_in, lmdep)
     402  IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
     403       // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
     404       ' enregistrements.'
     405
     406!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
     407  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
     408  IF(extrp) ALLOCATE(work(imdep, jmdep))
     409
     410  WRITE(lunout, *)
     411  WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
     412       ' CHAMPS.'
     413  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
     414  DO l=1, lmdep
     415    ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
     416    CALL ncerr(ierr, fnam)
     417    CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
     418         champ, ibar)
     419
     420    IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
     421         work)
    412422
    413423    IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
    414424      IF(l==1) THEN
    415         WRITE(lunout,*)                                                        &
     425        WRITE(lunout, *)                                                      &
    416426  '-------------------------------------------------------------------------'
    417         WRITE(lunout,*)                                                        &
    418   '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
    419         WRITE(lunout,*)                                                        &
     427        WRITE(lunout, *)                                                     &
     428  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
     429        WRITE(lunout, *)                                                      &
    420430  '-------------------------------------------------------------------------'
    421431      END IF
    … …  
    433443        CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
    434444                                    iim, jjp1, rlonv, rlatu, champint)
    435         CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
     445        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
    436446                                    iim, jjp1, rlonv, rlatu, champint)
    437447      END SELECT
    438448    END IF
    439     champtime(:,:,l)=champint
     449    champtime(:, :, l)=champint
    440450  END DO
    441   ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
    442 
    443   DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
     451  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
     452
     453  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
    444454  IF(extrp) DEALLOCATE(work)
    445455
    446 !--- TIME INTERPOLATION --------------------------------------------------------
    447   WRITE(lunout,*)
    448   WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
    449   WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
    450   WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
    451   ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
    452   DO j=1,jjp1
    453     DO i=1,iim
    454       CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
    455       DO k=1,ndays
    456         time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
    457         CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
    458         champan(i,j,k) = by
    459       END DO
     456!--- TIME INTERPOLATION ------------------------------------------------------
     457  WRITE(lunout, *)
     458  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
     459  WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
     460  WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
     461  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
     462  skip = .false.
     463  n_extrap = 0
     464  DO j=1, jjp1
     465    DO i=1, iim
     466      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
     467           vc_beg=0., vc_end=0.)
     468      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
     469           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
     470      if (ierr < 0) stop 1
     471      n_extrap = n_extrap + ierr
    460472    END DO
    461473  END DO
    462   champan(iip1,:,:)=champan(1,:,:)
    463   DEALLOCATE(yder,champtime,timeyear)
     474  if (n_extrap /= 0) then
     475     print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
     476  end if
     477  champan(iip1, :, :)=champan(1, :, :)
     478  DEALLOCATE(yder, champtime, timeyear)
    464479
    465480!--- Checking the result
    466   DO j=1,jjp1
    467     CALL minmax(iip1,champan(1,j,10),chmin,chmax)
    468     WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
     481  DO j=1, jjp1
     482    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
     483    WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
    469484  END DO
    470485
    471 !--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
     486!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
    472487  IF(mode=='SST') THEN
    473     WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
     488    WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
    474489    WHERE(champan<271.38) champan=271.38
    475490  END IF
    476491
    477 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
     492!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
    478493  IF(mode=='SIC') THEN
    479     WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
    480     IF(.NOT.lCPL) champan(:,:,:)=champan(:,:,:)/100.
    481     champan(iip1,:,:)=champan(1,:,:)
     494    WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
     495    IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
     496    champan(iip1, :, :)=champan(1, :, :)
    482497    WHERE(champan>1.0) champan=1.0
    483498    WHERE(champan<0.0) champan=0.0
    484499  END IF
    485500
    486 !--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
    487   ALLOCATE(champo(klon,ndays))
    488   DO k=1,ndays
    489     CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k))
     501!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
     502  ALLOCATE(champo(klon, ndays))
     503  DO k=1, ndays
     504    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
    490505  END DO
    491506  DEALLOCATE(champan)
Note: See TracChangeset for help on using the changeset viewer.