Changeset 1425 for LMDZ4/trunk


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.

Location:
LMDZ4/trunk
Files:
18 added
4 deleted
14 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/arch/arch-SX8_BRODIE.fcm

    r1279 r1425  
    22%LINK                sxmpif90
    33%AR                  sxar
    4 %MAKE                make
     4%MAKE                sxgmake
    55%FPP_FLAGS           -P -traditional
    66%FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_MATHKEISAN

  • LMDZ4/trunk/arch/arch-linux-32bit.fcm

    r1146 r1425  
    1 %COMPILER            pgf90
    2 %LINK                pgf90
     1%COMPILER            pgf95
     2%LINK                pgf95
    33%AR                  ar
    44%MAKE                make
    … …  
    77%BASE_FFLAGS         
    88%PROD_FFLAGS         -fast
    9 %DEV_FFLAGS          -g
    10 %DEBUG_FFLAGS        -g
     9%DEV_FFLAGS          -g -O1
     10%DEBUG_FFLAGS        -g -O0 -Kieee -Ktrap=fp -Mbounds
    1111%MPI_FFLAGS
    1212%OMP_FFLAGS         
    13 %BASE_LD             -Wl,-Bstatic -L/usr/lib/gcc-lib/i386-linux/2.95.2
     13%BASE_LD             
    1414%MPI_LD
    1515%OMP_LD             

  • LMDZ4/trunk/bld.cfg

    r1327 r1425  
    103103# extension for module output
    104104bld::outfile_ext::mod .mod
     105bld::tool::SHELL   /bin/bash

  • LMDZ4/trunk/libf/dyn3d/abort_gcm.F

    r1279 r1425  
    2626      character(len=*) message
    2727
    28 !      write(lunout,*) 'in abort_gcm'
    29       write(6,*) 'in abort_gcm'
     28      write(lunout,*) 'in abort_gcm'
    3029#ifdef CPP_IOIPSL
    3130      call histclo
    … …  
    3736c     call histclo(4)
    3837c     call histclo(5)
    39 c     write(lunout,*) 'Stopping in ', modname
    40 c     write(lunout,*) 'Reason = ',message
    41 c     if (ierr .eq. 0) then
    42 c       write(lunout,*) 'Everything is cool'
    43 c     else
    44 c       write(lunout,*) 'Houston, we have a problem ', ierr
    45 c     endif
    46       write(6,*) 'Stopping in ', modname
    47       write(6,*) 'Reason = ',message
     38      write(lunout,*) 'Stopping in ', modname
     39      write(lunout,*) 'Reason = ',message
    4840      if (ierr .eq. 0) then
    49         write(6,*) 'Everything is cool'
     41        write(lunout,*) 'Everything is cool'
     42        stop
    5043      else
    51         write(6,*) 'Houston, we have a problem ', ierr
     44        write(lunout,*) 'Houston, we have a problem ', ierr
     45        stop 1
    5246      endif
    53       STOP
    5447      END

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

    r1403 r1425  
    9393#endif
    9494! of #ifndef CPP_EARTH #else
    95   STOP
    9695
    9796END PROGRAM ce0l

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

    r1406 r1425  
    106106  INTEGER :: iflag_thermals_ed, iflag_thermals_optflux
    107107  REAL    :: tau_thermals, solarlong0,  seuil_inversion
    108   INTEGER :: read_climoz ! read ozone climatology 
     108  INTEGER :: read_climoz ! read ozone climatology
    109109!  Allowed values are 0, 1 and 2
    110110!     0: do not read an ozone climatology
    … …  
    140140                   iflag_thermals_ed,iflag_thermals_optflux,            &
    141141                   iflag_coupl,iflag_clos,iflag_wake, read_climoz,      &
    142                    alp_offset )
     142                   alp_offset)
    143143
    144144! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
    … …  
    291291  q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
    292292
    293 !     Parameterization of ozone chemistry:
    294 !     Look for ozone tracer:
    295       i = 1
    296       DO
    297          found = tname(i)=="O3" .OR. tname(i)=="o3"
    298          if (found .or. i == nqtot) exit
    299          i = i + 1
    300       end do
    301       if (found) then
    302          call regr_lat_time_coefoz
    303          call press_coefoz
    304          call regr_pr_o3(p3d, q3d(:, :, :, i))
    305 !     Convert from mole fraction to mass fraction:
    306          q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
    307       end if
     293! Parameterization of ozone chemistry:
     294! Look for ozone tracer:
     295  i = 1
     296  DO
     297    found = tname(i)=="O3" .OR. tname(i)=="o3"
     298    if (found .or. i == nqtot) exit
     299    i = i + 1
     300  end do
     301  if (found) then
     302    call regr_lat_time_coefoz
     303    call press_coefoz
     304    call regr_pr_o3(p3d, q3d(:, :, :, i))
     305!   Convert from mole fraction to mass fraction:
     306    q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
     307  end if
    308308
    309309!--- OZONE CLIMATOLOGY

  • 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)

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

    r1323 r1425  
    643643!-------------------------------------------------------------------------------
    644644!
    645 SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
    646                           lat_in2, pls_in, var3d, ibar)
    647 !
    648 !-------------------------------------------------------------------------------
     645SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
     646     lon_in2, lat_in2, pls_in, var3d, ibar)
     647
     648  use pchsp_95_m, only: pchsp_95
     649  use pchfe_95_m, only: pchfe_95
     650
    649651! Arguments:
    650652  CHARACTER(LEN=*),             INTENT(IN)    :: varname
    … …  
    655657  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
    656658  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
    657   REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
    658   REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
     659  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
     660  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
    659661  LOGICAL,                      INTENT(IN)    :: ibar
    660 !-------------------------------------------------------------------------------
     662!----------------------------------------------------------------------------
    661663! Local variables:
    662664#include "iniprint.h"
    663   LOGICAL               :: check=.TRUE.
    664   REAL                  :: bx, by, chmin, chmax
    665   INTEGER               :: ii, ij, il
    666   REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
     665  LOGICAL:: check=.TRUE., skip
     666  REAL                  chmin, chmax
     667  INTEGER ii, ij, il, ierr
     668  integer n_extrap ! number of extrapolated points
     669  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
    667670  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
    668   REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
    669   INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
    670 !-------------------------------------------------------------------------------
    671   IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
    672   IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
    673   IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
    674 
    675   IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
    676   CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
     671  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
     672
     673!---------------------------------------------------------------------------
     674  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
     675  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
     676       ttm_dyn
     677  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
     678       llm_dyn
     679
     680  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
     681  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
     682       var_ana3d)
    677683
    678684!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
    679   ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
    680   lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
    681   lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
     685  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
     686  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
     687  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
    682688
    683689!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
    684   ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
    685   CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
     690  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
     691  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
    686692                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
    687   DEALLOCATE(lon_ini,lat_ini)
     693  DEALLOCATE(lon_ini, lat_ini)
    688694
    689695!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
    690   ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
    691   DO il=1,llm_dyn
    692     CALL interp_startvar(varname, ibar, il==1,                                 &
    693       iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
    694       lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
     696  DO il=1, llm_dyn
     697    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
     698         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
     699         lon_in2, lat_in2, var_tmp3d(:, :, il))
    695700  END DO
    696   DEALLOCATE(lon_rad,lat_rad)
    697 
    698   ALLOCATE(lind(llm_dyn))
    699   DO il=1,llm_dyn
    700     lind(il) = llm_dyn-il+1
     701  DEALLOCATE(lon_rad, lat_rad)
     702
     703!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
     704  ax = lev_dyn(llm_dyn:1:-1)
     705  skip = .false.
     706  n_extrap = 0
     707  DO ij=1, jml
     708    DO ii=1, iml-1
     709      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
     710      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
     711      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
     712           var3d(ii, ij, lml:1:-1), ierr)
     713      if (ierr < 0) stop 1
     714      n_extrap = n_extrap + ierr
     715    END DO
    701716  END DO
    702 
    703 !--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
    704   ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
    705   DO ij=1,jml
    706     DO ii=1,iml-1
    707       ax(:)=lev_dyn(lind(:))
    708       ay(:)=var_tmp3d(ii,ij,lind(:))
    709       CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
    710       DO il=1,lml
    711         bx=pls_in(ii,ij,il)
    712         CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
    713         var3d(ii,ij,il) = by
    714       END DO
    715     END DO
    716     var3d(iml,ij,:) = var3d(1,ij,:)
     717  if (n_extrap /= 0) then
     718     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
     719  end if
     720  var3d(iml, :, :) = var3d(1, :, :)
     721
     722  DO il=1, lml
     723    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
     724    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
    717725  END DO
    718   DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
    719 
    720   DO il=1,lml
    721     CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
    722     WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
    723   END DO
    724 
    725   RETURN
    726726
    727727END SUBROUTINE start_inter_3d

  • LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F

    r1279 r1425  
    2323C         ierr    = severity of situation ( = 0 normal )
    2424
    25       character (len=*) :: modname
     25      character(len=*) modname
    2626      integer ierr
    27       character (len=*) :: message
     27      character(len=*) message
    2828
    2929      write(lunout,*) 'in abort_gcm'
    … …  
    4545      if (ierr .eq. 0) then
    4646        write(lunout,*) 'Everything is cool'
     47        stop
    4748      else
    4849        write(lunout,*) 'Houston, we have a problem ', ierr
    49       STOP
     50        stop 1
    5051      endif
    5152      END

  • LMDZ4/trunk/libf/dyn3dpar/ce0l.F90

    r1403 r1425  
    55!
    66PROGRAM ce0l
    7 !
    87!-------------------------------------------------------------------------------
    98! Purpose: Calls etat0, creates initial states and limit_netcdf
    … …  
    104103#endif
    105104! of #ifndef CPP_EARTH #else
    106   STOP
    107105
    108106END PROGRAM ce0l

  • LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90

    r1406 r1425  
    1111! Note: This routine is designed to work for Earth
    1212!-------------------------------------------------------------------------------
     13  USE control_mod
    1314#ifdef CPP_EARTH
    1415  USE startvar
    … …  
    2829  USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
    2930#endif
    30   USE control_mod
    3131  IMPLICIT NONE
    3232!-------------------------------------------------------------------------------
    … …  
    114114!-------------------------------------------------------------------------------
    115115  REAL    :: alp_offset
    116   LOGICAL found
     116  logical found
    117117
    118118!--- Constants
    … …  
    505505  wake_cstar(:) = 0.
    506506  wake_fip(:) = 0.
     507  wake_pe = 0.
     508  fm_therm = 0.
     509  entr_therm = 0.
     510  detr_therm = 0.
     511
    507512  CALL fonte_neige_init(run_off_lic_0)
    508513  CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )

  • 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)

  • LMDZ4/trunk/libf/dyn3dpar/startvar.F90

    r1323 r1425  
    643643!-------------------------------------------------------------------------------
    644644!
    645 SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
    646                           lat_in2, pls_in, var3d, ibar)
    647 !
    648 !-------------------------------------------------------------------------------
     645SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
     646     lon_in2, lat_in2, pls_in, var3d, ibar)
     647
     648  use pchsp_95_m, only: pchsp_95
     649  use pchfe_95_m, only: pchfe_95
     650
    649651! Arguments:
    650652  CHARACTER(LEN=*),             INTENT(IN)    :: varname
    … …  
    655657  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
    656658  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
    657   REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
    658   REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
     659  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
     660  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
    659661  LOGICAL,                      INTENT(IN)    :: ibar
    660 !-------------------------------------------------------------------------------
     662!----------------------------------------------------------------------------
    661663! Local variables:
    662664#include "iniprint.h"
    663   LOGICAL               :: check=.TRUE.
    664   REAL                  :: bx, by, chmin, chmax
    665   INTEGER               :: ii, ij, il
    666   REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
     665  LOGICAL:: check=.TRUE., skip
     666  REAL                  chmin, chmax
     667  INTEGER ii, ij, il, ierr
     668  integer n_extrap ! number of extrapolated points
     669  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
    667670  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
    668   REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
    669   INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
    670 !-------------------------------------------------------------------------------
    671   IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
    672   IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
    673   IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
    674 
    675   IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
    676   CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
     671  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
     672
     673!---------------------------------------------------------------------------
     674  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
     675  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
     676       ttm_dyn
     677  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
     678       llm_dyn
     679
     680  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
     681  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
     682       var_ana3d)
    677683
    678684!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
    679   ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
    680   lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
    681   lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
     685  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
     686  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
     687  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
    682688
    683689!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
    684   ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
    685   CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
     690  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
     691  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
    686692                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
    687   DEALLOCATE(lon_ini,lat_ini)
     693  DEALLOCATE(lon_ini, lat_ini)
    688694
    689695!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
    690   ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
    691   DO il=1,llm_dyn
    692     CALL interp_startvar(varname, ibar, il==1,                                 &
    693       iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
    694       lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
     696  DO il=1, llm_dyn
     697    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
     698         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
     699         lon_in2, lat_in2, var_tmp3d(:, :, il))
    695700  END DO
    696   DEALLOCATE(lon_rad,lat_rad)
    697 
    698   ALLOCATE(lind(llm_dyn))
    699   DO il=1,llm_dyn
    700     lind(il) = llm_dyn-il+1
     701  DEALLOCATE(lon_rad, lat_rad)
     702
     703!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
     704  ax = lev_dyn(llm_dyn:1:-1)
     705  skip = .false.
     706  n_extrap = 0
     707  DO ij=1, jml
     708    DO ii=1, iml-1
     709      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
     710      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
     711      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
     712           var3d(ii, ij, lml:1:-1), ierr)
     713      if (ierr < 0) stop 1
     714      n_extrap = n_extrap + ierr
     715    END DO
    701716  END DO
    702 
    703 !--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
    704   ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
    705   DO ij=1,jml
    706     DO ii=1,iml-1
    707       ax(:)=lev_dyn(lind(:))
    708       ay(:)=var_tmp3d(ii,ij,lind(:))
    709       CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
    710       DO il=1,lml
    711         bx=pls_in(ii,ij,il)
    712         CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
    713         var3d(ii,ij,il) = by
    714       END DO
    715     END DO
    716     var3d(iml,ij,:) = var3d(1,ij,:)
     717  if (n_extrap /= 0) then
     718     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
     719  end if
     720  var3d(iml, :, :) = var3d(1, :, :)
     721
     722  DO il=1, lml
     723    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
     724    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
    717725  END DO
    718   DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
    719 
    720   DO il=1,lml
    721     CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
    722     WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
    723   END DO
    724 
    725   RETURN
    726726
    727727END SUBROUTINE start_inter_3d

  • LMDZ4/trunk/libf/phylmd/physiq.F

    r1422 r1425  
    11271127
    11281128      integer, save::  read_climoz ! read ozone climatology
     1129C     (let it keep the default OpenMP shared attribute)
    11291130C     Allowed values are 0, 1 and 2
    11301131C     0: do not read an ozone climatology
    … …  
    11341135
    11351136      integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
     1137C     (let it keep the default OpenMP shared attribute)
    11361138
    11371139      real, pointer, save:: press_climoz(:)
     1140C     (let it keep the default OpenMP shared attribute)
    11381141!     edges of pressure intervals for ozone climatologies, in Pa, in strictly
    11391142!     ascending order
Note: See TracChangeset for help on using the changeset viewer.