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/dyn3dpar/limit_netcdf.F90

    r1404 r1425  
    2020!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
    2121!-------------------------------------------------------------------------------
     22  USE control_mod
    2223#ifdef CPP_EARTH
    2324  USE dimphy
    … …  
    2728                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
    2829                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
    29                    NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED
     30                   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
    3031  USE inter_barxy_m, only: inter_barxy
    3132#endif
    32   USE control_mod
    3333  IMPLICIT NONE
    3434!-------------------------------------------------------------------------------
    … …  
    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)
    … …  
    347352    CASE('ALB'); varname='ALBEDO'; title='Albedo'
    348353  END SELECT
    349   extrp=.FALSE.
    350   IF ( PRESENT(flag) ) THEN
    351     IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
    352   END IF 
    353 
    354 
    355 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
    356   ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
    357   ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
    358   ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
     354  extrp=.FALSE.
     355  IF ( PRESENT(flag) ) THEN
     356    IF ( flag .AND. mode=='SST' ) extrp=.TRUE.
     357  END IF
     358
     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)
    359363
    360364!--- Longitude
    361   ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
    362   CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
    363   ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
    364   ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
    365   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
    366370
    367371!--- Latitude
    368   ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
    369   CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
    370   ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
    371   ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
    372   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
    373377
    374378!--- Time (variable is not needed - it is rebuilt - but calendar is)
    375   ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
    376   CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
    377   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)
    378382  cal_in=' '
    379   ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
     383  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
    380384  IF(ierr/=NF90_NOERR) THEN
    381385    SELECT CASE(mode)
    382       CASE('RUG','ALB'); cal_in='360d'
    383       CASE('SIC','SST'); cal_in='gregorian'
     386      CASE('RUG', 'ALB'); cal_in='360d'
     387      CASE('SIC', 'SST'); cal_in='gregorian'
    384388    END SELECT
    385     WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
    386      &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.'
    387391  END IF
    388   WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
    389 
    390 !--- 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 --------------------
    391396  !--- Determining input file number of days, depending on calendar
    392   ndays_in=year_len(anneeref,cal_in)
     397  ndays_in=year_len(anneeref, cal_in)
    393398
    394399!--- Time vector reconstruction (time vector from file is not trusted)
    395400!--- If input records are not monthly, time sampling has to be constant !
    396   timeyear=mid_months(anneeref,cal_in,lmdep)
    397   IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode)   &
    398     //' ne comportent pas 12, mais ',lmdep,' enregistrements.'
    399 
    400 !--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
    401   ALLOCATE(champ(imdep,jmdep),champtime(iim,jjp1,lmdep))
    402   IF(extrp) ALLOCATE(work(imdep,jmdep))
    403 
    404   WRITE(lunout,*)
    405   WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
    406   ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
    407   DO l=1,lmdep
    408     ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
    409     CALL ncerr(ierr,fnam)
    410     CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
    411 
    412     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)
    413422
    414423    IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
    415424      IF(l==1) THEN
    416         WRITE(lunout,*)                                                        &
     425        WRITE(lunout, *)                                                      &
    417426  '-------------------------------------------------------------------------'
    418         WRITE(lunout,*)                                                        &
    419   '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
    420         WRITE(lunout,*)                                                        &
     427        WRITE(lunout, *)                                                     &
     428  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
     429        WRITE(lunout, *)                                                      &
    421430  '-------------------------------------------------------------------------'
    422431      END IF
    … …  
    434443        CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
    435444                                    iim, jjp1, rlonv, rlatu, champint)
    436         CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
     445        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
    437446                                    iim, jjp1, rlonv, rlatu, champint)
    438447      END SELECT
    439448    END IF
    440     champtime(:,:,l)=champint
     449    champtime(:, :, l)=champint
    441450  END DO
    442   ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
    443 
    444   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)
    445454  IF(extrp) DEALLOCATE(work)
    446455
    447 !--- TIME INTERPOLATION --------------------------------------------------------
    448   WRITE(lunout,*)
    449   WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
    450   WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
    451   WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
    452   ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
    453   DO j=1,jjp1
    454     DO i=1,iim
    455       CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
    456       DO k=1,ndays
    457         time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
    458         CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
    459         champan(i,j,k) = by
    460       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
    461472    END DO
    462473  END DO
    463   champan(iip1,:,:)=champan(1,:,:)
    464   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)
    465479
    466480!--- Checking the result
    467   DO j=1,jjp1
    468     CALL minmax(iip1,champan(1,j,10),chmin,chmax)
    469     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
    470484  END DO
    471485
    472 !--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
     486!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
    473487  IF(mode=='SST') THEN
    474     WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
     488    WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
    475489    WHERE(champan<271.38) champan=271.38
    476490  END IF
    477491
    478 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
     492!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
    479493  IF(mode=='SIC') THEN
    480     WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
    481     IF(.NOT.lCPL) champan(:,:,:)=champan(:,:,:)/100.
    482     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, :, :)
    483497    WHERE(champan>1.0) champan=1.0
    484498    WHERE(champan<0.0) champan=0.0
    485499  END IF
    486500
    487 !--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
    488   ALLOCATE(champo(klon,ndays))
    489   DO k=1,ndays
    490     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))
    491505  END DO
    492506  DEALLOCATE(champan)
Note: See TracChangeset for help on using the changeset viewer.