Index: LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F (revision 1421) +++ LMDZ4/trunk/libf/dyn3dpar/abort_gcm.F (revision 1425) @@ -23,7 +23,7 @@ C ierr = severity of situation ( = 0 normal ) - character (len=*) :: modname + character(len=*) modname integer ierr - character (len=*) :: message + character(len=*) message write(lunout,*) 'in abort_gcm' @@ -45,7 +45,8 @@ if (ierr .eq. 0) then write(lunout,*) 'Everything is cool' + stop else write(lunout,*) 'Houston, we have a problem ', ierr - STOP + stop 1 endif END Index: LMDZ4/trunk/libf/dyn3dpar/ce0l.F90 =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/ce0l.F90 (revision 1421) +++ LMDZ4/trunk/libf/dyn3dpar/ce0l.F90 (revision 1425) @@ -5,5 +5,4 @@ ! PROGRAM ce0l -! !------------------------------------------------------------------------------- ! Purpose: Calls etat0, creates initial states and limit_netcdf @@ -104,5 +103,4 @@ #endif ! of #ifndef CPP_EARTH #else - STOP END PROGRAM ce0l Index: LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90 =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90 (revision 1421) +++ LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90 (revision 1425) @@ -11,4 +11,5 @@ ! Note: This routine is designed to work for Earth !------------------------------------------------------------------------------- + USE control_mod #ifdef CPP_EARTH USE startvar @@ -28,5 +29,4 @@ USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR #endif - USE control_mod IMPLICIT NONE !------------------------------------------------------------------------------- @@ -114,5 +114,5 @@ !------------------------------------------------------------------------------- REAL :: alp_offset - LOGICAL found + logical found !--- Constants @@ -505,4 +505,9 @@ wake_cstar(:) = 0. wake_fip(:) = 0. + wake_pe = 0. + fm_therm = 0. + entr_therm = 0. + detr_therm = 0. + CALL fonte_neige_init(run_off_lic_0) CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil ) Index: LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F90 =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F90 (revision 1421) +++ LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F90 (revision 1425) @@ -20,4 +20,5 @@ ! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs) !------------------------------------------------------------------------------- + USE control_mod #ifdef CPP_EARTH USE dimphy @@ -27,8 +28,7 @@ NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT, & NF90_NOERR, NF90_NOWRITE, NF90_DOUBLE, NF90_GLOBAL, & - NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED + NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED, NF90_FLOAT USE inter_barxy_m, only: inter_barxy #endif - USE control_mod IMPLICIT NONE !------------------------------------------------------------------------------- @@ -267,5 +267,4 @@ #endif ! of #ifdef CPP_EARTH - STOP @@ -281,5 +280,5 @@ SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL) ! -!------------------------------------------------------------------------------- +!----------------------------------------------------------------------------- ! Comments: ! There are two assumptions concerning the NetCDF files, that are satisfied @@ -287,11 +286,15 @@ ! 1) The last dimension of the variables used is the time record. ! 2) Dimensional variables have the same names as corresponding dimensions. -!------------------------------------------------------------------------------- - USE netcdf, ONLY : NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, & - NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, & - NF90_GET_VAR, NF90_GET_ATT +!----------------------------------------------------------------------------- + USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, & + NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, & + NF90_GET_ATT USE dimphy, ONLY : klon USE phys_state_var_mod, ONLY : pctsrf USE control_mod + use pchsp_95_m, only: pchsp_95 + use pchfe_95_m, only: pchfe_95 + use arth_m, only: arth + IMPLICIT NONE #include "dimensions.h" @@ -300,5 +303,5 @@ #include "indicesol.h" #include "iniprint.h" -!------------------------------------------------------------------------------- +!----------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name @@ -306,12 +309,12 @@ LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels INTEGER, INTENT(IN) :: ndays ! current year number of days - REAL, POINTER, DIMENSION(:,:) :: champo ! output field = f(t) - LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) - REAL, OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask + REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) + LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) + REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE) -!------------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Local variables: !--- NetCDF - INTEGER :: ierr, ncid, varid ! NetCDF identifiers + INTEGER :: ncid, varid ! NetCDF identifiers CHARACTER(LEN=30) :: dnam ! dimension name CHARACTER(LEN=80) :: varname ! NetCDF variable name @@ -323,10 +326,9 @@ !--- fields INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ' - REAL, ALLOCATABLE, DIMENSION(:,:) :: champ ! wanted field on initial grid + REAL, ALLOCATABLE, DIMENSION(:, :) :: champ ! wanted field on initial grid REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear - REAL, DIMENSION(iim,jjp1) :: champint ! interpolated field - REAL, ALLOCATABLE, DIMENSION(:,:,:) :: champtime - REAL, ALLOCATABLE, DIMENSION(:,:,:) :: champan - REAL :: time, by + REAL, DIMENSION(iim, jjp1) :: champint ! interpolated field + REAL, ALLOCATABLE, DIMENSION(:, :, :) :: champtime + REAL, ALLOCATABLE, DIMENSION(:, :, :) :: champan !--- input files CHARACTER(LEN=20) :: cal_in ! calendar @@ -334,10 +336,13 @@ !--- misc INTEGER :: i, j, k, l ! loop counters - REAL, ALLOCATABLE, DIMENSION(:,:) :: work ! used for extrapolation + REAL, ALLOCATABLE, DIMENSION(:, :) :: work ! used for extrapolation CHARACTER(LEN=25) :: title ! for messages LOGICAL :: extrp ! flag for extrapolation REAL :: chmin, chmax -!------------------------------------------------------------------------------- -!---Variables depending on keyword 'mode' -------------------------------------- + INTEGER ierr + integer n_extrap ! number of extrapolated points + logical skip +!------------------------------------------------------------------------------ +!---Variables depending on keyword 'mode' ------------------------------------- NULLIFY(champo) SELECT CASE(mode) @@ -347,76 +352,80 @@ CASE('ALB'); varname='ALBEDO'; title='Albedo' END SELECT - extrp=.FALSE. - IF ( PRESENT(flag) ) THEN - IF ( flag .AND. mode=='SST' ) extrp=.TRUE. - END IF - - -!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------ - ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid); CALL ncerr(ierr,fnam) - ierr=NF90_INQ_VARID(ncid,varname,varid); CALL ncerr(ierr,fnam) - ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids); CALL ncerr(ierr,fnam) + extrp=.FALSE. + IF ( PRESENT(flag) ) THEN + IF ( flag .AND. mode=='SST' ) extrp=.TRUE. + END IF + +!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- + ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam) + ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) + ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam) !--- Longitude - ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep) - CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep)) - ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam) - ierr=NF90_GET_VAR(ncid,varid,dlon_ini); CALL ncerr(ierr,fnam) - WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep + ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep) + CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep)) + ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) + ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam) + WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep !--- Latitude - ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep) - CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep)) - ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam) - ierr=NF90_GET_VAR(ncid,varid,dlat_ini); CALL ncerr(ierr,fnam) - WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep + ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep) + CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep)) + ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) + ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam) + WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep !--- Time (variable is not needed - it is rebuilt - but calendar is) - ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep) - CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep)) - ierr=NF90_INQ_VARID(ncid,dnam,varid); CALL ncerr(ierr,fnam) + ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep) + CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep)) + ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) cal_in=' ' - ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in) + ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in) IF(ierr/=NF90_NOERR) THEN SELECT CASE(mode) - CASE('RUG','ALB'); cal_in='360d' - CASE('SIC','SST'); cal_in='gregorian' + CASE('RUG', 'ALB'); cal_in='360d' + CASE('SIC', 'SST'); cal_in='gregorian' END SELECT - WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d& - &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.' + WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & + // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' END IF - WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in - -!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ---------------------- + WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & + cal_in + +!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- !--- Determining input file number of days, depending on calendar - ndays_in=year_len(anneeref,cal_in) + ndays_in=year_len(anneeref, cal_in) !--- Time vector reconstruction (time vector from file is not trusted) !--- If input records are not monthly, time sampling has to be constant ! - timeyear=mid_months(anneeref,cal_in,lmdep) - IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode) & - //' ne comportent pas 12, mais ',lmdep,' enregistrements.' - -!--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------ - ALLOCATE(champ(imdep,jmdep),champtime(iim,jjp1,lmdep)) - IF(extrp) ALLOCATE(work(imdep,jmdep)) - - WRITE(lunout,*) - WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.' - ierr=NF90_INQ_VARID(ncid,varname,varid); CALL ncerr(ierr,fnam) - DO l=1,lmdep - ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/)) - CALL ncerr(ierr,fnam) - CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar) - - IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) + timeyear=mid_months(anneeref, cal_in, lmdep) + IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' & + // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, & + ' enregistrements.' + +!--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- + ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep)) + IF(extrp) ALLOCATE(work(imdep, jmdep)) + + WRITE(lunout, *) + WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, & + ' CHAMPS.' + ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) + DO l=1, lmdep + ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/)) + CALL ncerr(ierr, fnam) + CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, & + champ, ibar) + + IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, & + work) IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN IF(l==1) THEN - WRITE(lunout,*) & + WRITE(lunout, *) & '-------------------------------------------------------------------------' - WRITE(lunout,*) & - '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' - WRITE(lunout,*) & + WRITE(lunout, *) & + 'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' + WRITE(lunout, *) & '-------------------------------------------------------------------------' END IF @@ -434,59 +443,64 @@ CASE('SIC'); CALL sea_ice (imdep, jmdep, dlon, dlat, champ, & iim, jjp1, rlonv, rlatu, champint) - CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, & + CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, & iim, jjp1, rlonv, rlatu, champint) END SELECT END IF - champtime(:,:,l)=champint + champtime(:, :, l)=champint END DO - ierr=NF90_CLOSE(ncid); CALL ncerr(ierr,fnam) - - DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ) + ierr=NF90_CLOSE(ncid); CALL ncerr(ierr, fnam) + + DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ) IF(extrp) DEALLOCATE(work) -!--- TIME INTERPOLATION -------------------------------------------------------- - WRITE(lunout,*) - WRITE(lunout,*)'INTERPOLATION TEMPORELLE.' - WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear - WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays - ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays)) - DO j=1,jjp1 - DO i=1,iim - CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder) - DO k=1,ndays - time=FLOAT((k-1)*ndays_in)/FLOAT(ndays) - CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by) - champan(i,j,k) = by - END DO +!--- TIME INTERPOLATION ------------------------------------------------------ + WRITE(lunout, *) + WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' + WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear + WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays + ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) + skip = .false. + n_extrap = 0 + DO j=1, jjp1 + DO i=1, iim + yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, & + vc_beg=0., vc_end=0.) + CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, & + arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr) + if (ierr < 0) stop 1 + n_extrap = n_extrap + ierr END DO END DO - champan(iip1,:,:)=champan(1,:,:) - DEALLOCATE(yder,champtime,timeyear) + if (n_extrap /= 0) then + print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap + end if + champan(iip1, :, :)=champan(1, :, :) + DEALLOCATE(yder, champtime, timeyear) !--- Checking the result - DO j=1,jjp1 - CALL minmax(iip1,champan(1,j,10),chmin,chmax) - WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j + DO j=1, jjp1 + CALL minmax(iip1, champan(1, j, 10), chmin, chmax) + WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j END DO -!--- SPECIAL FILTER FOR SST: SST>271.38 ---------------------------------------- +!--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- IF(mode=='SST') THEN - WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38' + WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' WHERE(champan<271.38) champan=271.38 END IF -!--- SPECIAL FILTER FOR SIC: 0.01.0) champan=1.0 WHERE(champan<0.0) champan=0.0 END IF -!--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------ - ALLOCATE(champo(klon,ndays)) - DO k=1,ndays - CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k)) +!--- DYNAMICAL TO PHYSICAL GRID ---------------------------------------------- + ALLOCATE(champo(klon, ndays)) + DO k=1, ndays + CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k)) END DO DEALLOCATE(champan) Index: LMDZ4/trunk/libf/dyn3dpar/spline.F =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/spline.F (revision 1421) +++ (revision ) @@ -1,79 +1,0 @@ -! -! $Header$ -! - subroutine spline(x,y,n,yp1,ypn,y2) - -c - -c Routine to set up the interpolating function for a cubic spline - -c interpolation (see "Numerical Recipes" for details). - -c - implicit real (a-h,o-z) - implicit integer (i-n) - - parameter(nllm=4096) - - dimension x(n),y(n),y2(n),u(nllm) - -c -c write(6,*)(x(i),i=1,n) -c write(6,*)(y(i),i=1,n) - - if(yp1.gt.0.99E30) then -c the lower boundary condition is set - y2(1)=0. -c either to be "natural" - u(1)=0. - - else -c or else to have a specified first - y2(1)=-0.5 -c derivative - u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) - - end if - - do 11 i=2,n-1 -c decomposition loop of the tridiagonal - sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) -c algorithm. Y2 and U are used - p=sig*y2(i-1)+2. -c for temporary storage of the decompo- - y2(i)=(sig-1.)/p -c sed factors - u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) - - . /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p - - 11 continue - - if(ypn.gt.0.99E30) then -c the upper boundary condition is set - qn=0. -c either to be "natural" - un=0. - - else -c or else to have a specified first - qn=0.5 -c derivative - un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) - - end if - - y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) - - do 12 k=n-1,1,-1 -c this is the backsubstitution loop of - y2(k)=y2(k)*y2(k+1)+u(k) -c the tridiagonal algorithm - 12 continue - -c - - return - - end - Index: LMDZ4/trunk/libf/dyn3dpar/splint.F =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/splint.F (revision 1421) +++ (revision ) @@ -1,56 +1,0 @@ -! -! $Header$ -! - - SUBROUTINE splint(xa,ya,y2a,n,x,y) - -c -c Routine to compute a cubic-spline interpolated value Y given the -c value of X, the arrays XA, YA and the 2nd derivative array Y2A -c computed by SUBROUTINE SPLINE. See "Numerical Recipes" for details -c - - IMPLICIT REAL (a-h,o-z) - IMPLICIT INTEGER (i-n) - DIMENSION xa(n),ya(n),y2a(n) - - kl0=1 - - khi=n -c means of bisection - 1 IF(khi-kl0.gt.1) THEN - - k=(khi+kl0)/2 - - IF(xa(k).gt.x) THEN - - khi=k - - ELSE - - kl0=k - - END IF - - GO TO 1 - - END IF -c KL0 and KHI now bracket the X - h=xa(khi)-xa(kl0) - - IF(h.eq.0.0) STOP - a=(xa(khi)-x)/h -c evaluation of cubic spline polynomial - b=(x-xa(kl0))/h - - y=a*ya(kl0)+b*ya(khi)+((a**3-a)*y2a(kl0)+(b**3-b)*y2a(khi))*(h**2) - - ./6. - -c - - RETURN - - END - - Index: LMDZ4/trunk/libf/dyn3dpar/startvar.F90 =================================================================== --- LMDZ4/trunk/libf/dyn3dpar/startvar.F90 (revision 1421) +++ LMDZ4/trunk/libf/dyn3dpar/startvar.F90 (revision 1425) @@ -643,8 +643,10 @@ !------------------------------------------------------------------------------- ! -SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,& - lat_in2, pls_in, var3d, ibar) -! -!------------------------------------------------------------------------------- +SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, & + lon_in2, lat_in2, pls_in, var3d, ibar) + + use pchsp_95_m, only: pchsp_95 + use pchfe_95_m, only: pchfe_95 + ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: varname @@ -655,73 +657,71 @@ REAL, DIMENSION(iml), INTENT(IN) :: lon_in2 REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2 - REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls_in - REAL, DIMENSION(iml,jml,lml), INTENT(OUT) :: var3d + REAL, DIMENSION(iml, jml, lml), INTENT(IN) :: pls_in + REAL, DIMENSION(iml, jml, lml), INTENT(OUT) :: var3d LOGICAL, INTENT(IN) :: ibar -!------------------------------------------------------------------------------- +!---------------------------------------------------------------------------- ! Local variables: #include "iniprint.h" - LOGICAL :: check=.TRUE. - REAL :: bx, by, chmin, chmax - INTEGER :: ii, ij, il - REAL, DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d + LOGICAL:: check=.TRUE., skip + REAL chmin, chmax + INTEGER ii, ij, il, ierr + integer n_extrap ! number of extrapolated points + REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini - REAL, DIMENSION(:), ALLOCATABLE :: lev_dyn, ax, ay, yder - INTEGER, DIMENSION(:), ALLOCATABLE :: lind -!------------------------------------------------------------------------------- - IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D field.' - IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn - IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn - - IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn)) - CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) + REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder + +!--------------------------------------------------------------------------- + IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D field.' + IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, & + ttm_dyn + IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, & + llm_dyn + + IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) + CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, & + var_ana3d) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS - ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) - lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad - lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad + ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) + lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad + lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS - ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn)) - CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, & + ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn)) + CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, & levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar) - DEALLOCATE(lon_ini,lat_ini) + DEALLOCATE(lon_ini, lat_ini) !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro - ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) - DO il=1,llm_dyn - CALL interp_startvar(varname, ibar, il==1, & - iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2, & - lon_in, lat_in, lon_in2, lat_in2, var_tmp3d(:,:,il)) + DO il=1, llm_dyn + CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, & + lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, & + lon_in2, lat_in2, var_tmp3d(:, :, il)) END DO - DEALLOCATE(lon_rad,lat_rad) - - ALLOCATE(lind(llm_dyn)) - DO il=1,llm_dyn - lind(il) = llm_dyn-il+1 + DEALLOCATE(lon_rad, lat_rad) + +!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND + ax = lev_dyn(llm_dyn:1:-1) + skip = .false. + n_extrap = 0 + DO ij=1, jml + DO ii=1, iml-1 + ay = var_tmp3d(ii, ij, llm_dyn:1:-1) + yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.) + CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & + var3d(ii, ij, lml:1:-1), ierr) + if (ierr < 0) stop 1 + n_extrap = n_extrap + ierr + END DO END DO - -!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND - ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) - DO ij=1,jml - DO ii=1,iml-1 - ax(:)=lev_dyn(lind(:)) - ay(:)=var_tmp3d(ii,ij,lind(:)) - CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder) - DO il=1,lml - bx=pls_in(ii,ij,il) - CALL SPLINT(ax, ay, yder, llm_dyn, bx, by) - var3d(ii,ij,il) = by - END DO - END DO - var3d(iml,ij,:) = var3d(1,ij,:) + if (n_extrap /= 0) then + print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap + end if + var3d(iml, :, :) = var3d(1, :, :) + + DO il=1, lml + CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax) + WRITE(lunout, *)' '//TRIM(varname)//' min max l ', il, chmin, chmax END DO - DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind) - - DO il=1,lml - CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax) - WRITE(lunout,*)' '//TRIM(varname)//' min max l ',il,chmin,chmax - END DO - - RETURN END SUBROUTINE start_inter_3d