| 1 | !------------------------------------------------------------------------------- |
|---|
| 2 | ! |
|---|
| 3 | PROGRAM make_sso |
|---|
| 4 | ! |
|---|
| 5 | !------------------------------------------------------------------------------- |
|---|
| 6 | ! Purpose: Project ETOPO file (GMT4 axes conventions) on spherical harmonics. |
|---|
| 7 | !------------------------------------------------------------------------------- |
|---|
| 8 | USE netcdf |
|---|
| 9 | ! USE sphpack |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | !------------------------------------------------------------------------------- |
|---|
| 12 | ! Local variables: |
|---|
| 13 | REAL, PARAMETER :: xpi=4.*ATAN(1.), & !--- pi |
|---|
| 14 | hp=xpi/2., d2r=xpi/180., & !--- pi, pi/2, pi/180 |
|---|
| 15 | REarth=2.E7/xpi, & !--- Earth radius |
|---|
| 16 | eps=1.E-6 !--- Tolerance for angles |
|---|
| 17 | !=== Transition length for spectral filters (to avoid side lobes) |
|---|
| 18 | ! REAL, PARAMETER :: del=2000./REarth !--- On unit sphere |
|---|
| 19 | !=== .TRUE./.FALSE.: Fast/slow routine with heavy O(N^3)/light O(N^2) storage |
|---|
| 20 | LOGICAL, PARAMETER :: lfast=.TRUE. !--- Fast/slow routine |
|---|
| 21 | |
|---|
| 22 | !--- From the input arguments |
|---|
| 23 | CHARACTER(LEN=256), ALLOCATABLE :: args(:) !--- All the arguments |
|---|
| 24 | CHARACTER(LEN=2048):: call_seq !--- Calling sequence |
|---|
| 25 | INTEGER :: nargs !--- Input args number |
|---|
| 26 | CHARACTER(LEN=256) :: f_in !--- Input file name |
|---|
| 27 | CHARACTER(LEN=256) :: vnam !--- Height variable name |
|---|
| 28 | CHARACTER(LEN=256) :: res_in !--- Input resolution |
|---|
| 29 | CHARACTER(LEN=256) :: fmsk !--- Mask file or keyword |
|---|
| 30 | !=== Filter on output field, of resolution=ffactor*output resolution |
|---|
| 31 | LOGICAL :: lfilter=.TRUE. !--- Filter activation |
|---|
| 32 | REAL :: ffactor !--- Filter resol. factor |
|---|
| 33 | REAL :: tfactor !--- Transition factor |
|---|
| 34 | INTEGER :: nlon_ou, nlat_ou !--- Output grid resol. |
|---|
| 35 | |
|---|
| 36 | !--- Quantities related to input field grid |
|---|
| 37 | CHARACTER(LEN=256) :: lonn, latn !--- Axes vars names |
|---|
| 38 | CHARACTER(LEN=256) :: lonu, latu !--- Angular units |
|---|
| 39 | REAL, ALLOCATABLE :: lon_in(:), lat_in(:) !--- Cells centers |
|---|
| 40 | INTEGER :: nlon_in, nlat_in !--- Input grid dims |
|---|
| 41 | INTEGER :: igrid_in(2) !--- Grid identification |
|---|
| 42 | REAL :: dlon, dlat !--- Domain sizes |
|---|
| 43 | REAL :: D1, del1 !--- Resol. at equator |
|---|
| 44 | REAL :: al !--- Longitude ratio |
|---|
| 45 | LOGICAL :: lcyc !--- T: redundant longitude |
|---|
| 46 | LOGICAL :: lXdec, lYinc !--- T: dec/increasing lon/lat |
|---|
| 47 | REAL, ALLOCATABLE, DIMENSION(:,:) ::v,h , a , b,mski!--- Heights, coeffs, mask |
|---|
| 48 | REAL, ALLOCATABLE, DIMENSION(:,:,:) ::hs, as, bs, & !--- Same, filtered |
|---|
| 49 | dxZ, dyZ, dxZ2, dyZ2, dxyZ, Xk, Xl, Xm, Xp, Xq, dH!--- Grads, correlations |
|---|
| 50 | |
|---|
| 51 | !--- Quantities related to output grid |
|---|
| 52 | CHARACTER(LEN=256) :: f_ou !--- Output file name |
|---|
| 53 | CHARACTER(LEN=256) :: res_ou !--- Output resolution |
|---|
| 54 | INTEGER(KIND=1), ALLOCATABLE :: mask(:,:) !--- Boolean out land mask |
|---|
| 55 | REAL, ALLOCATABLE :: lon_ou(:), lat_ou(:), tmp1(:) !--- Output grid coords |
|---|
| 56 | INTEGER :: igrid_ou(2) !--- Grid identification |
|---|
| 57 | REAL :: D2, del2 !--- Resol. at equator |
|---|
| 58 | REAL, ALLOCATABLE, DIMENSION(:,:) :: msko, & !--- Fractional land mask |
|---|
| 59 | h0, & !--- Large scale orography |
|---|
| 60 | Zphi, Zsig, Zthe, Zgam, mu, tmp2 !--- Geopo, SSO parameters |
|---|
| 61 | |
|---|
| 62 | !--- Miscellanous Netcdf-related variables |
|---|
| 63 | INTEGER :: fID, xID, loID, phiID, sigID, theID, picID, dIDs(4), nn |
|---|
| 64 | INTEGER :: vID, yID, laID, meaID, gamID, mskID, valID, muID, ix180 |
|---|
| 65 | LOGICAL :: ll |
|---|
| 66 | |
|---|
| 67 | !--- Miscellanous SPHEREPACK-specific quantities |
|---|
| 68 | INTEGER :: l1, idh, mdab, lsh, la1, la2, lwa, lsvmin, lsave, ierr, intl |
|---|
| 69 | INTEGER :: l2, jdh, ndab, lvh, lb1, lb2, lwb, lwkmin, lwork, ldwork |
|---|
| 70 | REAL, ALLOCATABLE :: wsha(:), wsave(:), wvh (:) |
|---|
| 71 | REAL, ALLOCATABLE :: wshs(:), dwork(:), work(:) |
|---|
| 72 | INTEGER :: k, n, nt, isym |
|---|
| 73 | |
|---|
| 74 | CHARACTER(LEN=256) :: sub, fnam, msg, arg |
|---|
| 75 | |
|---|
| 76 | !------------------------------------------------------------------------------- |
|---|
| 77 | sub='make_sso' |
|---|
| 78 | |
|---|
| 79 | !=== DEAL WITH INPUT ARGUMENTS ================================================= |
|---|
| 80 | call_seq='./'//TRIM(sub)//'.e' |
|---|
| 81 | nargs=COMMAND_ARGUMENT_COUNT(); ALLOCATE(args(nargs)) |
|---|
| 82 | !--- MANUAL DISPLAYED IF NO ARGUMENTS ARE PROVIDED |
|---|
| 83 | IF(nargs==0) THEN |
|---|
| 84 | WRITE(*,*)'' |
|---|
| 85 | WRITE(*,*)'General calling sequence:' |
|---|
| 86 | WRITE(*,*)'' |
|---|
| 87 | WRITE(*,*)TRIM(call_seq)//' -i <file> -v <var> -res <nlon> <nlat>' & |
|---|
| 88 | //' -m <mask> -f <ffac> -t <tfac>' |
|---|
| 89 | WRITE(*,*)'' |
|---|
| 90 | WRITE(*,*)'Where:' |
|---|
| 91 | WRITE(*,*)'' |
|---|
| 92 | WRITE(*,*)' <file>: input file name' |
|---|
| 93 | WRITE(*,*)' <var> : height variable name' |
|---|
| 94 | WRITE(*,*)' <nlon>: number of distinct longitude points of output grid' |
|---|
| 95 | WRITE(*,*)' <nlat>: number of distinct latitude points of output grid' |
|---|
| 96 | WRITE(*,*)' <mask>: can be: 1. a mask file (o2a.nc)' |
|---|
| 97 | WRITE(*,*)' 2. "noro": computed with grid_noro' |
|---|
| 98 | WRITE(*,*)' default: 3. "spec": computed internally' |
|---|
| 99 | WRITE(*,*)' <ffac>: filtering width: ffac*Do (Do:output resol.)' |
|---|
| 100 | WRITE(*,*)' <tfac>: transition width: Di+tfac*(Do-Di) (Di: input resol.)' |
|---|
| 101 | WRITE(*,*) |
|---|
| 102 | WRITE(*,*)'Note that latitudes grid contains both poles.' |
|---|
| 103 | WRITE(*,*) |
|---|
| 104 | STOP |
|---|
| 105 | END IF |
|---|
| 106 | !--- ARGUMENTS PARSING |
|---|
| 107 | DO k=1,nargs |
|---|
| 108 | CALL GET_COMMAND_ARGUMENT(NUMBER=k,VALUE=args(k)) |
|---|
| 109 | call_seq=TRIM(call_seq)//' '//TRIM(args(k)) |
|---|
| 110 | END DO |
|---|
| 111 | k=1 |
|---|
| 112 | f_in=''; nlon_ou=0; nlat_ou=0; vnam=''; fmsk=''; ffactor=2.0; tfactor=0. |
|---|
| 113 | DO |
|---|
| 114 | k=k+1; ll=.FALSE.; arg=args(k-1) |
|---|
| 115 | SELECT CASE(arg) |
|---|
| 116 | CASE('-i'); f_in=args(k); msg='Missing file "'//TRIM(f_in)//'".' |
|---|
| 117 | ll=NF90_OPEN(f_in,NF90_NOWRITE,fID)/=NF90_NOERR |
|---|
| 118 | IF(.NOT.ll) n=NF90_CLOSE(fID) |
|---|
| 119 | CASE('-res'); nlon_ou=str2int(args(k)); k=k+1 |
|---|
| 120 | nlat_ou=str2int(args(k)) |
|---|
| 121 | CASE('-v'); vnam=args(k) |
|---|
| 122 | CASE('-m'); fmsk=args(k) |
|---|
| 123 | CASE('-f'); ffactor=str2real(args(k)) |
|---|
| 124 | CASE('-t'); tfactor=str2real(args(k)) |
|---|
| 125 | CASE DEFAULT; msg='Unrecognized argument "'//TRIM(arg)//'".'; ll=.TRUE. |
|---|
| 126 | END SELECT |
|---|
| 127 | IF(ll) THEN; WRITE(*,*)msg; STOP; END IF |
|---|
| 128 | k=k+1 |
|---|
| 129 | IF(k>nargs) EXIT |
|---|
| 130 | END DO |
|---|
| 131 | !--- CHECK OUTPUT GRID IS DEFINED |
|---|
| 132 | IF(TRIM(fmsk)=='') fmsk='spec' |
|---|
| 133 | IF(ALL(['noro','spec']/=fmsk)) THEN |
|---|
| 134 | msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file' |
|---|
| 135 | CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=NF90_NOERR,msg) |
|---|
| 136 | CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean") !--- MASK ID |
|---|
| 137 | CALL nc(NF90_INQUIRE_VARIABLE(fID,vID,dimids=dIDs)) !--- DIMS IDS |
|---|
| 138 | CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(1),len=nlon_ou),'x')!--- NB LONG |
|---|
| 139 | CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(2),len=nlat_ou),'y')!--- NB LAT |
|---|
| 140 | CALL nc(NF90_CLOSE(fID)) |
|---|
| 141 | END IF |
|---|
| 142 | IF(nlon_ou<=0.OR.nlat_ou<=0) THEN |
|---|
| 143 | WRITE(*,*)'Missing or wrong "-res" arguments for output grid'; STOP |
|---|
| 144 | END IF |
|---|
| 145 | |
|---|
| 146 | !=== READ THE INPUT FIELD ======================================================= |
|---|
| 147 | CALL nc(NF90_OPEN(f_in,NF90_NOWRITE,fID)) |
|---|
| 148 | WRITE(*,*)'>> Reading variable "'//TRIM(vnam)//'" from "'//TRIM(f_in)//'"...' |
|---|
| 149 | |
|---|
| 150 | !--- GET HEIGHT VARIABLE AND ITS DIMENSIONS |
|---|
| 151 | CALL nc(NF90_INQ_VARID(fID,vnam,vID)) |
|---|
| 152 | CALL nc(NF90_INQUIRE_VARIABLE (fID,vID,dimids=dIDs)) |
|---|
| 153 | CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(1),len=nlon_in,name=lonn)) |
|---|
| 154 | CALL nc(NF90_INQUIRE_DIMENSION(fID,dIDs(2),len=nlat_in,name=latn)) |
|---|
| 155 | WRITE(*,*)TRIM(vnam)//' is '//TRIM(lonn)//'('//TRIM(int2str(nlon_in))//')*' & |
|---|
| 156 | //TRIM(latn)//'('//TRIM(int2str(nlat_in))//')' |
|---|
| 157 | nlon_in=nlon_in-1 !--- DISTINCT POINTS NUMBER |
|---|
| 158 | ALLOCATE(lon_in(nlon_in+1),lat_in(nlat_in),h(nlon_in+1,nlat_in)) |
|---|
| 159 | CALL nc(NF90_INQ_VARID(fID,lonn,loID) ,lonn) |
|---|
| 160 | CALL nc(NF90_INQ_VARID(fID,latn,laID) ,latn) |
|---|
| 161 | CALL nc(NF90_GET_VAR (fID,loID,lon_in) ,lonn) |
|---|
| 162 | CALL nc(NF90_GET_VAR (fID,laID,lat_in) ,latn) |
|---|
| 163 | CALL nc(NF90_GET_ATT (fID,loID,'units',lonu),lonn) |
|---|
| 164 | CALL nc(NF90_GET_ATT (fID,laID,'units',latu),latn) |
|---|
| 165 | CALL nc(NF90_GET_VAR (fID,vID,h(:,:),[1,1],[nlon_in+1,nlat_in]),vnam) |
|---|
| 166 | CALL nc(NF90_CLOSE(fID)) |
|---|
| 167 | |
|---|
| 168 | !--- CHECK WETHER GRID IS CORRECT (GLOBAL DOMAIN, IDENTIFIED UNITS...) |
|---|
| 169 | dlon=ABS(lon_in(nlon_in+1)-lon_in(1)); dlat=ABS(lat_in(nlat_in)-lat_in(1)) |
|---|
| 170 | CALL err(lonu(1:6)/='degree','Invalid longitudes unit: "'//TRIM(lonu)//'"') |
|---|
| 171 | CALL err(latu(1:6)/='degree','Invalid latitudes unit: "' //TRIM(latu)//'"') |
|---|
| 172 | CALL err(ABS(dlon-360.)>eps,'Longitudes domain is not global.') |
|---|
| 173 | CALL err(ABS(dlat-180.)>eps, 'Latitudes domain is not global.') |
|---|
| 174 | |
|---|
| 175 | !--- ENSURE LONGITUDES ARE INCREASING AND LATITUDES ARE DECREASING |
|---|
| 176 | lXdec=lon_in(1)>lon_in(nlon_in) |
|---|
| 177 | lYinc=lat_in(1)<lat_in(nlat_in) |
|---|
| 178 | IF(lXdec) THEN; lon_in=lon_in(nlon_in+1:1:-1); h(:,:)=h(nlon_in+1:1:-1,:); END IF |
|---|
| 179 | IF(lYinc) THEN; lat_in=lat_in(nlat_in :1:-1); h(:,:)=h(:,nlat_in :1:-1); END IF |
|---|
| 180 | |
|---|
| 181 | !--- ENSURE LONGITUDES ARE BETWEEN -180 and 180deg WITHOUT REDUNDANT END POINT |
|---|
| 182 | DO ix180=1,nlon_in+1; IF(lon_in(ix180)>=180.) EXIT; END DO |
|---|
| 183 | IF(ix180>=nlon_in+1) THEN |
|---|
| 184 | lon_in(:)=[lon_in(ix180:nlon_in)-360.,lon_in(1:ix180)] |
|---|
| 185 | ALLOCATE(tmp2(nlon_in+1,nlat_in)); nn=nlon_in-ix180+1 |
|---|
| 186 | tmp2(1:nn,:)=h(ix180:nlon_in,:); tmp2(nn+1:nlon_in+1,:)=h(1:ix180,:) |
|---|
| 187 | DEALLOCATE(tmp2) |
|---|
| 188 | END IF |
|---|
| 189 | |
|---|
| 190 | !--- REMOVE OCEAN ; FLIP FIELD (SPHEREPACK: LAT FIRST) ; INPUT BOOLEAN MASK |
|---|
| 191 | WHERE(h(:,:)<0.) h(:,:)=0 |
|---|
| 192 | idh=nlat_in; jdh=nlon_in; CALL flip(h) |
|---|
| 193 | ALLOCATE(mski(idh,jdh)); mski(:,:)=0.; WHERE(h(:,:)>0.) mski(:,:)=1. |
|---|
| 194 | |
|---|
| 195 | !=== SPHERICAL HARMONICS ANALYSIS =============================================== |
|---|
| 196 | !--- Allocate work vector (min. length depends on used routines) |
|---|
| 197 | nt=2 ! MAX nt |
|---|
| 198 | l1=MIN(nlat_in,(nlon_in+2)/2) |
|---|
| 199 | IF(MODULO(nlon_in,2)==1) l1=MIN(nlat_in,(nlon_in+1)/2) |
|---|
| 200 | l2=(nlat_in+1)/2 |
|---|
| 201 | IF(lfast) THEN |
|---|
| 202 | lwork=2*(nlat_in+1) ! shaeci/shseci |
|---|
| 203 | lwork=MAX(lwork,nlat_in*(nt*nlon_in+MAX(3*l2,nlon_in))) ! shaec /shaes |
|---|
| 204 | lwork=MAX(lwork,4*(nlat_in+1)) ! vhseci |
|---|
| 205 | lwork=MAX(lwork,nlat_in*(2*nt*nlon_in+MAX(6*l2,nlon_in))+nlat_in*(2*l1*nt+1)) ! gradec |
|---|
| 206 | lsh=2*nlat_in*l2+3*( (l1-2 )*(nlat_in+nlat_in-l1-1))/2+nlon_in+15 ! shaec /shsec |
|---|
| 207 | lvh=4*nlat_in*l2+3*MAX(l1-2,0)*(nlat_in+nlat_in-l1-1) +nlon_in+15 ! vhseci/gradec |
|---|
| 208 | ELSE |
|---|
| 209 | lwork=5*nlat_in*l2+3*((l1-2)*(nlat_in+nlat_in-l1-1))/2 ! shaesi/shsesi |
|---|
| 210 | lwork=MAX(lwork,(nt+1)*nlat_in*nlon_in) ! shaes /shses |
|---|
| 211 | lwork=MAX(lwork,3*(MAX(l1-2,0)*(nlat_in+nlat_in-l1-1))/2+5*l2*nlat_in) ! vhsesi |
|---|
| 212 | lwork=MAX(lwork,nlat_in*((2*nt+1)*nlon_in+2*l1*nt+1)) ! grades |
|---|
| 213 | lsh=(l1*l2*(nlat_in+nlat_in-l1+1))/2+nlon_in+15 ! shaes /shses |
|---|
| 214 | lvh= l1*l2*(nlat_in+nlat_in-l1+1) +nlon_in+15 ! vhsesi/grades |
|---|
| 215 | END IF |
|---|
| 216 | ALLOCATE(work(lwork)) |
|---|
| 217 | ALLOCATE(wsha(lsh),wshs(lsh),wvh(lvh)) |
|---|
| 218 | |
|---|
| 219 | !--- Allocate spectral coefficients vectors |
|---|
| 220 | mdab=l1; ndab=nlat_in; IF(.NOT.lfast) mdab=MIN(nlat_in,(nlon_in+2)/2) |
|---|
| 221 | ALLOCATE(a(mdab,ndab),b(mdab,ndab)) |
|---|
| 222 | isym=0; nt=1 !--- Symetries, fields nb |
|---|
| 223 | |
|---|
| 224 | !--- FIELD PROJECTION ON THE SPHERICAL HARMONICS BASE. |
|---|
| 225 | WRITE(*,*)'>> Projecting fields on the spherical harmonics...' |
|---|
| 226 | IF(lfast) THEN |
|---|
| 227 | CALL shaeci(nlat_in,nlon_in,wsha,lsh,work,lwork,ierr) |
|---|
| 228 | CALL error('shaeci',ierr) |
|---|
| 229 | CALL shaec(nlat_in,nlon_in,isym,nt,h,idh,jdh,a,b,mdab,ndab,wsha,lsh,work,lwork,ierr) |
|---|
| 230 | CALL error('shaec',ierr) |
|---|
| 231 | ELSE |
|---|
| 232 | CALL shaesi(nlat_in,nlon_in,wsha,lsh,work,lwork,ierr) |
|---|
| 233 | CALL error('shaesi',ierr) |
|---|
| 234 | CALL shaes(nlat_in,nlon_in,isym,nt,h,idh,jdh,a,b,mdab,ndab,wsha,lsh,work,lwork,ierr) |
|---|
| 235 | CALL error('shaes',ierr) |
|---|
| 236 | END IF |
|---|
| 237 | |
|---|
| 238 | !--- ALIASING DAMPLING FILTER |
|---|
| 239 | ! WRITE(*,*)'>> Bi-Laplacian filtering to avoid Gibbs phenomenon...' |
|---|
| 240 | ! CALL diffusion_filter(a,b) |
|---|
| 241 | |
|---|
| 242 | !=== INPUT/OUTPUT RESOLUTIONS AT EQUATOR (normalized with Earth radius) |
|---|
| 243 | !--- The area of a cell at (lam,phi) of size(dl,dp) is: dl.SIN(dp).COS(phi) |
|---|
| 244 | !--- At equator, Res=SQRT(dl.SIN(dp)) ; with dl=2.pi/nlon and dp=pi/nlat: |
|---|
| 245 | !--- Res=SQRT(2.pi.SIN(pi/nlat)/nlon)~pi*SQRT(2/(nlon*nlat)) |
|---|
| 246 | !--- Note: nlat is the slices number = latitude points number -1 |
|---|
| 247 | D1=SQRT(2.*xpi*SIN(xpi/(nlat_in-1))/nlon_in) !--- Input grid |
|---|
| 248 | D2=SQRT(2.*xpi*SIN(xpi/(nlat_ou-1))/nlon_ou) !--- Output grid |
|---|
| 249 | |
|---|
| 250 | !=== FILTER FACTOR EQUAL TO ZERO => NO FILTERING |
|---|
| 251 | lfilter=ffactor/=0. |
|---|
| 252 | |
|---|
| 253 | !=== DEFAULT TRANSITION SCALE: INPUT GRID RESOLUTION |
|---|
| 254 | WRITE(*,*)'tfactor=',tfactor |
|---|
| 255 | del1=D1; del2=D1+(D2-D1)*tfactor |
|---|
| 256 | WRITE(*,*)'2*D1,del1,D2,del2=',2*D1,del1,D2,del2 |
|---|
| 257 | |
|---|
| 258 | !=== SCALES SEPARATION AT THE OUTPUT GRID RESOLUTION ============================ |
|---|
| 259 | nt=2 |
|---|
| 260 | WRITE(*,*)'>> Scales separation at '//TRIM(real2str(D2/d2r))//' degrees (deta'& |
|---|
| 261 | //'ils below '//TRIM(real2str(60*D1/d2r))//' minutes removed)...' |
|---|
| 262 | ALLOCATE(as(mdab,ndab,nt),bs(mdab,ndab,nt)) |
|---|
| 263 | DO n=1,nt; as(:,:,n)=a; bs(:,:,n)=b; END DO |
|---|
| 264 | CALL lowpass_filter(as(:,:,1),bs(:,:,1), D2,del2) !--- LARGE SCALES |
|---|
| 265 | CALL bandpass_filter(as(:,:,2),bs(:,:,2),2*D1,del1,D2,del2) !--- SMALL SCALES |
|---|
| 266 | |
|---|
| 267 | !--- BACKWARD TRANSFORM FOR LARGE AND SUB-CELL SCALES OROGRAPHYS |
|---|
| 268 | WRITE(*,*)'>> Backward transform of large and sub-cell scales orographies...' |
|---|
| 269 | ALLOCATE(hs(idh,jdh,nt)) |
|---|
| 270 | IF(lfast) THEN |
|---|
| 271 | CALL shseci(nlat_in,nlon_in,wshs,lsh,work,lwork,ierr) |
|---|
| 272 | CALL error('shseci',ierr) |
|---|
| 273 | CALL shsec(nlat_in,nlon_in,isym,nt,hs,idh,jdh,as,bs,mdab,ndab,wshs,lsh,work,lwork,ierr) |
|---|
| 274 | CALL error('shsec',ierr) |
|---|
| 275 | ELSE |
|---|
| 276 | CALL shsesi(nlat_in,nlon_in,wshs,lsh,work,lwork,ierr) |
|---|
| 277 | CALL error('shsesi',ierr) |
|---|
| 278 | CALL shses(nlat_in,nlon_in,isym,nt,hs,idh,jdh,as,bs,mdab,ndab,wshs,lsh,work,lwork,ierr) |
|---|
| 279 | CALL error('shses',ierr) |
|---|
| 280 | END IF |
|---|
| 281 | WHERE(mski==0.); hs(:,:,1)=0.; hs(:,:,2)=0.; END WHERE |
|---|
| 282 | hs(:,:,2)=hs(:,:,2)**2 |
|---|
| 283 | |
|---|
| 284 | !=== COMPUTE GRADIENT ; n=1: large scale ; n=2: sub-cell grid =================== |
|---|
| 285 | WRITE(*,*)'>> Computing large and sub-cell scales gradients...' |
|---|
| 286 | ALLOCATE(dxZ(idh,jdh,nt),dyZ(idh,jdh,nt)) |
|---|
| 287 | IF(lfast) THEN |
|---|
| 288 | CALL vhseci(nlat_in,nlon_in,wvh,lvh,work,lwork,ierr) |
|---|
| 289 | CALL error('vhseci',ierr) |
|---|
| 290 | DO n=1,2 |
|---|
| 291 | CALL gradec(nlat_in,nlon_in,isym,1,dyZ(:,:,n),dxZ(:,:,n),idh,jdh, & |
|---|
| 292 | as(:,:,n),bs(:,:,n),mdab,ndab,wvh,lvh,work,lwork,ierr) |
|---|
| 293 | END DO |
|---|
| 294 | CALL error('gradec',ierr) |
|---|
| 295 | ELSE |
|---|
| 296 | CALL vhsesi(nlat_in,nlon_in,wvh,lvh,work,lwork,ierr) |
|---|
| 297 | CALL error('vhsesi',ierr) |
|---|
| 298 | DO n=1,2 |
|---|
| 299 | CALL gradec(nlat_in,nlon_in,isym,1,dyZ(:,:,n),dxZ(:,:,n),idh,jdh,as(:,:,n),& |
|---|
| 300 | bs(:,:,n),mdab,ndab,wvh,lvh,work,lwork,ierr) |
|---|
| 301 | END DO |
|---|
| 302 | CALL error('grades',ierr) |
|---|
| 303 | END IF |
|---|
| 304 | DEALLOCATE(as,bs) |
|---|
| 305 | dxZ=dxZ/REarth; dyZ=dyZ/REarth |
|---|
| 306 | IF(lat_in(nlat_in)>lat_in(1)) dyZ=-dyZ |
|---|
| 307 | WHERE(mski==0.) |
|---|
| 308 | dxZ(:,:,1)=0.; dyZ(:,:,1)=0. |
|---|
| 309 | dxZ(:,:,2)=0.; dyZ(:,:,2)=0. |
|---|
| 310 | END WHERE |
|---|
| 311 | ALLOCATE(dH(idh,jdh,6)) |
|---|
| 312 | dH(:,:,1)=dxZ(:,:,1)*dxZ(:,:,1); dH(:,:,4)=dxZ(:,:,2)*dxZ(:,:,2) |
|---|
| 313 | dH(:,:,2)=dyZ(:,:,1)*dyZ(:,:,1); dH(:,:,5)=dyZ(:,:,2)*dyZ(:,:,2) |
|---|
| 314 | dH(:,:,3)=dxZ(:,:,1)*dyZ(:,:,1); dH(:,:,6)=dxZ(:,:,2)*dyZ(:,:,2) |
|---|
| 315 | DEALLOCATE(dxZ,dyZ) |
|---|
| 316 | |
|---|
| 317 | !=============================================================================== |
|---|
| 318 | IF(lfilter) THEN !=== FILTER TO REMOVE ALIASING =========================== |
|---|
| 319 | !=============================================================================== |
|---|
| 320 | WRITE(*,*)'>> Smoothing large and sub-cell scales fields at '// & |
|---|
| 321 | &TRIM(real2str(ffactor*D2/d2r))//' degrees...' |
|---|
| 322 | CALL spatial_filter_m(dH ,ffactor*D2,del2) |
|---|
| 323 | CALL spatial_filter_1(hs(:,:,2),ffactor*D2,del2) |
|---|
| 324 | END IF |
|---|
| 325 | DEALLOCATE(work) |
|---|
| 326 | !=============================================================================== |
|---|
| 327 | |
|---|
| 328 | !=== DEAL WITH OUTPUT GRID ===================================================== |
|---|
| 329 | intl=0 |
|---|
| 330 | !--- OUTPUT GRID (CONTAINS A LONGITUDE REDUNDANT POINT FOR NOW) |
|---|
| 331 | lon_ou = [(360.*DBLE(n)/DBLE(nlon_ou ), n = 0, nlon_ou)] + lon_in(1) |
|---|
| 332 | lat_ou = [(180.*DBLE(n)/DBLE(nlat_ou-1), n = nlat_ou-1, 0, -1)] - 90. |
|---|
| 333 | la1 = MIN(nlat_in,(nlon_in+1)/2); la2 = (nlat_in+1)/2 |
|---|
| 334 | lb1 = MIN(nlat_ou,(nlon_ou+1)/2); lb2 = (nlat_ou+1)/2 |
|---|
| 335 | lwa = 4*nlat_in*la2+3*MAX(la1-2,0)*(2*nlat_in-la1-1)+la2+nlon_in+15 |
|---|
| 336 | lwb = 4*nlat_ou*lb2+3*MAX(lb1-2,0)*(2*nlat_ou-lb1-1)+ nlon_ou+15 |
|---|
| 337 | lsave=lwa+lwb |
|---|
| 338 | l1 = MAX(nlat_in,nlat_ou) |
|---|
| 339 | l2 = MAX(nlon_in,nlon_ou) |
|---|
| 340 | lwork = 2*l1*(8*MIN(l1,(l2+2)/2)+4*l2+3) |
|---|
| 341 | ldwork = 2*l1*(l1+1)+1 |
|---|
| 342 | ALLOCATE(wsave(lsave),work(lwork),dwork(ldwork)) |
|---|
| 343 | |
|---|
| 344 | !--- Latitude grid: equally spaced partition of [-pi/2,pi/2] |
|---|
| 345 | igrid_in(:)=[-1,1] !--- NP->SP, nlat*nlon |
|---|
| 346 | igrid_ou(:)=[-1,0] !--- NP->SP, nlon*nlat |
|---|
| 347 | |
|---|
| 348 | !--- BUILD FRACTIONAL LAND MASK msko AND GEOPOTENTIAL ('noro' CASE) |
|---|
| 349 | SELECT CASE(fmsk) |
|---|
| 350 | CASE('noro') |
|---|
| 351 | WRITE(*,*)'>> Computing mask & geop (grid_noro compatible)...' |
|---|
| 352 | ALLOCATE(msko(nlon_ou+1,nlat_ou), Zphi(nlon_ou+1,nlat_ou)) |
|---|
| 353 | CALL flip(h) |
|---|
| 354 | CALL grid_noro0(lon_in*d2r, lat_in*d2r, h,lon_ou*d2r, lat_ou*d2r, Zphi, msko) |
|---|
| 355 | msko = msko(1:nlon_ou,:) |
|---|
| 356 | Zphi = Zphi(1:nlon_ou,:) |
|---|
| 357 | CASE('spec') |
|---|
| 358 | WRITE(*,*)'>> Computing mask in spectral space...' |
|---|
| 359 | ALLOCATE(msko(nlon_ou,nlat_ou)) |
|---|
| 360 | CALL trssph(intl, igrid_in, nlon_in, nlat_in, mski, & |
|---|
| 361 | igrid_ou, nlon_ou, nlat_ou, msko, wsave, lsave, lsvmin,& |
|---|
| 362 | work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 363 | intl=1 |
|---|
| 364 | CASE DEFAULT |
|---|
| 365 | WRITE(*,*)'>> Trying to read mask from "'//TRIM(fmsk)//'" file' |
|---|
| 366 | ALLOCATE(msko(nlon_ou,nlat_ou)) |
|---|
| 367 | msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file' |
|---|
| 368 | CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=NF90_NOERR,msg) |
|---|
| 369 | CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean") !--- MASK ID |
|---|
| 370 | CALL nc(NF90_GET_VAR(fID,vID,msko(:,:))) !--- MASK |
|---|
| 371 | CALL nc(NF90_CLOSE(fID)) |
|---|
| 372 | msko(:,:)=1.0-msko(:,:) |
|---|
| 373 | END SELECT |
|---|
| 374 | DEALLOCATE(mski,h) |
|---|
| 375 | |
|---|
| 376 | !--- GET RID OF REDUNDANT LONGITUDE (USELESS FOR SPHEREPACK) |
|---|
| 377 | lon_ou=lon_ou(1:nlon_ou) |
|---|
| 378 | |
|---|
| 379 | !=== REGRIDDING ON THE TARGET OUTPUT GRID ====================================== |
|---|
| 380 | WRITE(*,*)'>> Regridding on the target grid...' |
|---|
| 381 | |
|---|
| 382 | !--- COMPUTE BOOLEAN OUTPUT MASK |
|---|
| 383 | ALLOCATE(mask(nlon_ou,nlat_ou)); mask=0_1; WHERE(msko>=0.5) mask=1_1 |
|---|
| 384 | |
|---|
| 385 | !--- MEAN OROGRAPHY |
|---|
| 386 | ALLOCATE(h0(nlon_ou,nlat_ou)) |
|---|
| 387 | CALL trssph(intl, igrid_in, nlon_in, nlat_in, hs(:,:,1), & |
|---|
| 388 | igrid_ou, nlon_ou, nlat_ou, h0, & |
|---|
| 389 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 390 | WHERE(h0<0..OR.mask==0_1) h0(:,:)=0. |
|---|
| 391 | |
|---|
| 392 | !--- STANDARD DEVIATION OF SUB-CELL SCALES OROGRAPHY |
|---|
| 393 | ALLOCATE(mu(nlon_ou,nlat_ou)) |
|---|
| 394 | CALL trssph(1, igrid_in, nlon_in, nlat_in, hs(:,:,2), & |
|---|
| 395 | igrid_ou, nlon_ou, nlat_ou, mu, & |
|---|
| 396 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 397 | DEALLOCATE(hs) |
|---|
| 398 | WHERE(mu(:,:)<0..OR.mask==0_1) mu(:,:)=0.; mu=SQRT(mu) |
|---|
| 399 | |
|---|
| 400 | !--- SQUARE OF LARGE AND SUB-CELL SCALES ZONAL SLOPE |
|---|
| 401 | ALLOCATE(dxZ2(nlon_ou,nlat_ou,2)) |
|---|
| 402 | CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,1), & |
|---|
| 403 | igrid_ou, nlon_ou, nlat_ou, dxZ2(:,:,1), & |
|---|
| 404 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 405 | CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,4), & |
|---|
| 406 | igrid_ou, nlon_ou, nlat_ou, dxZ2(:,:,2), & |
|---|
| 407 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 408 | WHERE(mask==0_1); dxZ2(:,:,1)=0.; dxZ2(:,:,2)=0.; END WHERE |
|---|
| 409 | |
|---|
| 410 | !--- SQUARE OF LARGE AND SUB-CELL SCALES MERIDIONAL SLOPE |
|---|
| 411 | ALLOCATE(dyZ2(nlon_ou,nlat_ou,2)) |
|---|
| 412 | CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,2), & |
|---|
| 413 | igrid_ou, nlon_ou, nlat_ou, dyZ2(:,:,1), & |
|---|
| 414 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 415 | CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,5), & |
|---|
| 416 | igrid_ou, nlon_ou, nlat_ou, dyZ2(:,:,2), & |
|---|
| 417 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 418 | WHERE(mask==0_1); dyZ2(:,:,1)=0.; dyZ2(:,:,2)=0.; END WHERE |
|---|
| 419 | |
|---|
| 420 | !--- PRODUCT OF LARGE AND SUB-CELL SCALES MERIDIONAL AND ZONAL SLOPES |
|---|
| 421 | ALLOCATE(dxyZ(nlon_ou,nlat_ou,2)) |
|---|
| 422 | CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,3), & |
|---|
| 423 | igrid_ou, nlon_ou, nlat_ou, dxyZ(:,:,1), & |
|---|
| 424 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 425 | CALL trssph(1, igrid_in, nlon_in, nlat_in, dH (:,:,6), & |
|---|
| 426 | igrid_ou, nlon_ou, nlat_ou, dxyZ(:,:,2), & |
|---|
| 427 | wsave, lsave, lsvmin, work, lwork, lwkmin, dwork, ldwork, ierr) |
|---|
| 428 | WHERE(mask==0_1); dxyZ(:,:,1)=0.; dxyZ(:,:,2)=0.; END WHERE |
|---|
| 429 | DEALLOCATE(dH) |
|---|
| 430 | |
|---|
| 431 | !--- COMPUTE LARGE AND SUB-CELL SCALES DERIVED QUANTITIES |
|---|
| 432 | ALLOCATE(Xk(nlon_ou,nlat_ou,2)); Xk=(dxZ2+dyZ2)/2. |
|---|
| 433 | ALLOCATE(Xl(nlon_ou,nlat_ou,2)); Xl=(dxZ2-dyZ2)/2. |
|---|
| 434 | DEALLOCATE(dxZ2,dyZ2) |
|---|
| 435 | ALLOCATE(Xm(nlon_ou,nlat_ou,2)); Xm=dxyZ |
|---|
| 436 | DEALLOCATE(dxyZ) |
|---|
| 437 | ALLOCATE(Xp(nlon_ou,nlat_ou,2),Xq(nlon_ou,nlat_ou,2)) |
|---|
| 438 | Xp=Xk; Xq=Xk; Xk=SQRT(Xl**2+Xm**2); Xp=Xp-Xk; Xq=Xq+Xk |
|---|
| 439 | DEALLOCATE(Xk) |
|---|
| 440 | WHERE(ABS(Xm)<=1e-16) Xm=1e-16*SIGN(1.,Xm) |
|---|
| 441 | ALLOCATE(Zthe(nlon_ou,nlat_ou)); Zthe=0.5*ATAN2(Xm(:,:,1),Xl(:,:,1))/d2r |
|---|
| 442 | DEALLOCATE(Xm,Xl) |
|---|
| 443 | ALLOCATE(Zgam(nlon_ou,nlat_ou)); Zgam=0. |
|---|
| 444 | WHERE(Xq(:,:,1)>0..AND.Xp(:,:,1)>1.E-8) Zgam(:,:)=Xp(:,:,1)/Xq(:,:,1) |
|---|
| 445 | DEALLOCATE(Xp) |
|---|
| 446 | ALLOCATE(Zsig(nlon_ou,nlat_ou)); Zsig=SQRT(Xq(:,:,2)) |
|---|
| 447 | WHERE(Xq(:,:,2)<=1.E-16) Zsig(:,:)=1.E-8 |
|---|
| 448 | DEALLOCATE(Xq) |
|---|
| 449 | WHERE(mask==0_1); Zthe=0.; Zgam=0.; Zsig=0.; END WHERE |
|---|
| 450 | |
|---|
| 451 | !--- Save the data |
|---|
| 452 | WRITE(*,*)'>> Saving fields in a netcdf file...' |
|---|
| 453 | res_in=TRIM(int2str(nlon_in))//'x'//TRIM(int2str(nlat_in)) |
|---|
| 454 | res_ou=TRIM(int2str(nlon_ou))//'x'//TRIM(int2str(nlat_ou)) |
|---|
| 455 | f_ou='make_sso_'//TRIM(res_ou)//'_'//TRIM(f_in) |
|---|
| 456 | fnam=f_ou |
|---|
| 457 | CALL nc(NF90_CREATE(f_ou,NF90_CLOBBER,fID)) |
|---|
| 458 | |
|---|
| 459 | CALL nc(NF90_DEF_DIM(fID,'x',nlon_ou,xID)) |
|---|
| 460 | CALL nc(NF90_DEF_VAR(fID,'x',NF90_REAL,xID,loID) ,'x') |
|---|
| 461 | CALL nc(NF90_PUT_ATT(fID,loID,'long_name','Longitude'),'x') |
|---|
| 462 | CALL nc(NF90_PUT_ATT(fID,loID,'units','degrees_east') ,'x') |
|---|
| 463 | |
|---|
| 464 | CALL nc(NF90_DEF_DIM(fID,'y',nlat_ou,yID)) |
|---|
| 465 | CALL nc(NF90_DEF_VAR(fID,'y',NF90_REAL,yID,laID) ,'y') |
|---|
| 466 | CALL nc(NF90_PUT_ATT(fID,laID,'long_name','Latitude') ,'y') |
|---|
| 467 | CALL nc(NF90_PUT_ATT(fID,laID,'units','degrees_north'),'y') |
|---|
| 468 | |
|---|
| 469 | CALL nc(NF90_DEF_VAR(fID,'mask',NF90_REAL,[xID,yID],mskID),'mask') |
|---|
| 470 | IF(fmsk=='noro') & |
|---|
| 471 | CALL nc(NF90_DEF_VAR(fID,'Zphi',NF90_REAL,[xID,yID],phiID),'Zphi') |
|---|
| 472 | CALL nc(NF90_DEF_VAR(fID,'Zmea',NF90_REAL,[xID,yID],meaID),'Zmea') |
|---|
| 473 | CALL nc(NF90_DEF_VAR(fID,'mu' ,NF90_REAL,[xID,yID], muID),'mu' ) |
|---|
| 474 | CALL nc(NF90_DEF_VAR(fID,'Zsig',NF90_REAL,[xID,yID],sigID),'Zsig') |
|---|
| 475 | CALL nc(NF90_DEF_VAR(fID,'Zgam',NF90_REAL,[xID,yID],gamID),'Zgam') |
|---|
| 476 | CALL nc(NF90_DEF_VAR(fID,'Zthe',NF90_REAL,[xID,yID],theID),'Zthe') |
|---|
| 477 | CALL nc(NF90_DEF_VAR(fID,'Zpic',NF90_REAL,[xID,yID],picID),'Zpic') |
|---|
| 478 | CALL nc(NF90_DEF_VAR(fID,'Zval',NF90_REAL,[xID,yID],valID),'Zval') |
|---|
| 479 | |
|---|
| 480 | CALL nc(NF90_PUT_ATT(fID,mskID,'long_name','Fractional land mask' ),'mask') |
|---|
| 481 | IF(fmsk=='noro') & |
|---|
| 482 | CALL nc(NF90_PUT_ATT(fID,phiID,'long_name','Geopotential' ),'Zphi') |
|---|
| 483 | CALL nc(NF90_PUT_ATT(fID,meaID,'long_name','Mean orography' ),'Zmea') |
|---|
| 484 | CALL nc(NF90_PUT_ATT(fID, muID,'long_name','Std deviation of sub-cell scales orography'),'mu' ) |
|---|
| 485 | CALL nc(NF90_PUT_ATT(fID,sigID,'long_name','Slope along principal axis' ),'Zsig') |
|---|
| 486 | CALL nc(NF90_PUT_ATT(fID,gamID,'long_name','Anisotropy (aspect ratio)' ),'Zgam') |
|---|
| 487 | CALL nc(NF90_PUT_ATT(fID,theID,'long_name','Orientation (principal axis)' ),'Zthe') |
|---|
| 488 | CALL nc(NF90_PUT_ATT(fID,picID,'long_name','Maximum height' ),'Zpic') |
|---|
| 489 | CALL nc(NF90_PUT_ATT(fID,valID,'long_name','Minimum height' ),'Zval') |
|---|
| 490 | |
|---|
| 491 | CALL nc(NF90_PUT_ATT(fID,mskID,'units','none' ),'mask') |
|---|
| 492 | IF(fmsk=='noro') & |
|---|
| 493 | CALL nc(NF90_PUT_ATT(fID,phiID,'units','m' ),'Zphi') |
|---|
| 494 | CALL nc(NF90_PUT_ATT(fID,meaID,'units','m' ),'Zmea') |
|---|
| 495 | CALL nc(NF90_PUT_ATT(fID, muID,'units','m' ),'mu' ) |
|---|
| 496 | CALL nc(NF90_PUT_ATT(fID,sigID,'units','m/m' ),'Zsig') |
|---|
| 497 | CALL nc(NF90_PUT_ATT(fID,gamID,'units','none' ),'Zgam') |
|---|
| 498 | CALL nc(NF90_PUT_ATT(fID,theID,'units','degrees'),'Zthe') |
|---|
| 499 | CALL nc(NF90_PUT_ATT(fID,picID,'units','m' ),'Zpic') |
|---|
| 500 | CALL nc(NF90_PUT_ATT(fID,valID,'units','m' ),'Zval') |
|---|
| 501 | |
|---|
| 502 | CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Conventions','COARDS/CF-1.0')) |
|---|
| 503 | CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Initial_Grid',TRIM(res_in))) |
|---|
| 504 | CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'history',TRIM(call_seq))) |
|---|
| 505 | CALL nc(NF90_ENDDEF(fID)) |
|---|
| 506 | |
|---|
| 507 | CALL nc(NF90_PUT_VAR(fID, loID,lon_ou),'x' ) |
|---|
| 508 | CALL nc(NF90_PUT_VAR(fID, laID,lat_ou),'y' ) |
|---|
| 509 | CALL nc(NF90_PUT_VAR(fID,mskID,msko),'mask') |
|---|
| 510 | IF(fmsk=='noro') & |
|---|
| 511 | CALL nc(NF90_PUT_VAR(fID,phiID,Zphi),'Zphi') |
|---|
| 512 | CALL nc(NF90_PUT_VAR(fID,meaID,h0 ),'Zmea') |
|---|
| 513 | CALL nc(NF90_PUT_VAR(fID, muID,mu ),'mu' ) |
|---|
| 514 | CALL nc(NF90_PUT_VAR(fID,sigID,Zsig),'Zsig') |
|---|
| 515 | CALL nc(NF90_PUT_VAR(fID,gamID,Zgam),'Zgam') |
|---|
| 516 | CALL nc(NF90_PUT_VAR(fID,theID,Zthe),'Zthe') |
|---|
| 517 | CALL nc(NF90_PUT_VAR(fID,picID, h0+2*mu ),'Zpic') |
|---|
| 518 | CALL nc(NF90_PUT_VAR(fID,valID,MAX(0.,h0-2*mu)),'Zval') |
|---|
| 519 | CALL nc(NF90_CLOSE(fID)) |
|---|
| 520 | WRITE(*,*)'Finished.' |
|---|
| 521 | |
|---|
| 522 | CONTAINS |
|---|
| 523 | |
|---|
| 524 | |
|---|
| 525 | !------------------------------------------------------------------------------- |
|---|
| 526 | ! |
|---|
| 527 | SUBROUTINE diffusion_filter(a,b) |
|---|
| 528 | ! |
|---|
| 529 | !------------------------------------------------------------------------------- |
|---|
| 530 | ! Purpose: ECMWF bi-laplacian diffusion filter to limit aliasing. |
|---|
| 531 | ! Note that one third of the frequencies is skipped to allow some accuracy for |
|---|
| 532 | ! derived quantities computation. |
|---|
| 533 | !------------------------------------------------------------------------------- |
|---|
| 534 | IMPLICIT NONE |
|---|
| 535 | !------------------------------------------------------------------------------- |
|---|
| 536 | ! Arguments: |
|---|
| 537 | REAL, INTENT(INOUT) :: a(:,:), b(:,:) |
|---|
| 538 | !------------------------------------------------------------------------------- |
|---|
| 539 | ! Local variables: |
|---|
| 540 | INTEGER :: m |
|---|
| 541 | REAL, ALLOCATABLE :: fn(:) |
|---|
| 542 | !------------------------------------------------------------------------------- |
|---|
| 543 | ALLOCATE(fn(ndab)); fn=[(REAL(m*(m-1)),m=1,ndab)] |
|---|
| 544 | a(1,:)=a(1,:)/(1+4.*fn(:)/fn(ndab))**2 |
|---|
| 545 | b(1,:)=b(1,:)/(1+4.*fn(:)/fn(ndab))**2 |
|---|
| 546 | DO m=1,ndab |
|---|
| 547 | IF(m<=2*ndab/3) THEN |
|---|
| 548 | a(m,:)=a(m,:)/(1+5.*fn(:)/fn(ndab))**2 |
|---|
| 549 | b(m,:)=b(m,:)/(1+5.*fn(:)/fn(ndab))**2 |
|---|
| 550 | ELSE |
|---|
| 551 | a(m,:)=0. |
|---|
| 552 | b(m,:)=0. |
|---|
| 553 | END IF |
|---|
| 554 | END DO |
|---|
| 555 | |
|---|
| 556 | END SUBROUTINE diffusion_filter |
|---|
| 557 | ! |
|---|
| 558 | !------------------------------------------------------------------------------- |
|---|
| 559 | |
|---|
| 560 | |
|---|
| 561 | !------------------------------------------------------------------------------- |
|---|
| 562 | ! |
|---|
| 563 | SUBROUTINE lowpass_filter(a,b,d,dd) |
|---|
| 564 | ! |
|---|
| 565 | !------------------------------------------------------------------------------- |
|---|
| 566 | ! Purpose: ECMWF lowpass filter. |
|---|
| 567 | !------------------------------------------------------------------------------- |
|---|
| 568 | IMPLICIT NONE |
|---|
| 569 | !------------------------------------------------------------------------------- |
|---|
| 570 | ! Arguments: |
|---|
| 571 | REAL, INTENT(INOUT) :: a(:,:), b(:,:) |
|---|
| 572 | REAL, INTENT(IN) :: d, dd |
|---|
| 573 | !------------------------------------------------------------------------------- |
|---|
| 574 | ! Local variables: |
|---|
| 575 | INTEGER :: m, na |
|---|
| 576 | REAL, ALLOCATABLE :: h(:) |
|---|
| 577 | !------------------------------------------------------------------------------- |
|---|
| 578 | na=SIZE(a,1) |
|---|
| 579 | IF(ALLOCATED(h)) THEN; IF(SIZE(h)/=na) DEALLOCATE(h); END IF |
|---|
| 580 | IF(.NOT.ALLOCATED(h)) ALLOCATE(h(na)) |
|---|
| 581 | CALL lowpass(d,dd,h) |
|---|
| 582 | DO m=1,SIZE(h); a(m,:)=a(m,:)*h; b(m,:)=b(m,:)*h; END DO |
|---|
| 583 | |
|---|
| 584 | END SUBROUTINE lowpass_filter |
|---|
| 585 | ! |
|---|
| 586 | !------------------------------------------------------------------------------- |
|---|
| 587 | |
|---|
| 588 | |
|---|
| 589 | !------------------------------------------------------------------------------- |
|---|
| 590 | ! |
|---|
| 591 | SUBROUTINE bandpass_filter(a,b,d1,dd1,d2,dd2) |
|---|
| 592 | ! |
|---|
| 593 | !------------------------------------------------------------------------------- |
|---|
| 594 | ! Purpose: ECMWF bandpass filter. Calls lowpass filter with 2 different scales. |
|---|
| 595 | !------------------------------------------------------------------------------- |
|---|
| 596 | IMPLICIT NONE |
|---|
| 597 | !------------------------------------------------------------------------------- |
|---|
| 598 | ! Arguments: |
|---|
| 599 | REAL, INTENT(INOUT) :: a(:,:), b(:,:) |
|---|
| 600 | REAL, INTENT(IN) :: d1, dd1, d2, dd2 |
|---|
| 601 | !------------------------------------------------------------------------------- |
|---|
| 602 | ! Local variables: |
|---|
| 603 | INTEGER :: m, na |
|---|
| 604 | REAL, ALLOCATABLE :: h(:), hh(:) |
|---|
| 605 | !------------------------------------------------------------------------------- |
|---|
| 606 | na=SIZE(a,1) |
|---|
| 607 | IF(ALLOCATED(hh)) THEN; IF(SIZE(hh)/=na) DEALLOCATE(hh); END IF |
|---|
| 608 | IF(ALLOCATED(h )) THEN; IF(SIZE(h )/=na) DEALLOCATE(h ); END IF |
|---|
| 609 | IF(.NOT.ALLOCATED(hh)) ALLOCATE(hh(na)) |
|---|
| 610 | IF(.NOT.ALLOCATED(h )) ALLOCATE(h (na)) |
|---|
| 611 | CALL lowpass(d1,dd1,hh) |
|---|
| 612 | CALL lowpass(d2,dd2,h ); hh=hh-h |
|---|
| 613 | DO m=1,SIZE(h); a(m,:)=a(m,:)*hh; b(m,:)=b(m,:)*hh; END DO |
|---|
| 614 | |
|---|
| 615 | END SUBROUTINE bandpass_filter |
|---|
| 616 | ! |
|---|
| 617 | !------------------------------------------------------------------------------- |
|---|
| 618 | |
|---|
| 619 | |
|---|
| 620 | !------------------------------------------------------------------------------- |
|---|
| 621 | ! |
|---|
| 622 | SUBROUTINE lowpass(d,dd,h) |
|---|
| 623 | ! |
|---|
| 624 | !------------------------------------------------------------------------------- |
|---|
| 625 | IMPLICIT NONE |
|---|
| 626 | !------------------------------------------------------------------------------- |
|---|
| 627 | ! Arguments: |
|---|
| 628 | REAL, INTENT(IN) :: d, dd ! d is DELTA/REarth ; dd is delta/REarth |
|---|
| 629 | REAL, INTENT(INOUT) :: h(:) |
|---|
| 630 | !------------------------------------------------------------------------------- |
|---|
| 631 | INTEGER :: m, n |
|---|
| 632 | REAL :: hd |
|---|
| 633 | !------------------------------------------------------------------------------- |
|---|
| 634 | n=SIZE(h); hd=0.5*d |
|---|
| 635 | h(:)=[((SIN(m*(hd-dd))/m+(COS(hp+m*hd)*SIN(hp+m*dd))/(hp/dd+m) & |
|---|
| 636 | + SIN(m*(hd+dd))/m+(COS(hp-m*hd)*SIN(hp-m*dd))/(hp/dd-m))**2/d/d,m=1,n)] |
|---|
| 637 | |
|---|
| 638 | END SUBROUTINE lowpass |
|---|
| 639 | ! |
|---|
| 640 | !------------------------------------------------------------------------------- |
|---|
| 641 | |
|---|
| 642 | |
|---|
| 643 | !------------------------------------------------------------------------------- |
|---|
| 644 | ! |
|---|
| 645 | SUBROUTINE spatial_filter_1(u,D,dd) |
|---|
| 646 | ! |
|---|
| 647 | !------------------------------------------------------------------------------- |
|---|
| 648 | ! Purpose: Spatial filter ; using spectral lowpass filter. |
|---|
| 649 | !------------------------------------------------------------------------------- |
|---|
| 650 | IMPLICIT NONE |
|---|
| 651 | !------------------------------------------------------------------------------- |
|---|
| 652 | ! Arguments: |
|---|
| 653 | REAL, INTENT(INOUT) :: u(:,:) |
|---|
| 654 | REAL, INTENT(IN) :: D, dd |
|---|
| 655 | !------------------------------------------------------------------------------- |
|---|
| 656 | ! Local variables: |
|---|
| 657 | REAL, ALLOCATABLE :: uu(:,:,:) |
|---|
| 658 | !------------------------------------------------------------------------------- |
|---|
| 659 | ALLOCATE(uu(SIZE(u,1),SIZE(u,2),1)); uu(:,:,1)=u |
|---|
| 660 | CALL spatial_filter_m(uu,D,dd) |
|---|
| 661 | u=uu(:,:,1) |
|---|
| 662 | DEALLOCATE(uu) |
|---|
| 663 | |
|---|
| 664 | END SUBROUTINE spatial_filter_1 |
|---|
| 665 | ! |
|---|
| 666 | !------------------------------------------------------------------------------- |
|---|
| 667 | |
|---|
| 668 | |
|---|
| 669 | !------------------------------------------------------------------------------- |
|---|
| 670 | ! |
|---|
| 671 | SUBROUTINE spatial_filter_m(u,D,dd) |
|---|
| 672 | ! |
|---|
| 673 | !------------------------------------------------------------------------------- |
|---|
| 674 | ! Purpose: Spatial filter ; using spectral lowpass filter. |
|---|
| 675 | !------------------------------------------------------------------------------- |
|---|
| 676 | IMPLICIT NONE |
|---|
| 677 | !------------------------------------------------------------------------------- |
|---|
| 678 | ! Arguments: |
|---|
| 679 | REAL, INTENT(INOUT) :: u(:,:,:) |
|---|
| 680 | REAL, INTENT(IN) :: D, dd |
|---|
| 681 | !------------------------------------------------------------------------------- |
|---|
| 682 | ! Local variables: |
|---|
| 683 | INTEGER :: ku |
|---|
| 684 | REAL, ALLOCATABLE :: aa(:,:,:),bb(:,:,:) |
|---|
| 685 | !------------------------------------------------------------------------------- |
|---|
| 686 | ku=SIZE(u,3) !--- NUMBER OF FIELDS TO SMOOTH |
|---|
| 687 | ALLOCATE(aa(mdab,ndab,ku),bb(mdab,ndab,ku)) |
|---|
| 688 | |
|---|
| 689 | !--- PROJECTION OF THE SUB-CELL SCALES FIELDS ON THE SPHERICAL HARMONICS BASE. |
|---|
| 690 | IF(lfast) THEN |
|---|
| 691 | CALL shaec(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wsha,lsh,work,lwork,ierr) |
|---|
| 692 | CALL error('shaec',ierr) |
|---|
| 693 | ELSE |
|---|
| 694 | CALL shaes(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wsha,lsh,work,lwork,ierr) |
|---|
| 695 | CALL error('shaes',ierr) |
|---|
| 696 | END IF |
|---|
| 697 | |
|---|
| 698 | !--- LOW PASS FILTER AT THE OUTPUT GRID RESOLUTION |
|---|
| 699 | DO n=1,ku; CALL lowpass_filter(aa(:,:,n),bb(:,:,n),D,dd); END DO |
|---|
| 700 | |
|---|
| 701 | !--- BACKWARD TRANSFORM FOR MEAN SUB-CELL SCALES PARAMETERS |
|---|
| 702 | IF(lfast) THEN |
|---|
| 703 | CALL shsec(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wshs,lsh,work,lwork,ierr) |
|---|
| 704 | CALL error('shsec',ierr) |
|---|
| 705 | ELSE |
|---|
| 706 | CALL shses(nlat_in,nlon_in,0,ku,u,idh,jdh,aa,bb,mdab,ndab,wshs,lsh,work,lwork,ierr) |
|---|
| 707 | CALL error('shses',ierr) |
|---|
| 708 | END IF |
|---|
| 709 | DEALLOCATE(aa,bb) |
|---|
| 710 | |
|---|
| 711 | END SUBROUTINE spatial_filter_m |
|---|
| 712 | ! |
|---|
| 713 | !------------------------------------------------------------------------------- |
|---|
| 714 | |
|---|
| 715 | |
|---|
| 716 | !------------------------------------------------------------------------------- |
|---|
| 717 | ! |
|---|
| 718 | SUBROUTINE grid_noro0(xin,yin,zin,xou,you,zphi,mask) |
|---|
| 719 | ! |
|---|
| 720 | !------------------------------------------------------------------------------- |
|---|
| 721 | ! Purpose: Sub-cell scales orographic parameters coomputation. Angles in radians |
|---|
| 722 | ! Notes: Input coordinates contain a longitude periodic redundant point. |
|---|
| 723 | !------------------------------------------------------------------------------- |
|---|
| 724 | IMPLICIT NONE |
|---|
| 725 | !------------------------------------------------------------------------------- |
|---|
| 726 | ! Arguments: |
|---|
| 727 | REAL, INTENT(IN) :: xin(:), yin(:) !--- INPUT COORDINATES (nxi) (nyi) |
|---|
| 728 | REAL, INTENT(IN) :: zin(:,:) !--- INPUT FIELD (nxi, nyi) |
|---|
| 729 | REAL, INTENT(IN) :: xou(:), you(:) !--- OUTPUT COORDINATES (nxo+1) (nyo) |
|---|
| 730 | REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (nxo+1, nyo) |
|---|
| 731 | REAL, INTENT(OUT) :: mask(:,:) !--- FRACTIONAL MASK (nxo+1, nyo) |
|---|
| 732 | !------------------------------------------------------------------------------- |
|---|
| 733 | ! Local variables: |
|---|
| 734 | REAL, ALLOCATABLE :: xi(:), xo_w(:), yo_s(:), zi(:,:), num_tot(:,:) |
|---|
| 735 | REAL, ALLOCATABLE :: yi(:), xo_e(:), yo_n(:), w(:,:), num_lan(:,:) |
|---|
| 736 | REAL :: dx, dy, wx, wy, zx, zy, ex, ey |
|---|
| 737 | INTEGER :: nxi, nxo, i, ii, nn |
|---|
| 738 | INTEGER :: nyi, nyo, j, jj, di |
|---|
| 739 | !------------------------------------------------------------------------------- |
|---|
| 740 | nxi=SIZE(xin); nxo=SIZE(xou)-1; ex=SIGN(1.,xin(nxi)-xin(1)) |
|---|
| 741 | nyi=SIZE(yin); nyo=SIZE(you); ey=SIGN(1.,yin(nyi)-yin(1)) |
|---|
| 742 | |
|---|
| 743 | !--- MARGIN TO ENSURE OUTPUT CELLS HAVE MATCHING INPUT CELLS |
|---|
| 744 | di=CEILING(1.+MAXVAL([xou(2:nxo+1)-xou(1:nxo)])/MINVAL([xin(2:nxi)-xin(1:nxi-1)])) |
|---|
| 745 | |
|---|
| 746 | !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES |
|---|
| 747 | ALLOCATE(xi(nxi+2*di)); xi(:)=[xin(nxi-di:nxi-1)-2.*xpi,xin,xin(2:di+1)+2.*xpi] |
|---|
| 748 | ALLOCATE(yi(nyi+2)); yi(:)=[2.*yin(1)-yin(2),yin,2.*yin(nyi)-yin(nyi-1)] |
|---|
| 749 | |
|---|
| 750 | ALLOCATE(zi(nxi+2*di,nyi+2)) |
|---|
| 751 | zi(1 +di:nxi+ di,2:nyi+1)=zin(:,:) |
|---|
| 752 | zi(1 : di,2:nyi+1)=zin(nxi-di:nxi-1,:) |
|---|
| 753 | zi(1+nxi +di:nxi+2*di,2:nyi+1)=zin(2 :di+1 ,:) |
|---|
| 754 | zi(1 :nxi/2+di, 1)=zi (1+nxi/2:nxi +di, 2) |
|---|
| 755 | zi(1+nxi/2+di:nxi+2*di, 1)=zi (1 :nxi/2+di, 2) |
|---|
| 756 | zi(1 :nxi/2+di, nyi+2)=zi (1+nxi/2:nxi +di,nyi+1) |
|---|
| 757 | zi(1+nxi/2+di:nxi+2*di, nyi+2)=zi (1 :nxi/2+di,nyi+1) |
|---|
| 758 | |
|---|
| 759 | ! write(*,'(a,2f10.5,a,2f10.5)')'xin=',xin(1:2),'...',xin(nxi-1:nxi) |
|---|
| 760 | ! write(*,'(a,2f10.5,a,2f10.5)')'yin=',yin(1:2),'...',yin(nyi-1:nyi) |
|---|
| 761 | ! write(*,'(a,2f10.5,a,2f10.5)')'xou=',xou(1:2),'...',xou(nxo:nxo+1) |
|---|
| 762 | ! write(*,'(a,2f10.5,a,2f10.5)')'you=',you(1:2),'...',you(nyo-1:nyo) |
|---|
| 763 | ! write(*,*)'xi=',xi(1:di+1),'...',xi(nxi+di:nxi+2*di) |
|---|
| 764 | ! write(*,'(a,3f10.5,a,3f10.5)')'yi=',yi(1:3),'...',yi(nyi:nyi+2) |
|---|
| 765 | |
|---|
| 766 | !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) bdac |
|---|
| 767 | ALLOCATE(xo_w(nxo+1)); xo_w=(xou+[2.*xou( 1)-xou( 2),xou(1:nxo )])/2.0 |
|---|
| 768 | ALLOCATE(yo_s(nyo )); yo_s=(you+[2.*you( 1)-you( 2),you(1:nyo-1)])/2.0 |
|---|
| 769 | ALLOCATE(xo_e(nxo+1)); xo_e=(xou+[xou(2:nxo+1),2.*xou(nxo+1)-xou(nxo)])/2.0 |
|---|
| 770 | ALLOCATE(yo_n(nyo )); yo_n=(you+[you(2:nyo ),2.*you(nyo)-you(nyo-1)])/2.0 |
|---|
| 771 | |
|---|
| 772 | !--- INITIALIZATIONS: |
|---|
| 773 | ALLOCATE(w(nxo+1,nyo)); w(:,:)=0.0; zphi(:,:)=0.0 |
|---|
| 774 | |
|---|
| 775 | !--- SUMMATION OVER GRIDPOINT AREA |
|---|
| 776 | zx=2.*xpi/REAL(nxi) !--- LON CELL LENGTH (1-SPHERE) |
|---|
| 777 | dx=zx/2. !--- HALF " |
|---|
| 778 | ALLOCATE(num_tot(nxo+1,nyo)); num_tot(:,:)=0. |
|---|
| 779 | ALLOCATE(num_lan(nxo+1,nyo)); num_lan(:,:)=0. |
|---|
| 780 | DO ii=1,nxo+1 |
|---|
| 781 | DO jj=1,nyo |
|---|
| 782 | DO j=2,nyi+1 |
|---|
| 783 | zy=xpi/REAL(nyi)*COS(yi(j)) !--- LAT CELL LENGTH (1-SPHERE) |
|---|
| 784 | dy=zy/2. |
|---|
| 785 | wy=MAX(0.,MIN(zy,MINVAL(dy+ey*[yo_n(jj)-yi(j),yi(j)-yo_s(jj)]))) |
|---|
| 786 | IF(wy==0.) CYCLE |
|---|
| 787 | DO i=2,nxi+2*di-1 |
|---|
| 788 | wx=MAX(0.,MIN(zx,MINVAL(dx+ex*[xo_e(ii)-xi(i),xi(i)-xo_w(ii)])*COS(yi(j)))) |
|---|
| 789 | IF(wx==0.) CYCLE |
|---|
| 790 | num_tot(ii,jj)=num_tot(ii,jj)+1.0 |
|---|
| 791 | IF(zi(i,j)>=1.) num_lan(ii,jj)=num_lan(ii,jj)+1.0 |
|---|
| 792 | w (ii,jj)= w(ii,jj)+wx*wy |
|---|
| 793 | zphi(ii,jj)=zphi(ii,jj)+wx*wy*zi(i,j) |
|---|
| 794 | END DO |
|---|
| 795 | END DO |
|---|
| 796 | END DO |
|---|
| 797 | END DO |
|---|
| 798 | |
|---|
| 799 | !--- COMPUTE FRACTIONAL MASK AND GEOPOTENTIAL |
|---|
| 800 | WHERE(w(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:) |
|---|
| 801 | nn=COUNT(w(:,:)==0.0) |
|---|
| 802 | IF(nn/=0) WRITE(*,*)'Problem with weight ; vanishing occurrences: ',nn |
|---|
| 803 | WHERE(w/=0.0) zphi=zphi/w |
|---|
| 804 | WHERE(w==0.0.OR.mask<0.1) zphi=0.0 |
|---|
| 805 | zphi(nxo+1,:)=zphi(1,:) |
|---|
| 806 | zphi(:, 1)=SUM(w(1:nxo, 1)*zphi(1:nxo, 1))/SUM(w(1:nxo, 1)) |
|---|
| 807 | zphi(:,nyo)=SUM(w(1:nxo,nyo)*zphi(1:nxo,nyo))/SUM(w(1:nxo,nyo)) |
|---|
| 808 | |
|---|
| 809 | END SUBROUTINE grid_noro0 |
|---|
| 810 | ! |
|---|
| 811 | !------------------------------------------------------------------------------- |
|---|
| 812 | |
|---|
| 813 | |
|---|
| 814 | !------------------------------------------------------------------------------- |
|---|
| 815 | ! |
|---|
| 816 | SUBROUTINE nc(ncres,var) |
|---|
| 817 | ! |
|---|
| 818 | !------------------------------------------------------------------------------- |
|---|
| 819 | ! Purpose: NetCDF errors handling. |
|---|
| 820 | !------------------------------------------------------------------------------- |
|---|
| 821 | IMPLICIT NONE |
|---|
| 822 | !------------------------------------------------------------------------------- |
|---|
| 823 | ! Arguments: |
|---|
| 824 | INTEGER, INTENT(IN) :: ncres |
|---|
| 825 | CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: var |
|---|
| 826 | !------------------------------------------------------------------------------- |
|---|
| 827 | ! Local variables: |
|---|
| 828 | CHARACTER(LEN=256) :: msg |
|---|
| 829 | !------------------------------------------------------------------------------- |
|---|
| 830 | IF(ncres/=NF90_NoErr) THEN |
|---|
| 831 | msg='Error in routine '//TRIM(sub) |
|---|
| 832 | IF(fnam/='') msg=TRIM(msg)//' for file "'//TRIM(fnam)//'"' |
|---|
| 833 | IF(PRESENT(var)) msg=TRIM(msg)//'" and variable "'//TRIM(var)//'"' |
|---|
| 834 | WRITE(*,*)TRIM(msg)//': '//NF90_STRERROR(ncres); STOP |
|---|
| 835 | END IF |
|---|
| 836 | |
|---|
| 837 | END SUBROUTINE nc |
|---|
| 838 | ! |
|---|
| 839 | !------------------------------------------------------------------------------- |
|---|
| 840 | |
|---|
| 841 | |
|---|
| 842 | !------------------------------------------------------------------------------- |
|---|
| 843 | ! |
|---|
| 844 | ELEMENTAL FUNCTION int2str(val) |
|---|
| 845 | ! |
|---|
| 846 | !------------------------------------------------------------------------------- |
|---|
| 847 | IMPLICIT NONE |
|---|
| 848 | !------------------------------------------------------------------------------- |
|---|
| 849 | ! Purpose: convert an integer into a string. |
|---|
| 850 | !------------------------------------------------------------------------------- |
|---|
| 851 | ! Arguments: |
|---|
| 852 | CHARACTER(LEN=128) :: int2str |
|---|
| 853 | INTEGER, INTENT(IN) :: val |
|---|
| 854 | !------------------------------------------------------------------------------- |
|---|
| 855 | ! Local variable: |
|---|
| 856 | CHARACTER(LEN=128) :: s1 |
|---|
| 857 | !------------------------------------------------------------------------------- |
|---|
| 858 | WRITE(s1,*)val; int2str=TRIM(ADJUSTL(s1)) |
|---|
| 859 | |
|---|
| 860 | END FUNCTION int2str |
|---|
| 861 | !------------------------------------------------------------------------------- |
|---|
| 862 | |
|---|
| 863 | |
|---|
| 864 | !------------------------------------------------------------------------------- |
|---|
| 865 | ! |
|---|
| 866 | SUBROUTINE err(ierr,str) |
|---|
| 867 | ! |
|---|
| 868 | !------------------------------------------------------------------------------- |
|---|
| 869 | IMPLICIT NONE |
|---|
| 870 | !------------------------------------------------------------------------------- |
|---|
| 871 | ! Argument: |
|---|
| 872 | LOGICAL, INTENT(IN) :: ierr |
|---|
| 873 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 874 | !------------------------------------------------------------------------------- |
|---|
| 875 | IF(.NOT.ierr) RETURN |
|---|
| 876 | WRITE(*,*)TRIM(str) |
|---|
| 877 | STOP |
|---|
| 878 | |
|---|
| 879 | END SUBROUTINE err |
|---|
| 880 | ! |
|---|
| 881 | !------------------------------------------------------------------------------- |
|---|
| 882 | |
|---|
| 883 | |
|---|
| 884 | !------------------------------------------------------------------------------- |
|---|
| 885 | ! |
|---|
| 886 | SUBROUTINE error(sub,ierr) |
|---|
| 887 | ! |
|---|
| 888 | !------------------------------------------------------------------------------- |
|---|
| 889 | IMPLICIT NONE |
|---|
| 890 | !------------------------------------------------------------------------------- |
|---|
| 891 | ! Arguments: |
|---|
| 892 | CHARACTER(LEN=*), INTENT(IN) :: sub |
|---|
| 893 | INTEGER, INTENT(IN) :: ierr |
|---|
| 894 | !------------------------------------------------------------------------------- |
|---|
| 895 | IF(ierr==0) RETURN |
|---|
| 896 | WRITE(*,*)'In routine "'//TRIM(sub)//'": error in the specification of argume& |
|---|
| 897 | &nt nr. '//TRIM(int2str(ierr)) |
|---|
| 898 | STOP |
|---|
| 899 | END SUBROUTINE error |
|---|
| 900 | ! |
|---|
| 901 | !------------------------------------------------------------------------------- |
|---|
| 902 | |
|---|
| 903 | |
|---|
| 904 | !------------------------------------------------------------------------------- |
|---|
| 905 | ! |
|---|
| 906 | ELEMENTAL FUNCTION real2str(val) |
|---|
| 907 | ! |
|---|
| 908 | !------------------------------------------------------------------------------- |
|---|
| 909 | IMPLICIT NONE |
|---|
| 910 | !------------------------------------------------------------------------------- |
|---|
| 911 | ! Purpose: convert an integer into a string. |
|---|
| 912 | !------------------------------------------------------------------------------- |
|---|
| 913 | ! Arguments: |
|---|
| 914 | CHARACTER(LEN=128) :: real2str |
|---|
| 915 | REAL, INTENT(IN) :: val |
|---|
| 916 | INTEGER :: r, s |
|---|
| 917 | !------------------------------------------------------------------------------- |
|---|
| 918 | s=NINT(val); r=NINT(1000.*(val-REAL(s))) |
|---|
| 919 | IF(r<0) THEN; s= INT(val); r= INT(1000.*(val-REAL(s))); END IF |
|---|
| 920 | real2str=TRIM(ADJUSTL(int2str(s)))//'.'//TRIM(ADJUSTL(int2str(r))) |
|---|
| 921 | |
|---|
| 922 | END FUNCTION real2str |
|---|
| 923 | !------------------------------------------------------------------------------- |
|---|
| 924 | |
|---|
| 925 | |
|---|
| 926 | !------------------------------------------------------------------------------- |
|---|
| 927 | ! |
|---|
| 928 | SUBROUTINE flip(a) |
|---|
| 929 | ! |
|---|
| 930 | !------------------------------------------------------------------------------- |
|---|
| 931 | IMPLICIT NONE |
|---|
| 932 | !------------------------------------------------------------------------------- |
|---|
| 933 | ! Arguments: |
|---|
| 934 | REAL, ALLOCATABLE, INTENT(INOUT) :: a(:,:) |
|---|
| 935 | !------------------------------------------------------------------------------- |
|---|
| 936 | ! Local variables: |
|---|
| 937 | INTEGER :: k, n |
|---|
| 938 | REAL, ALLOCATABLE :: tmp(:,:) |
|---|
| 939 | !------------------------------------------------------------------------------- |
|---|
| 940 | n=SIZE(a,2) |
|---|
| 941 | ALLOCATE(tmp(n,SIZE(a,1))); DO k=1,n; tmp(k,:)=a(:,k); END DO |
|---|
| 942 | CALL MOVE_ALLOC(FROM=tmp,TO=a) |
|---|
| 943 | |
|---|
| 944 | END SUBROUTINE flip |
|---|
| 945 | ! |
|---|
| 946 | !------------------------------------------------------------------------------- |
|---|
| 947 | |
|---|
| 948 | |
|---|
| 949 | !------------------------------------------------------------------------------- |
|---|
| 950 | ELEMENTAL LOGICAL FUNCTION is_int(str) RESULT(out) |
|---|
| 951 | !------------------------------------------------------------------------------- |
|---|
| 952 | IMPLICIT NONE |
|---|
| 953 | !------------------------------------------------------------------------------- |
|---|
| 954 | ! Purpose: tests if str is a valid integer number |
|---|
| 955 | !------------------------------------------------------------------------------- |
|---|
| 956 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 957 | INTEGER :: k, k0 |
|---|
| 958 | !------------------------------------------------------------------------------- |
|---|
| 959 | out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN |
|---|
| 960 | k0=1; IF(INDEX('+-',str(1:1))/=0) k0=2 |
|---|
| 961 | DO k=k0,LEN_TRIM(str); IF(INDEX('0123456789',str(k:k))==0) out=.FALSE.; END DO |
|---|
| 962 | END FUNCTION is_int |
|---|
| 963 | !------------------------------------------------------------------------------- |
|---|
| 964 | |
|---|
| 965 | !------------------------------------------------------------------------------- |
|---|
| 966 | ELEMENTAL LOGICAL FUNCTION is_dec(str) RESULT(out) |
|---|
| 967 | !------------------------------------------------------------------------------- |
|---|
| 968 | IMPLICIT NONE |
|---|
| 969 | !------------------------------------------------------------------------------- |
|---|
| 970 | ! Purpose: tests if str is a valid decimal number or an integer |
|---|
| 971 | !------------------------------------------------------------------------------- |
|---|
| 972 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 973 | INTEGER :: ix |
|---|
| 974 | !------------------------------------------------------------------------------- |
|---|
| 975 | out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN |
|---|
| 976 | |
|---|
| 977 | !--- TEST IF THE NUMBER IS A PURE INTEGER |
|---|
| 978 | IF(is_int(str)) RETURN |
|---|
| 979 | |
|---|
| 980 | !--- CHECK FOR POINT SEPARATOR (SHOULD HAVE A STRING INT BEFORE AND/OR AFTER |
|---|
| 981 | ix=INDEX(str,'.'); out=ix>0.AND.LEN_TRIM(str)/=1; IF(.NOT.out) RETURN |
|---|
| 982 | |
|---|
| 983 | !--- CHECK SURROUNDING STRINGS ARE INTEGERS |
|---|
| 984 | out=is_int(str(1:ix-1)); IF(ix==LEN_TRIM(str)) RETURN |
|---|
| 985 | out=out.AND.is_int(str(ix+1:LEN_TRIM(str))) |
|---|
| 986 | END FUNCTION is_dec |
|---|
| 987 | !------------------------------------------------------------------------------- |
|---|
| 988 | |
|---|
| 989 | !------------------------------------------------------------------------------- |
|---|
| 990 | ELEMENTAL LOGICAL FUNCTION is_num(str) RESULT(out) |
|---|
| 991 | !------------------------------------------------------------------------------- |
|---|
| 992 | IMPLICIT NONE |
|---|
| 993 | !------------------------------------------------------------------------------- |
|---|
| 994 | ! Purpose: test if str is a valid integer/decimal/scientific notation number |
|---|
| 995 | !------------------------------------------------------------------------------- |
|---|
| 996 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 997 | INTEGER :: ix, ls |
|---|
| 998 | !------------------------------------------------------------------------------- |
|---|
| 999 | out=LEN_TRIM(str)/=0; IF(.NOT.out) RETURN; ls=LEN_TRIM(str) |
|---|
| 1000 | |
|---|
| 1001 | !--- TEST IF THE NUMBER IS A PURE INTEGER OR A PURE DECIMAL NUMBER |
|---|
| 1002 | IF(is_int(str)) RETURN |
|---|
| 1003 | IF(is_dec(str)) RETURN |
|---|
| 1004 | |
|---|
| 1005 | !--- CHECK FOR "D" OR "E" SEPARATOR |
|---|
| 1006 | ix=MAX(INDEX(strLow(str),'d'),INDEX(strLow(str),'e')); out=ix>0 |
|---|
| 1007 | IF(.NOT.out) RETURN |
|---|
| 1008 | |
|---|
| 1009 | !--- CHECK SURROUNDING STRINGS ARE NUMBERS |
|---|
| 1010 | out=is_dec(str(1:ix-1)).AND.is_int(str(ix+1:ls)) |
|---|
| 1011 | END FUNCTION is_num |
|---|
| 1012 | !------------------------------------------------------------------------------- |
|---|
| 1013 | |
|---|
| 1014 | !------------------------------------------------------------------------------- |
|---|
| 1015 | INTEGER FUNCTION str2int(str) RESULT(out) |
|---|
| 1016 | !------------------------------------------------------------------------------- |
|---|
| 1017 | IMPLICIT NONE |
|---|
| 1018 | !------------------------------------------------------------------------------- |
|---|
| 1019 | ! Purpose: convert a string into an integer |
|---|
| 1020 | !------------------------------------------------------------------------------- |
|---|
| 1021 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 1022 | INTEGER :: k |
|---|
| 1023 | !------------------------------------------------------------------------------- |
|---|
| 1024 | CALL err(.NOT.is_int(str),'Invalid integer "'//TRIM(str)//'"') |
|---|
| 1025 | out=0 |
|---|
| 1026 | DO k=1,LEN_TRIM(str); out=10*out+INT(INDEX('0123456789',str(k:k))-1); END DO |
|---|
| 1027 | END FUNCTION str2int |
|---|
| 1028 | !------------------------------------------------------------------------------- |
|---|
| 1029 | |
|---|
| 1030 | !------------------------------------------------------------------------------- |
|---|
| 1031 | REAL FUNCTION str2real(str) RESULT(out) |
|---|
| 1032 | !------------------------------------------------------------------------------- |
|---|
| 1033 | IMPLICIT NONE |
|---|
| 1034 | !------------------------------------------------------------------------------- |
|---|
| 1035 | ! Purpose: Convert a string into an real |
|---|
| 1036 | !------------------------------------------------------------------------------- |
|---|
| 1037 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 1038 | INTEGER :: ix, l |
|---|
| 1039 | !------------------------------------------------------------------------------- |
|---|
| 1040 | CALL err(.NOT.is_num(str),'Invalid number "'//TRIM(str)//'"') |
|---|
| 1041 | ix=MAX(INDEX(strLow(str),'d'),INDEX(strLow(str),'e')); l=LEN_TRIM(str) |
|---|
| 1042 | IF(ix==0) out=dec2real(str) |
|---|
| 1043 | IF(ix/=0) out=dec2real(str(1:ix-1))*10.**str2int(str(ix+1:l)) |
|---|
| 1044 | END FUNCTION str2real |
|---|
| 1045 | !------------------------------------------------------------------------------- |
|---|
| 1046 | |
|---|
| 1047 | !------------------------------------------------------------------------------- |
|---|
| 1048 | ELEMENTAL REAL FUNCTION dec2real(str) RESULT(out) |
|---|
| 1049 | !------------------------------------------------------------------------------- |
|---|
| 1050 | IMPLICIT NONE |
|---|
| 1051 | !------------------------------------------------------------------------------- |
|---|
| 1052 | ! Purpose: convert a valid (already checked) string into an real |
|---|
| 1053 | !------------------------------------------------------------------------------- |
|---|
| 1054 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 1055 | INTEGER :: k, k0, ix |
|---|
| 1056 | !------------------------------------------------------------------------------- |
|---|
| 1057 | out=0. |
|---|
| 1058 | k0=1; IF(INDEX('+-',str(1:1))/=0) k0=2; ix=INDEX(str,'.') !--- CHECK SIGN |
|---|
| 1059 | DO k=k0,LEN_TRIM(str); IF(k==ix) CYCLE !--- SKIP .+- |
|---|
| 1060 | out=out*10.+REAL(INDEX('0123456789',str(k:k))-1) |
|---|
| 1061 | END DO |
|---|
| 1062 | IF(ix/=0) out=out*10.**(ix-LEN_TRIM(str)) !--- COMA |
|---|
| 1063 | IF(k0==2) out=-out !--- SIGN |
|---|
| 1064 | END FUNCTION dec2real |
|---|
| 1065 | !------------------------------------------------------------------------------- |
|---|
| 1066 | |
|---|
| 1067 | !------------------------------------------------------------------------------- |
|---|
| 1068 | ! |
|---|
| 1069 | ELEMENTAL CHARACTER(LEN=256) FUNCTION strLow(str) RESULT(out) |
|---|
| 1070 | ! |
|---|
| 1071 | !------------------------------------------------------------------------------- |
|---|
| 1072 | IMPLICIT NONE |
|---|
| 1073 | !------------------------------------------------------------------------------- |
|---|
| 1074 | ! Purpose: Lower case conversion. |
|---|
| 1075 | !------------------------------------------------------------------------------- |
|---|
| 1076 | ! Parameters: |
|---|
| 1077 | CHARACTER(LEN=*), INTENT(IN) :: str |
|---|
| 1078 | !------------------------------------------------------------------------------- |
|---|
| 1079 | ! Local variable: |
|---|
| 1080 | INTEGER :: k, ix |
|---|
| 1081 | !------------------------------------------------------------------------------- |
|---|
| 1082 | out=str |
|---|
| 1083 | DO k=1,LEN(str) |
|---|
| 1084 | ix=IACHAR(str(k:k)); IF(64<ix.AND.ix<91) out(k:k)=ACHAR(ix+32) |
|---|
| 1085 | END DO |
|---|
| 1086 | |
|---|
| 1087 | END FUNCTION strLow |
|---|
| 1088 | ! |
|---|
| 1089 | !------------------------------------------------------------------------------- |
|---|
| 1090 | |
|---|
| 1091 | END PROGRAM make_sso |
|---|
| 1092 | ! |
|---|
| 1093 | !------------------------------------------------------------------------------- |
|---|
