Index: LMDZ6/trunk/tools/make_sso/README =================================================================== --- LMDZ6/trunk/tools/make_sso/README (revision 4167) +++ LMDZ6/trunk/tools/make_sso/README (revision 4167) @@ -0,0 +1,26 @@ +make_sso + +routine to compute subfgrid scale orography parameters +from a spherical harmonics method + +Requirements: +- SPHEREPACK library +- FFTW ? +- ifortran +- netcdf + +Options: +-i input_relief_file +-res nlon nlat : of the output grid +-v : orography variable name of input file +-m : option + +Example to compile then run make_sso from a high-resolution Relief.nc file to +a 180x90 horizontal grid. + +./compile +file=Relief.nc +./make_sso.e -i $file -res 180 91 -v z + +Output file is ${file}_180x91 + Index: LMDZ6/trunk/tools/make_sso/compile =================================================================== --- LMDZ6/trunk/tools/make_sso/compile (revision 4167) +++ LMDZ6/trunk/tools/make_sso/compile (revision 4167) @@ -0,0 +1,48 @@ +case ${HOSTNAME:0:5} in + irene) + module unload netcdf-fortran + module unload netcdf-c + module load netcdf-c + module load netcdf-fortran + module load fftw3/gnu + netcdf="-L$NETCDF_LIBDIR" + ;; + cicla) + module load intel + module load netcdf4/4.3.3.1-ifort + NETCDF=$(which ncdump); NETCDF=${NETCDF%%/bin*} + NETCDFFORTRAN_INCDIR=$NETCDF/include + NETCDFFORTRAN_LIBDIR=$NETCDF/lib + netcdf="" + ;; +esac + +if [ ! -f Relief.nc ]; then +echo Attention missing Relief.nc file +echo I download it +wget https://web.lmd.jussieu.fr/~lmdz/pub/3DInputData/Limit/Relief.nc +fi + +if [ ! -d SPHEREPACK ]; then +echo Missing SPHEREPACK library. I download it +wget https://web.lmd.jussieu.fr/~lmdz/pub/import/spherepack3.2.tgz +tar -xvf spherepack3.2.tgz +mv SPHEREPACK/make.inc SPHEREPACK/make.inc.old +mv make.inc SPHEREPACK/. +rm spherepack3.2.tgz +fi + + + +if [ ! -f SPHEREPACK/lib/libspherepack.a ]; then + cd SPHEREPACK ; gmake ; cd - +fi +SRC="make_sso_SpherePack.f90" +mode=prod + +case $mode in + prod) OPTS="-O3" ;; + debug) OPTS="-check all -g" ;; +esac + +ifort -assume realloc_lhs $OPTS -o make_sso.e make_sso_SpherePack.f90 -I$NETCDFFORTRAN_INCDIR $netcdf -L$NETCDFFORTRAN_LIBDIR -lnetcdff -lnetcdf -L./SPHEREPACK/lib -lspherepack -L$FFTW3_LIBDIR -lfftw3 Index: LMDZ6/trunk/tools/make_sso/make.inc =================================================================== --- LMDZ6/trunk/tools/make_sso/make.inc (revision 4167) +++ LMDZ6/trunk/tools/make_sso/make.inc (revision 4167) @@ -0,0 +1,11 @@ + +LIB=../lib/libspherepack.a + +F90 := ifort -I../lib +CPP := ifort -E + + +MAKE := gmake +AR := /usr/bin/ar + + Index: LMDZ6/trunk/tools/make_sso/make_sso_SpherePack.f90 =================================================================== --- LMDZ6/trunk/tools/make_sso/make_sso_SpherePack.f90 (revision 4167) +++ LMDZ6/trunk/tools/make_sso/make_sso_SpherePack.f90 (revision 4167) @@ -0,0 +1,1093 @@ +!------------------------------------------------------------------------------- +! +PROGRAM make_sso +! +!------------------------------------------------------------------------------- +! Purpose: Project ETOPO file (GMT4 axes conventions) on spherical harmonics. +!------------------------------------------------------------------------------- + USE netcdf +! USE sphpack + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Local variables: + REAL, PARAMETER :: xpi=4.*ATAN(1.), & !--- pi + hp=xpi/2., d2r=xpi/180., & !--- pi, pi/2, pi/180 + REarth=2.E7/xpi, & !--- Earth radius + eps=1.E-6 !--- Tolerance for angles + !=== Transition length for spectral filters (to avoid side lobes) +! REAL, PARAMETER :: del=2000./REarth !--- On unit sphere + !=== .TRUE./.FALSE.: Fast/slow routine with heavy O(N^3)/light O(N^2) storage + LOGICAL, PARAMETER :: lfast=.TRUE. !--- Fast/slow routine + +!--- From the input arguments + CHARACTER(LEN=256), ALLOCATABLE :: args(:) !--- All the arguments + CHARACTER(LEN=2048):: call_seq !--- Calling sequence + INTEGER :: nargs !--- Input args number + CHARACTER(LEN=256) :: f_in !--- Input file name + CHARACTER(LEN=256) :: vnam !--- Height variable name + CHARACTER(LEN=256) :: res_in !--- Input resolution + CHARACTER(LEN=256) :: fmsk !--- Mask file or keyword + !=== Filter on output field, of resolution=ffactor*output resolution + LOGICAL :: lfilter=.TRUE. !--- Filter activation + REAL :: ffactor !--- Filter resol. factor + REAL :: tfactor !--- Transition factor + INTEGER :: nlon_ou, nlat_ou !--- Output grid resol. + +!--- Quantities related to input field grid + CHARACTER(LEN=256) :: lonn, latn !--- Axes vars names + CHARACTER(LEN=256) :: lonu, latu !--- Angular units + REAL, ALLOCATABLE :: lon_in(:), lat_in(:) !--- Cells centers + INTEGER :: nlon_in, nlat_in !--- Input grid dims + INTEGER :: igrid_in(2) !--- Grid identification + REAL :: dlon, dlat !--- Domain sizes + REAL :: D1, del1 !--- Resol. at equator + REAL :: al !--- Longitude ratio + LOGICAL :: lcyc !--- T: redundant longitude + LOGICAL :: lXdec, lYinc !--- T: dec/increasing lon/lat + REAL, ALLOCATABLE, DIMENSION(:,:) ::v,h , a , b,mski!--- Heights, coeffs, mask + REAL, ALLOCATABLE, DIMENSION(:,:,:) ::hs, as, bs, & !--- Same, filtered + dxZ, dyZ, dxZ2, dyZ2, dxyZ, Xk, Xl, Xm, Xp, Xq, dH!--- Grads, correlations + +!--- Quantities related to output grid + CHARACTER(LEN=256) :: f_ou !--- Output file name + CHARACTER(LEN=256) :: res_ou !--- Output resolution + INTEGER(KIND=1), ALLOCATABLE :: mask(:,:) !--- Boolean out land mask + REAL, ALLOCATABLE :: lon_ou(:), lat_ou(:), tmp1(:) !--- Output grid coords + INTEGER :: igrid_ou(2) !--- Grid identification + REAL :: D2, del2 !--- Resol. at equator + REAL, ALLOCATABLE, DIMENSION(:,:) :: msko, & !--- Fractional land mask + h0, & !--- Large scale orography + Zphi, Zsig, Zthe, Zgam, mu, tmp2 !--- Geopo, SSO parameters + +!--- Miscellanous Netcdf-related variables + INTEGER :: fID, xID, loID, phiID, sigID, theID, picID, dIDs(4), nn + INTEGER :: vID, yID, laID, meaID, gamID, mskID, valID, muID, ix180 + LOGICAL :: ll + +!--- Miscellanous SPHEREPACK-specific quantities + INTEGER :: l1, idh, mdab, lsh, la1, la2, lwa, lsvmin, lsave, ierr, intl + INTEGER :: l2, jdh, ndab, lvh, lb1, lb2, lwb, lwkmin, lwork, ldwork + REAL, ALLOCATABLE :: wsha(:), wsave(:), wvh (:) + REAL, ALLOCATABLE :: wshs(:), dwork(:), work(:) + INTEGER :: k, n, nt, isym + + CHARACTER(LEN=256) :: sub, fnam, msg, arg + +!------------------------------------------------------------------------------- + sub='make_sso' + +!=== DEAL WITH INPUT ARGUMENTS ================================================= + call_seq='./'//TRIM(sub)//'.e' + nargs=COMMAND_ARGUMENT_COUNT(); ALLOCATE(args(nargs)) +!--- MANUAL DISPLAYED IF NO ARGUMENTS ARE PROVIDED + IF(nargs==0) THEN + WRITE(*,*)'' + WRITE(*,*)'General calling sequence:' + WRITE(*,*)'' + WRITE(*,*)TRIM(call_seq)//' -i -v -res ' & + //' -m -f -t ' + WRITE(*,*)'' + WRITE(*,*)'Where:' + WRITE(*,*)'' + WRITE(*,*)' : input file name' + WRITE(*,*)' : height variable name' + WRITE(*,*)' : number of distinct longitude points of output grid' + WRITE(*,*)' : number of distinct latitude points of output grid' + WRITE(*,*)' : can be: 1. a mask file (o2a.nc)' + WRITE(*,*)' 2. "noro": computed with grid_noro' + WRITE(*,*)' default: 3. "spec": computed internally' + WRITE(*,*)' : filtering width: ffac*Do (Do:output resol.)' + WRITE(*,*)' : transition width: Di+tfac*(Do-Di) (Di: input resol.)' + WRITE(*,*) + WRITE(*,*)'Note that latitudes grid contains both poles.' + WRITE(*,*) + STOP + END IF +!--- ARGUMENTS PARSING + DO k=1,nargs + CALL GET_COMMAND_ARGUMENT(NUMBER=k,VALUE=args(k)) + call_seq=TRIM(call_seq)//' '//TRIM(args(k)) + END DO + k=1 + f_in=''; nlon_ou=0; nlat_ou=0; vnam=''; fmsk=''; ffactor=2.0; tfactor=0. + DO + k=k+1; ll=.FALSE.; arg=args(k-1) + SELECT CASE(arg) + CASE('-i'); f_in=args(k); msg='Missing file "'//TRIM(f_in)//'".' + ll=NF90_OPEN(f_in,NF90_NOWRITE,fID)/=NF90_NOERR + IF(.NOT.ll) n=NF90_CLOSE(fID) + CASE('-res'); nlon_ou=str2int(args(k)); k=k+1 + nlat_ou=str2int(args(k)) + CASE('-v'); vnam=args(k) + CASE('-m'); fmsk=args(k) + CASE('-f'); ffactor=str2real(args(k)) + CASE('-t'); tfactor=str2real(args(k)) + CASE DEFAULT; msg='Unrecognized argument "'//TRIM(arg)//'".'; ll=.TRUE. + END SELECT + IF(ll) THEN; WRITE(*,*)msg; STOP; END IF + k=k+1 + IF(k>nargs) EXIT + END DO +!--- CHECK OUTPUT GRID IS DEFINED + IF(TRIM(fmsk)=='') fmsk='spec' + IF(ALL(['noro','spec']/=fmsk)) THEN + msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file' + CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=NF90_NOERR,msg) + CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean") !--- MASK ID + CALL nc(NF90_INQUIRE_VARIABLE(fID,vID,dimids=dIDs)) !--- DIMS IDS + CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(1),len=nlon_ou),'x')!--- NB LONG + CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(2),len=nlat_ou),'y')!--- NB LAT + CALL nc(NF90_CLOSE(fID)) + END IF + IF(nlon_ou<=0.OR.nlat_ou<=0) THEN + WRITE(*,*)'Missing or wrong "-res" arguments for output grid'; STOP + END IF + +!=== READ THE INPUT FIELD ======================================================= + CALL nc(NF90_OPEN(f_in,NF90_NOWRITE,fID)) + WRITE(*,*)'>> Reading variable "'//TRIM(vnam)//'" from "'//TRIM(f_in)//'"...' + +!--- GET HEIGHT VARIABLE AND ITS DIMENSIONS + CALL nc(NF90_INQ_VARID(fID,vnam,vID)) + CALL nc(NF90_INQUIRE_VARIABLE (fID,vID,dimids=dIDs)) + CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(1),len=nlon_in,name=lonn)) + CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(2),len=nlat_in,name=latn)) + WRITE(*,*)TRIM(vnam)//' is '//TRIM(lonn)//'('//TRIM(int2str(nlon_in))//')*' & + //TRIM(latn)//'('//TRIM(int2str(nlat_in))//')' + nlon_in=nlon_in-1 !--- DISTINCT POINTS NUMBER + ALLOCATE(lon_in(nlon_in+1),lat_in(nlat_in),h(nlon_in+1,nlat_in)) + CALL nc(NF90_INQ_VARID(fID,lonn,loID) ,lonn) + CALL nc(NF90_INQ_VARID(fID,latn,laID) ,latn) + CALL nc(NF90_GET_VAR (fID,loID,lon_in) ,lonn) + CALL nc(NF90_GET_VAR (fID,laID,lat_in) ,latn) + CALL nc(NF90_GET_ATT (fID,loID,'units',lonu),lonn) + CALL nc(NF90_GET_ATT (fID,laID,'units',latu),latn) + CALL nc(NF90_GET_VAR (fID,vID,h(:,:),[1,1],[nlon_in+1,nlat_in]),vnam) + CALL nc(NF90_CLOSE(fID)) + +!--- CHECK WETHER GRID IS CORRECT (GLOBAL DOMAIN, IDENTIFIED UNITS...) + dlon=ABS(lon_in(nlon_in+1)-lon_in(1)); dlat=ABS(lat_in(nlat_in)-lat_in(1)) + CALL err(lonu(1:6)/='degree','Invalid longitudes unit: "'//TRIM(lonu)//'"') + CALL err(latu(1:6)/='degree','Invalid latitudes unit: "' //TRIM(latu)//'"') + CALL err(ABS(dlon-360.)>eps,'Longitudes domain is not global.') + CALL err(ABS(dlat-180.)>eps, 'Latitudes domain is not global.') + +!--- ENSURE LONGITUDES ARE INCREASING AND LATITUDES ARE DECREASING + lXdec=lon_in(1)>lon_in(nlon_in) + lYinc=lat_in(1)=180.) EXIT; END DO + IF(ix180>=nlon_in+1) THEN + lon_in(:)=[lon_in(ix180:nlon_in)-360.,lon_in(1:ix180)] + ALLOCATE(tmp2(nlon_in+1,nlat_in)); nn=nlon_in-ix180+1 + tmp2(1:nn,:)=h(ix180:nlon_in,:); tmp2(nn+1:nlon_in+1,:)=h(1:ix180,:) + DEALLOCATE(tmp2) + END IF + +!--- REMOVE OCEAN ; FLIP FIELD (SPHEREPACK: LAT FIRST) ; INPUT BOOLEAN MASK + WHERE(h(:,:)<0.) h(:,:)=0 + idh=nlat_in; jdh=nlon_in; CALL flip(h) + ALLOCATE(mski(idh,jdh)); mski(:,:)=0.; WHERE(h(:,:)>0.) mski(:,:)=1. + +!=== SPHERICAL HARMONICS ANALYSIS =============================================== +!--- Allocate work vector (min. length depends on used routines) + nt=2 ! MAX nt + l1=MIN(nlat_in,(nlon_in+2)/2) + IF(MODULO(nlon_in,2)==1) l1=MIN(nlat_in,(nlon_in+1)/2) + l2=(nlat_in+1)/2 + IF(lfast) THEN + lwork=2*(nlat_in+1) ! shaeci/shseci + lwork=MAX(lwork,nlat_in*(nt*nlon_in+MAX(3*l2,nlon_in))) ! shaec /shaes + lwork=MAX(lwork,4*(nlat_in+1)) ! vhseci + lwork=MAX(lwork,nlat_in*(2*nt*nlon_in+MAX(6*l2,nlon_in))+nlat_in*(2*l1*nt+1)) ! gradec + lsh=2*nlat_in*l2+3*( (l1-2 )*(nlat_in+nlat_in-l1-1))/2+nlon_in+15 ! shaec /shsec + lvh=4*nlat_in*l2+3*MAX(l1-2,0)*(nlat_in+nlat_in-l1-1) +nlon_in+15 ! vhseci/gradec + ELSE + lwork=5*nlat_in*l2+3*((l1-2)*(nlat_in+nlat_in-l1-1))/2 ! shaesi/shsesi + lwork=MAX(lwork,(nt+1)*nlat_in*nlon_in) ! shaes /shses + lwork=MAX(lwork,3*(MAX(l1-2,0)*(nlat_in+nlat_in-l1-1))/2+5*l2*nlat_in) ! vhsesi + lwork=MAX(lwork,nlat_in*((2*nt+1)*nlon_in+2*l1*nt+1)) ! grades + lsh=(l1*l2*(nlat_in+nlat_in-l1+1))/2+nlon_in+15 ! shaes /shses + lvh= l1*l2*(nlat_in+nlat_in-l1+1) +nlon_in+15 ! vhsesi/grades + END IF + ALLOCATE(work(lwork)) + ALLOCATE(wsha(lsh),wshs(lsh),wvh(lvh)) + +!--- Allocate spectral coefficients vectors + mdab=l1; ndab=nlat_in; IF(.NOT.lfast) mdab=MIN(nlat_in,(nlon_in+2)/2) + ALLOCATE(a(mdab,ndab),b(mdab,ndab)) + isym=0; nt=1 !--- Symetries, fields nb + +!--- FIELD PROJECTION ON THE SPHERICAL HARMONICS BASE. + WRITE(*,*)'>> Projecting fields on the spherical harmonics...' + IF(lfast) THEN + CALL shaeci(nlat_in,nlon_in,wsha,lsh,work,lwork,ierr) + CALL error('shaeci',ierr) + CALL shaec(nlat_in,nlon_in,isym,nt,h,idh,jdh,a,b,mdab,ndab,wsha,lsh,work,lwork,ierr) + CALL error('shaec',ierr) + ELSE + CALL shaesi(nlat_in,nlon_in,wsha,lsh,work,lwork,ierr) + CALL error('shaesi',ierr) + CALL shaes(nlat_in,nlon_in,isym,nt,h,idh,jdh,a,b,mdab,ndab,wsha,lsh,work,lwork,ierr) + CALL error('shaes',ierr) + END IF + +!--- ALIASING DAMPLING FILTER +! WRITE(*,*)'>> Bi-Laplacian filtering to avoid Gibbs phenomenon...' +! CALL diffusion_filter(a,b) + +!=== INPUT/OUTPUT RESOLUTIONS AT EQUATOR (normalized with Earth radius) +!--- The area of a cell at (lam,phi) of size(dl,dp) is: dl.SIN(dp).COS(phi) +!--- At equator, Res=SQRT(dl.SIN(dp)) ; with dl=2.pi/nlon and dp=pi/nlat: +!--- Res=SQRT(2.pi.SIN(pi/nlat)/nlon)~pi*SQRT(2/(nlon*nlat)) +!--- Note: nlat is the slices number = latitude points number -1 + D1=SQRT(2.*xpi*SIN(xpi/(nlat_in-1))/nlon_in) !--- Input grid + D2=SQRT(2.*xpi*SIN(xpi/(nlat_ou-1))/nlon_ou) !--- Output grid + +!=== FILTER FACTOR EQUAL TO ZERO => NO FILTERING + lfilter=ffactor/=0. + +!=== DEFAULT TRANSITION SCALE: INPUT GRID RESOLUTION + WRITE(*,*)'tfactor=',tfactor + del1=D1; del2=D1+(D2-D1)*tfactor + WRITE(*,*)'2*D1,del1,D2,del2=',2*D1,del1,D2,del2 + +!=== SCALES SEPARATION AT THE OUTPUT GRID RESOLUTION ============================ + nt=2 + WRITE(*,*)'>> Scales separation at '//TRIM(real2str(D2/d2r))//' degrees (deta'& + //'ils below '//TRIM(real2str(60*D1/d2r))//' minutes removed)...' + ALLOCATE(as(mdab,ndab,nt),bs(mdab,ndab,nt)) + DO n=1,nt; as(:,:,n)=a; bs(:,:,n)=b; END DO + CALL lowpass_filter(as(:,:,1),bs(:,:,1), D2,del2) !--- LARGE SCALES + CALL bandpass_filter(as(:,:,2),bs(:,:,2),2*D1,del1,D2,del2) !--- SMALL SCALES + +!--- BACKWARD TRANSFORM FOR LARGE AND SUB-CELL SCALES OROGRAPHYS + WRITE(*,*)'>> Backward transform of large and sub-cell scales orographies...' + ALLOCATE(hs(idh,jdh,nt)) + IF(lfast) THEN + CALL shseci(nlat_in,nlon_in,wshs,lsh,work,lwork,ierr) + CALL error('shseci',ierr) + CALL shsec(nlat_in,nlon_in,isym,nt,hs,idh,jdh,as,bs,mdab,ndab,wshs,lsh,work,lwork,ierr) + CALL error('shsec',ierr) + ELSE + CALL shsesi(nlat_in,nlon_in,wshs,lsh,work,lwork,ierr) + CALL error('shsesi',ierr) + CALL shses(nlat_in,nlon_in,isym,nt,hs,idh,jdh,as,bs,mdab,ndab,wshs,lsh,work,lwork,ierr) + CALL error('shses',ierr) + END IF + WHERE(mski==0.); hs(:,:,1)=0.; hs(:,:,2)=0.; END WHERE + hs(:,:,2)=hs(:,:,2)**2 + +!=== COMPUTE GRADIENT ; n=1: large scale ; n=2: sub-cell grid =================== + WRITE(*,*)'>> Computing large and sub-cell scales gradients...' + ALLOCATE(dxZ(idh,jdh,nt),dyZ(idh,jdh,nt)) + IF(lfast) THEN + CALL vhseci(nlat_in,nlon_in,wvh,lvh,work,lwork,ierr) + CALL error('vhseci',ierr) + DO n=1,2 + CALL gradec(nlat_in,nlon_in,isym,1,dyZ(:,:,n),dxZ(:,:,n),idh,jdh, & + as(:,:,n),bs(:,:,n),mdab,ndab,wvh,lvh,work,lwork,ierr) + END DO + CALL error('gradec',ierr) + ELSE + CALL vhsesi(nlat_in,nlon_in,wvh,lvh,work,lwork,ierr) + CALL error('vhsesi',ierr) + DO n=1,2 + CALL gradec(nlat_in,nlon_in,isym,1,dyZ(:,:,n),dxZ(:,:,n),idh,jdh,as(:,:,n),& + bs(:,:,n),mdab,ndab,wvh,lvh,work,lwork,ierr) + END DO + CALL error('grades',ierr) + END IF + DEALLOCATE(as,bs) + dxZ=dxZ/REarth; dyZ=dyZ/REarth + IF(lat_in(nlat_in)>lat_in(1)) dyZ=-dyZ + WHERE(mski==0.) + dxZ(:,:,1)=0.; dyZ(:,:,1)=0. + dxZ(:,:,2)=0.; dyZ(:,:,2)=0. + END WHERE + ALLOCATE(dH(idh,jdh,6)) + dH(:,:,1)=dxZ(:,:,1)*dxZ(:,:,1); dH(:,:,4)=dxZ(:,:,2)*dxZ(:,:,2) + dH(:,:,2)=dyZ(:,:,1)*dyZ(:,:,1); dH(:,:,5)=dyZ(:,:,2)*dyZ(:,:,2) + dH(:,:,3)=dxZ(:,:,1)*dyZ(:,:,1); dH(:,:,6)=dxZ(:,:,2)*dyZ(:,:,2) + DEALLOCATE(dxZ,dyZ) + +!=============================================================================== + IF(lfilter) THEN !=== FILTER TO REMOVE ALIASING =========================== +!=============================================================================== + WRITE(*,*)'>> Smoothing large and sub-cell scales fields at '// & + &TRIM(real2str(ffactor*D2/d2r))//' degrees...' + CALL spatial_filter_m(dH ,ffactor*D2,del2) + CALL spatial_filter_1(hs(:,:,2),ffactor*D2,del2) + END IF + DEALLOCATE(work) +!=============================================================================== + +!=== DEAL WITH OUTPUT GRID ===================================================== + intl=0 +!--- OUTPUT GRID (CONTAINS A LONGITUDE REDUNDANT POINT FOR NOW) + lon_ou = [(360.*DBLE(n)/DBLE(nlon_ou ), n = 0, nlon_ou)] + lon_in(1) + lat_ou = [(180.*DBLE(n)/DBLE(nlat_ou-1), n = nlat_ou-1, 0, -1)] - 90. + la1 = MIN(nlat_in,(nlon_in+1)/2); la2 = (nlat_in+1)/2 + lb1 = MIN(nlat_ou,(nlon_ou+1)/2); lb2 = (nlat_ou+1)/2 + lwa = 4*nlat_in*la2+3*MAX(la1-2,0)*(2*nlat_in-la1-1)+la2+nlon_in+15 + lwb = 4*nlat_ou*lb2+3*MAX(lb1-2,0)*(2*nlat_ou-lb1-1)+ nlon_ou+15 + lsave=lwa+lwb + l1 = MAX(nlat_in,nlat_ou) + l2 = MAX(nlon_in,nlon_ou) + lwork = 2*l1*(8*MIN(l1,(l2+2)/2)+4*l2+3) + ldwork = 2*l1*(l1+1)+1 + ALLOCATE(wsave(lsave),work(lwork),dwork(ldwork)) + +!--- Latitude grid: equally spaced partition of [-pi/2,pi/2] + igrid_in(:)=[-1,1] !--- NP->SP, nlat*nlon + igrid_ou(:)=[-1,0] !--- NP->SP, nlon*nlat + +!--- BUILD FRACTIONAL LAND MASK msko AND GEOPOTENTIAL ('noro' CASE) + SELECT CASE(fmsk) + CASE('noro') + WRITE(*,*)'>> Computing mask & geop (grid_noro compatible)...' + ALLOCATE(msko(nlon_ou+1,nlat_ou), Zphi(nlon_ou+1,nlat_ou)) + CALL flip(h) + CALL grid_noro0(lon_in*d2r, lat_in*d2r, h,lon_ou*d2r, lat_ou*d2r, Zphi, msko) + msko = msko(1:nlon_ou,:) + Zphi = Zphi(1:nlon_ou,:) + CASE('spec') + WRITE(*,*)'>> Computing mask in spectral space...' + ALLOCATE(msko(nlon_ou,nlat_ou)) + CALL trssph(intl, igrid_in, nlon_in, nlat_in, mski, & + igrid_ou, nlon_ou, nlat_ou, msko, wsave, lsave, lsvmin,& + work, lwork, lwkmin, dwork, ldwork, ierr) + intl=1 + CASE DEFAULT + WRITE(*,*)'>> Trying to read mask from "'//TRIM(fmsk)//'" file' + ALLOCATE(msko(nlon_ou,nlat_ou)) + msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file' + CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=NF90_NOERR,msg) + CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean") !--- MASK ID + CALL nc(NF90_GET_VAR(fID,vID,msko(:,:))) !--- MASK + CALL nc(NF90_CLOSE(fID)) + msko(:,:)=1.0-msko(:,:) + END SELECT + DEALLOCATE(mski,h) + +!--- GET RID OF REDUNDANT LONGITUDE (USELESS FOR SPHEREPACK) + lon_ou=lon_ou(1:nlon_ou) + +!=== REGRIDDING ON THE TARGET OUTPUT GRID ====================================== + WRITE(*,*)'>> Regridding on the target grid...' + +!--- COMPUTE BOOLEAN OUTPUT MASK + ALLOCATE(mask(nlon_ou,nlat_ou)); mask=0_1; WHERE(msko>=0.5) mask=1_1 + +!--- MEAN OROGRAPHY + ALLOCATE(h0(nlon_ou,nlat_ou)) + CALL trssph(intl, igrid_in, nlon_in, nlat_in, hs(:,:,1), & + igrid_ou, nlon_ou, nlat_ou, h0, & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + WHERE(h0<0..OR.mask==0_1) h0(:,:)=0. + +!--- STANDARD DEVIATION OF SUB-CELL SCALES OROGRAPHY + ALLOCATE(mu(nlon_ou,nlat_ou)) + CALL trssph(1, igrid_in, nlon_in, nlat_in, hs(:,:,2), & + igrid_ou, nlon_ou, nlat_ou, mu, & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + DEALLOCATE(hs) + WHERE(mu(:,:)<0..OR.mask==0_1) mu(:,:)=0.; mu=SQRT(mu) + +!--- SQUARE OF LARGE AND SUB-CELL SCALES ZONAL SLOPE + ALLOCATE(dxZ2(nlon_ou,nlat_ou,2)) + CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,1), & + igrid_ou, nlon_ou, nlat_ou, dxZ2(:,:,1), & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,4), & + igrid_ou, nlon_ou, nlat_ou, dxZ2(:,:,2), & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + WHERE(mask==0_1); dxZ2(:,:,1)=0.; dxZ2(:,:,2)=0.; END WHERE + +!--- SQUARE OF LARGE AND SUB-CELL SCALES MERIDIONAL SLOPE + ALLOCATE(dyZ2(nlon_ou,nlat_ou,2)) + CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,2), & + igrid_ou, nlon_ou, nlat_ou, dyZ2(:,:,1), & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,5), & + igrid_ou, nlon_ou, nlat_ou, dyZ2(:,:,2), & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + WHERE(mask==0_1); dyZ2(:,:,1)=0.; dyZ2(:,:,2)=0.; END WHERE + +!--- PRODUCT OF LARGE AND SUB-CELL SCALES MERIDIONAL AND ZONAL SLOPES + ALLOCATE(dxyZ(nlon_ou,nlat_ou,2)) + CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,3), & + igrid_ou, nlon_ou, nlat_ou, dxyZ(:,:,1), & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,6), & + igrid_ou, nlon_ou, nlat_ou, dxyZ(:,:,2), & + wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) + WHERE(mask==0_1); dxyZ(:,:,1)=0.; dxyZ(:,:,2)=0.; END WHERE + DEALLOCATE(dH) + +!--- COMPUTE LARGE AND SUB-CELL SCALES DERIVED QUANTITIES + ALLOCATE(Xk(nlon_ou,nlat_ou,2)); Xk=(dxZ2+dyZ2)/2. + ALLOCATE(Xl(nlon_ou,nlat_ou,2)); Xl=(dxZ2-dyZ2)/2. + DEALLOCATE(dxZ2,dyZ2) + ALLOCATE(Xm(nlon_ou,nlat_ou,2)); Xm=dxyZ + DEALLOCATE(dxyZ) + ALLOCATE(Xp(nlon_ou,nlat_ou,2),Xq(nlon_ou,nlat_ou,2)) + Xp=Xk; Xq=Xk; Xk=SQRT(Xl**2+Xm**2); Xp=Xp-Xk; Xq=Xq+Xk + DEALLOCATE(Xk) + WHERE(ABS(Xm)<=1e-16) Xm=1e-16*SIGN(1.,Xm) + ALLOCATE(Zthe(nlon_ou,nlat_ou)); Zthe=0.5*ATAN2(Xm(:,:,1),Xl(:,:,1))/d2r + DEALLOCATE(Xm,Xl) + ALLOCATE(Zgam(nlon_ou,nlat_ou)); Zgam=0. + WHERE(Xq(:,:,1)>0..AND.Xp(:,:,1)>1.E-8) Zgam(:,:)=Xp(:,:,1)/Xq(:,:,1) + DEALLOCATE(Xp) + ALLOCATE(Zsig(nlon_ou,nlat_ou)); Zsig=SQRT(Xq(:,:,2)) + WHERE(Xq(:,:,2)<=1.E-16) Zsig(:,:)=1.E-8 + DEALLOCATE(Xq) + WHERE(mask==0_1); Zthe=0.; Zgam=0.; Zsig=0.; END WHERE + +!--- Save the data + WRITE(*,*)'>> Saving fields in a netcdf file...' + res_in=TRIM(int2str(nlon_in))//'x'//TRIM(int2str(nlat_in)) + res_ou=TRIM(int2str(nlon_ou))//'x'//TRIM(int2str(nlat_ou)) + f_ou=TRIM(f_in)//'_'//TRIM(res_ou) + fnam=f_ou + CALL nc(NF90_CREATE(f_ou,NF90_CLOBBER,fID)) + + CALL nc(NF90_DEF_DIM(fID,'x',nlon_ou,xID)) + CALL nc(NF90_DEF_VAR(fID,'x',NF90_REAL,xID,loID) ,'x') + CALL nc(NF90_PUT_ATT(fID,loID,'long_name','Longitude'),'x') + CALL nc(NF90_PUT_ATT(fID,loID,'units','degrees_east') ,'x') + + CALL nc(NF90_DEF_DIM(fID,'y',nlat_ou,yID)) + CALL nc(NF90_DEF_VAR(fID,'y',NF90_REAL,yID,laID) ,'y') + CALL nc(NF90_PUT_ATT(fID,laID,'long_name','Latitude') ,'y') + CALL nc(NF90_PUT_ATT(fID,laID,'units','degrees_north'),'y') + + CALL nc(NF90_DEF_VAR(fID,'mask',NF90_REAL,[xID,yID],mskID),'mask') + IF(fmsk=='noro') & + CALL nc(NF90_DEF_VAR(fID,'Zphi',NF90_REAL,[xID,yID],phiID),'Zphi') + CALL nc(NF90_DEF_VAR(fID,'Zmea',NF90_REAL,[xID,yID],meaID),'Zmea') + CALL nc(NF90_DEF_VAR(fID,'mu' ,NF90_REAL,[xID,yID], muID),'mu' ) + CALL nc(NF90_DEF_VAR(fID,'Zsig',NF90_REAL,[xID,yID],sigID),'Zsig') + CALL nc(NF90_DEF_VAR(fID,'Zgam',NF90_REAL,[xID,yID],gamID),'Zgam') + CALL nc(NF90_DEF_VAR(fID,'Zthe',NF90_REAL,[xID,yID],theID),'Zthe') + CALL nc(NF90_DEF_VAR(fID,'Zpic',NF90_REAL,[xID,yID],picID),'Zpic') + CALL nc(NF90_DEF_VAR(fID,'Zval',NF90_REAL,[xID,yID],valID),'Zval') + + CALL nc(NF90_PUT_ATT(fID,mskID,'long_name','Fractional land mask' ),'mask') + IF(fmsk=='noro') & + CALL nc(NF90_PUT_ATT(fID,phiID,'long_name','Geopotential' ),'Zphi') + CALL nc(NF90_PUT_ATT(fID,meaID,'long_name','Mean orography' ),'Zmea') + CALL nc(NF90_PUT_ATT(fID, muID,'long_name','Std deviation of sub-cell scales orography'),'mu' ) + CALL nc(NF90_PUT_ATT(fID,sigID,'long_name','Slope along principal axis' ),'Zsig') + CALL nc(NF90_PUT_ATT(fID,gamID,'long_name','Anisotropy (aspect ratio)' ),'Zgam') + CALL nc(NF90_PUT_ATT(fID,theID,'long_name','Orientation (principal axis)' ),'Zthe') + CALL nc(NF90_PUT_ATT(fID,picID,'long_name','Maximum height' ),'Zpic') + CALL nc(NF90_PUT_ATT(fID,valID,'long_name','Minimum height' ),'Zval') + + CALL nc(NF90_PUT_ATT(fID,mskID,'units','none' ),'mask') + IF(fmsk=='noro') & + CALL nc(NF90_PUT_ATT(fID,phiID,'units','m' ),'Zphi') + CALL nc(NF90_PUT_ATT(fID,meaID,'units','m' ),'Zmea') + CALL nc(NF90_PUT_ATT(fID, muID,'units','m' ),'mu' ) + CALL nc(NF90_PUT_ATT(fID,sigID,'units','m/m' ),'Zsig') + CALL nc(NF90_PUT_ATT(fID,gamID,'units','none' ),'Zgam') + CALL nc(NF90_PUT_ATT(fID,theID,'units','degrees'),'Zthe') + CALL nc(NF90_PUT_ATT(fID,picID,'units','m' ),'Zpic') + CALL nc(NF90_PUT_ATT(fID,valID,'units','m' ),'Zval') + + CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Conventions','COARDS/CF-1.0')) + CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Initial_Grid',TRIM(res_in))) + CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'history',TRIM(call_seq))) + CALL nc(NF90_ENDDEF(fID)) + + CALL nc(NF90_PUT_VAR(fID, loID,lon_ou),'x' ) + CALL nc(NF90_PUT_VAR(fID, laID,lat_ou),'y' ) + CALL nc(NF90_PUT_VAR(fID,mskID,msko),'mask') + IF(fmsk=='noro') & + CALL nc(NF90_PUT_VAR(fID,phiID,Zphi),'Zphi') + CALL nc(NF90_PUT_VAR(fID,meaID,h0 ),'Zmea') + CALL nc(NF90_PUT_VAR(fID, muID,mu ),'mu' ) + CALL nc(NF90_PUT_VAR(fID,sigID,Zsig),'Zsig') + CALL nc(NF90_PUT_VAR(fID,gamID,Zgam),'Zgam') + CALL nc(NF90_PUT_VAR(fID,theID,Zthe),'Zthe') + CALL nc(NF90_PUT_VAR(fID,picID, h0+2*mu ),'Zpic') + CALL nc(NF90_PUT_VAR(fID,valID,MAX(0.,h0-2*mu)),'Zval') + CALL nc(NF90_CLOSE(fID)) + WRITE(*,*)'Finished.' + +CONTAINS + + +!------------------------------------------------------------------------------- +! +SUBROUTINE diffusion_filter(a,b) +! +!------------------------------------------------------------------------------- +! Purpose: ECMWF bi-laplacian diffusion filter to limit aliasing. +! Note that one third of the frequencies is skipped to allow some accuracy for +! derived quantities computation. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(INOUT) :: a(:,:), b(:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: m + REAL, ALLOCATABLE :: fn(:) +!------------------------------------------------------------------------------- + ALLOCATE(fn(ndab)); fn=[(REAL(m*(m-1)),m=1,ndab)] + a(1,:)=a(1,:)/(1+4.*fn(:)/fn(ndab))**2 + b(1,:)=b(1,:)/(1+4.*fn(:)/fn(ndab))**2 + DO m=1,ndab + IF(m<=2*ndab/3) THEN + a(m,:)=a(m,:)/(1+5.*fn(:)/fn(ndab))**2 + b(m,:)=b(m,:)/(1+5.*fn(:)/fn(ndab))**2 + ELSE + a(m,:)=0. + b(m,:)=0. + END IF + END DO + +END SUBROUTINE diffusion_filter +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE lowpass_filter(a,b,d,dd) +! +!------------------------------------------------------------------------------- +! Purpose: ECMWF lowpass filter. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(INOUT) :: a(:,:), b(:,:) + REAL, INTENT(IN) :: d, dd +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: m, na + REAL, ALLOCATABLE :: h(:) +!------------------------------------------------------------------------------- + na=SIZE(a,1) + IF(ALLOCATED(h)) THEN; IF(SIZE(h)/=na) DEALLOCATE(h); END IF + IF(.NOT.ALLOCATED(h)) ALLOCATE(h(na)) + CALL lowpass(d,dd,h) + DO m=1,SIZE(h); a(m,:)=a(m,:)*h; b(m,:)=b(m,:)*h; END DO + +END SUBROUTINE lowpass_filter +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE bandpass_filter(a,b,d1,dd1,d2,dd2) +! +!------------------------------------------------------------------------------- +! Purpose: ECMWF bandpass filter. Calls lowpass filter with 2 different scales. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(INOUT) :: a(:,:), b(:,:) + REAL, INTENT(IN) :: d1, dd1, d2, dd2 +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: m, na + REAL, ALLOCATABLE :: h(:), hh(:) +!------------------------------------------------------------------------------- + na=SIZE(a,1) + IF(ALLOCATED(hh)) THEN; IF(SIZE(hh)/=na) DEALLOCATE(hh); END IF + IF(ALLOCATED(h )) THEN; IF(SIZE(h )/=na) DEALLOCATE(h ); END IF + IF(.NOT.ALLOCATED(hh)) ALLOCATE(hh(na)) + IF(.NOT.ALLOCATED(h )) ALLOCATE(h (na)) + CALL lowpass(d1,dd1,hh) + CALL lowpass(d2,dd2,h ); hh=hh-h + DO m=1,SIZE(h); a(m,:)=a(m,:)*hh; b(m,:)=b(m,:)*hh; END DO + +END SUBROUTINE bandpass_filter +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE lowpass(d,dd,h) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(IN) :: d, dd ! d is DELTA/REarth ; dd is delta/REarth + REAL, INTENT(INOUT) :: h(:) +!------------------------------------------------------------------------------- + INTEGER :: m, n + REAL :: hd +!------------------------------------------------------------------------------- + n=SIZE(h); hd=0.5*d + h(:)=[((SIN(m*(hd-dd))/m+(COS(hp+m*hd)*SIN(hp+m*dd))/(hp/dd+m) & + + SIN(m*(hd+dd))/m+(COS(hp-m*hd)*SIN(hp-m*dd))/(hp/dd-m))**2/d/d,m=1,n)] + +END SUBROUTINE lowpass +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE spatial_filter_1(u,D,dd) +! +!------------------------------------------------------------------------------- +! Purpose: Spatial filter ; using spectral lowpass filter. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(INOUT) :: u(:,:) + REAL, INTENT(IN) :: D, dd +!------------------------------------------------------------------------------- +! Local variables: + REAL, ALLOCATABLE :: uu(:,:,:) +!------------------------------------------------------------------------------- + ALLOCATE(uu(SIZE(u,1),SIZE(u,2),1)); uu(:,:,1)=u + CALL spatial_filter_m(uu,D,dd) + u=uu(:,:,1) + DEALLOCATE(uu) + +END SUBROUTINE spatial_filter_1 +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE spatial_filter_m(u,D,dd) +! +!------------------------------------------------------------------------------- +! Purpose: Spatial filter ; using spectral lowpass filter. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(INOUT) :: u(:,:,:) + REAL, INTENT(IN) :: D, dd +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: ku + REAL, ALLOCATABLE :: aa(:,:,:),bb(:,:,:) +!------------------------------------------------------------------------------- + ku=SIZE(u,3) !--- NUMBER OF FIELDS TO SMOOTH + ALLOCATE(aa(mdab,ndab,ku),bb(mdab,ndab,ku)) + +!--- PROJECTION OF THE SUB-CELL SCALES FIELDS ON THE SPHERICAL HARMONICS BASE. + IF(lfast) THEN + CALL shaec(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wsha,lsh,work,lwork,ierr) + CALL error('shaec',ierr) + ELSE + CALL shaes(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wsha,lsh,work,lwork,ierr) + CALL error('shaes',ierr) + END IF + +!--- LOW PASS FILTER AT THE OUTPUT GRID RESOLUTION + DO n=1,ku; CALL lowpass_filter(aa(:,:,n),bb(:,:,n),D,dd); END DO + +!--- BACKWARD TRANSFORM FOR MEAN SUB-CELL SCALES PARAMETERS + IF(lfast) THEN + CALL shsec(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wshs,lsh,work,lwork,ierr) + CALL error('shsec',ierr) + ELSE + CALL shses(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wshs,lsh,work,lwork,ierr) + CALL error('shses',ierr) + END IF + DEALLOCATE(aa,bb) + +END SUBROUTINE spatial_filter_m +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE grid_noro0(xin,yin,zin,xou,you,zphi,mask) +! +!------------------------------------------------------------------------------- +! Purpose: Sub-cell scales orographic parameters coomputation. Angles in radians +! Notes: Input coordinates contain a longitude periodic redundant point. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(IN) :: xin(:), yin(:) !--- INPUT COORDINATES (nxi) (nyi) + REAL, INTENT(IN) :: zin(:,:) !--- INPUT FIELD (nxi, nyi) + REAL, INTENT(IN) :: xou(:), you(:) !--- OUTPUT COORDINATES (nxo+1) (nyo) + REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (nxo+1, nyo) + REAL, INTENT(OUT) :: mask(:,:) !--- FRACTIONAL MASK (nxo+1, nyo) +!------------------------------------------------------------------------------- +! Local variables: + REAL, ALLOCATABLE :: xi(:), xo_w(:), yo_s(:), zi(:,:), num_tot(:,:) + REAL, ALLOCATABLE :: yi(:), xo_e(:), yo_n(:), w(:,:), num_lan(:,:) + REAL :: dx, dy, wx, wy, zx, zy, ex, ey + INTEGER :: nxi, nxo, i, ii, nn + INTEGER :: nyi, nyo, j, jj, di +!------------------------------------------------------------------------------- + nxi=SIZE(xin); nxo=SIZE(xou)-1; ex=SIGN(1.,xin(nxi)-xin(1)) + nyi=SIZE(yin); nyo=SIZE(you); ey=SIGN(1.,yin(nyi)-yin(1)) + +!--- MARGIN TO ENSURE OUTPUT CELLS HAVE MATCHING INPUT CELLS + di=CEILING(1.+MAXVAL([xou(2:nxo+1)-xou(1:nxo)])/MINVAL([xin(2:nxi)-xin(1:nxi-1)])) + +!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES + ALLOCATE(xi(nxi+2*di)); xi(:)=[xin(nxi-di:nxi-1)-2.*xpi,xin,xin(2:di+1)+2.*xpi] + ALLOCATE(yi(nyi+2)); yi(:)=[2.*yin(1)-yin(2),yin,2.*yin(nyi)-yin(nyi-1)] + + ALLOCATE(zi(nxi+2*di,nyi+2)) + zi(1 +di:nxi+ di,2:nyi+1)=zin(:,:) + zi(1 : di,2:nyi+1)=zin(nxi-di:nxi-1,:) + zi(1+nxi +di:nxi+2*di,2:nyi+1)=zin(2 :di+1 ,:) + zi(1 :nxi/2+di, 1)=zi (1+nxi/2:nxi +di, 2) + zi(1+nxi/2+di:nxi+2*di, 1)=zi (1 :nxi/2+di, 2) + zi(1 :nxi/2+di, nyi+2)=zi (1+nxi/2:nxi +di,nyi+1) + zi(1+nxi/2+di:nxi+2*di, nyi+2)=zi (1 :nxi/2+di,nyi+1) + +! write(*,'(a,2f10.5,a,2f10.5)')'xin=',xin(1:2),'...',xin(nxi-1:nxi) +! write(*,'(a,2f10.5,a,2f10.5)')'yin=',yin(1:2),'...',yin(nyi-1:nyi) +! write(*,'(a,2f10.5,a,2f10.5)')'xou=',xou(1:2),'...',xou(nxo:nxo+1) +! write(*,'(a,2f10.5,a,2f10.5)')'you=',you(1:2),'...',you(nyo-1:nyo) +! write(*,*)'xi=',xi(1:di+1),'...',xi(nxi+di:nxi+2*di) +! write(*,'(a,3f10.5,a,3f10.5)')'yi=',yi(1:3),'...',yi(nyi:nyi+2) + +!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) bdac + ALLOCATE(xo_w(nxo+1)); xo_w=(xou+[2.*xou( 1)-xou( 2),xou(1:nxo )])/2.0 + ALLOCATE(yo_s(nyo )); yo_s=(you+[2.*you( 1)-you( 2),you(1:nyo-1)])/2.0 + ALLOCATE(xo_e(nxo+1)); xo_e=(xou+[xou(2:nxo+1),2.*xou(nxo+1)-xou(nxo)])/2.0 + ALLOCATE(yo_n(nyo )); yo_n=(you+[you(2:nyo ),2.*you(nyo)-you(nyo-1)])/2.0 + +!--- INITIALIZATIONS: + ALLOCATE(w(nxo+1,nyo)); w(:,:)=0.0; zphi(:,:)=0.0 + +!--- SUMMATION OVER GRIDPOINT AREA + zx=2.*xpi/REAL(nxi) !--- LON CELL LENGTH (1-SPHERE) + dx=zx/2. !--- HALF " + ALLOCATE(num_tot(nxo+1,nyo)); num_tot(:,:)=0. + ALLOCATE(num_lan(nxo+1,nyo)); num_lan(:,:)=0. + DO ii=1,nxo+1 + DO jj=1,nyo + DO j=2,nyi+1 + zy=xpi/REAL(nyi)*COS(yi(j)) !--- LAT CELL LENGTH (1-SPHERE) + dy=zy/2. + wy=MAX(0.,MIN(zy,MINVAL(dy+ey*[yo_n(jj)-yi(j),yi(j)-yo_s(jj)]))) + IF(wy==0.) CYCLE + DO i=2,nxi+2*di-1 + wx=MAX(0.,MIN(zx,MINVAL(dx+ex*[xo_e(ii)-xi(i),xi(i)-xo_w(ii)])*COS(yi(j)))) + IF(wx==0.) CYCLE + num_tot(ii,jj)=num_tot(ii,jj)+1.0 + IF(zi(i,j)>=1.) num_lan(ii,jj)=num_lan(ii,jj)+1.0 + w (ii,jj)= w(ii,jj)+wx*wy + zphi(ii,jj)=zphi(ii,jj)+wx*wy*zi(i,j) + END DO + END DO + END DO + END DO + +!--- COMPUTE FRACTIONAL MASK AND GEOPOTENTIAL + WHERE(w(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:) + nn=COUNT(w(:,:)==0.0) + IF(nn/=0) WRITE(*,*)'Problem with weight ; vanishing occurrences: ',nn + WHERE(w/=0.0) zphi=zphi/w + WHERE(w==0.0.OR.mask<0.1) zphi=0.0 + zphi(nxo+1,:)=zphi(1,:) + zphi(:, 1)=SUM(w(1:nxo, 1)*zphi(1:nxo, 1))/SUM(w(1:nxo, 1)) + zphi(:,nyo)=SUM(w(1:nxo,nyo)*zphi(1:nxo,nyo))/SUM(w(1:nxo,nyo)) + +END SUBROUTINE grid_noro0 +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE nc(ncres,var) +! +!------------------------------------------------------------------------------- +! Purpose: NetCDF errors handling. +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ncres + CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: var +!------------------------------------------------------------------------------- +! Local variables: + CHARACTER(LEN=256) :: msg +!------------------------------------------------------------------------------- + IF(ncres/=NF90_NoErr) THEN + msg='Error in routine '//TRIM(sub) + IF(fnam/='') msg=TRIM(msg)//' for file "'//TRIM(fnam)//'"' + IF(PRESENT(var)) msg=TRIM(msg)//'" and variable "'//TRIM(var)//'"' + WRITE(*,*)TRIM(msg)//': '//NF90_STRERROR(ncres); STOP + END IF + +END SUBROUTINE nc +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +ELEMENTAL FUNCTION int2str(val) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: convert an integer into a string. +!------------------------------------------------------------------------------- +! Arguments: + CHARACTER(LEN=128) :: int2str + INTEGER, INTENT(IN) :: val +!------------------------------------------------------------------------------- +! Local variable: + CHARACTER(LEN=128) :: s1 +!------------------------------------------------------------------------------- + WRITE(s1,*)val; int2str=TRIM(ADJUSTL(s1)) + +END FUNCTION int2str +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE err(ierr,str) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Argument: + LOGICAL, INTENT(IN) :: ierr + CHARACTER(LEN=*), INTENT(IN) :: str +!------------------------------------------------------------------------------- + IF(.NOT.ierr) RETURN + WRITE(*,*)TRIM(str) + STOP + +END SUBROUTINE err +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE error(sub,ierr) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + CHARACTER(LEN=*), INTENT(IN) :: sub + INTEGER, INTENT(IN) :: ierr +!------------------------------------------------------------------------------- + IF(ierr==0) RETURN + WRITE(*,*)'In routine "'//TRIM(sub)//'": error in the specification of argume& + &nt nr. '//TRIM(int2str(ierr)) + STOP +END SUBROUTINE error +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +ELEMENTAL FUNCTION real2str(val) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: convert an integer into a string. +!------------------------------------------------------------------------------- +! Arguments: + CHARACTER(LEN=128) :: real2str + REAL, INTENT(IN) :: val + INTEGER :: r, s +!------------------------------------------------------------------------------- + s=NINT(val); r=NINT(1000.*(val-REAL(s))) + IF(r<0) THEN; s= INT(val); r= INT(1000.*(val-REAL(s))); END IF + real2str=TRIM(ADJUSTL(int2str(s)))//'.'//TRIM(ADJUSTL(int2str(r))) + +END FUNCTION real2str +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE flip(a) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, ALLOCATABLE, INTENT(INOUT) :: a(:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: k, n + REAL, ALLOCATABLE :: tmp(:,:) +!------------------------------------------------------------------------------- + n=SIZE(a,2) + ALLOCATE(tmp(n,SIZE(a,1))); DO k=1,n; tmp(k,:)=a(:,k); END DO + CALL MOVE_ALLOC(FROM=tmp,TO=a) + +END SUBROUTINE flip +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +ELEMENTAL LOGICAL FUNCTION is_int(str) RESULT(out) +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: tests if str is a valid integer number +!------------------------------------------------------------------------------- + CHARACTER(LEN=*), INTENT(IN) :: str + INTEGER :: k, k0 +!------------------------------------------------------------------------------- + out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN + k0=1; IF(INDEX('+-',str(1:1))/=0) k0=2 + DO k=k0,LEN_TRIM(str); IF(INDEX('0123456789',str(k:k))==0) out=.FALSE.; END DO +END FUNCTION is_int +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +ELEMENTAL LOGICAL FUNCTION is_dec(str) RESULT(out) +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: tests if str is a valid decimal number or an integer +!------------------------------------------------------------------------------- + CHARACTER(LEN=*), INTENT(IN) :: str + INTEGER :: ix +!------------------------------------------------------------------------------- + out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN + + !--- TEST IF THE NUMBER IS A PURE INTEGER + IF(is_int(str)) RETURN + + !--- CHECK FOR POINT SEPARATOR (SHOULD HAVE A STRING INT BEFORE AND/OR AFTER + ix=INDEX(str,'.'); out=ix>0.AND.LEN_TRIM(str)/=1; IF(.NOT.out) RETURN + + !--- CHECK SURROUNDING STRINGS ARE INTEGERS + out=is_int(str(1:ix-1)); IF(ix==LEN_TRIM(str)) RETURN + out=out.AND.is_int(str(ix+1:LEN_TRIM(str))) +END FUNCTION is_dec +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +ELEMENTAL LOGICAL FUNCTION is_num(str) RESULT(out) +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: test if str is a valid integer/decimal/scientific notation number +!------------------------------------------------------------------------------- + CHARACTER(LEN=*), INTENT(IN) :: str + INTEGER :: ix, ls +!------------------------------------------------------------------------------- + out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN; ls=LEN_TRIM(str) + + !--- TEST IF THE NUMBER IS A PURE INTEGER OR A PURE DECIMAL NUMBER + IF(is_int(str)) RETURN + IF(is_dec(str)) RETURN + + !--- CHECK FOR "D" OR "E" SEPARATOR + ix=MAX(INDEX(strLow(str),'d'),INDEX(strLow(str),'e')); out=ix>0 + IF(.NOT.out) RETURN + + !--- CHECK SURROUNDING STRINGS ARE NUMBERS + out=is_dec(str(1:ix-1)).AND.is_int(str(ix+1:ls)) +END FUNCTION is_num +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +INTEGER FUNCTION str2int(str) RESULT(out) +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: convert a string into an integer +!------------------------------------------------------------------------------- + CHARACTER(LEN=*), INTENT(IN) :: str + INTEGER :: k +!------------------------------------------------------------------------------- + CALL err(.NOT.is_int(str),'Invalid integer "'//TRIM(str)//'"') + out=0 + DO k=1,LEN_TRIM(str); out=10*out+INT(INDEX('0123456789',str(k:k))-1); END DO +END FUNCTION str2int +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +REAL FUNCTION str2real(str) RESULT(out) +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: Convert a string into an real +!------------------------------------------------------------------------------- + CHARACTER(LEN=*), INTENT(IN) :: str + INTEGER :: ix, l +!------------------------------------------------------------------------------- + CALL err(.NOT.is_num(str),'Invalid number "'//TRIM(str)//'"') + ix=MAX(INDEX(strLow(str),'d'),INDEX(strLow(str),'e')); l=LEN_TRIM(str) + IF(ix==0) out=dec2real(str) + IF(ix/=0) out=dec2real(str(1:ix-1))*10.**str2int(str(ix+1:l)) +END FUNCTION str2real +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +ELEMENTAL REAL FUNCTION dec2real(str) RESULT(out) +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: convert a valid (already checked) string into an real +!------------------------------------------------------------------------------- + CHARACTER(LEN=*), INTENT(IN) :: str + INTEGER :: k, k0, ix +!------------------------------------------------------------------------------- + out=0. + k0=1; IF(INDEX('+-',str(1:1))/=0) k0=2; ix=INDEX(str,'.') !--- CHECK SIGN + DO k=k0,LEN_TRIM(str); IF(k==ix) CYCLE !--- SKIP .+- + out=out*10.+REAL(INDEX('0123456789',str(k:k))-1) + END DO + IF(ix/=0) out=out*10.**(ix-LEN_TRIM(str)) !--- COMA + IF(k0==2) out=-out !--- SIGN +END FUNCTION dec2real +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! +ELEMENTAL CHARACTER(LEN=256) FUNCTION strLow(str) RESULT(out) +! +!------------------------------------------------------------------------------- + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Purpose: Lower case conversion. +!------------------------------------------------------------------------------- +! Parameters: + CHARACTER(LEN=*), INTENT(IN) :: str +!------------------------------------------------------------------------------- +! Local variable: + INTEGER :: k, ix +!------------------------------------------------------------------------------- + out=str + DO k=1,LEN(str) + ix=IACHAR(str(k:k)); IF(64