program msshcorr c computes a long wavelength interannual correctional term c to the mean sea surface dependent on the period specified. c c mssh is computed over the period 01-01-93 - 32-12-04 (12 years). c c It is optional to correct for trend + interannual mean values. c It is optional to correcdt of intra-annual variations c The trend is computed in a joint estimation with the mean value. c The interannual terms are computed as the annual mean value c where the mean + trend + estimated annual signal has been removed. c The annual signal have been estiamated using the following proce c dure to take into account intra-annual variation. c January - December monthly means have been calculated in 2 by 2 c degree bins and the monthly mean is subsequently removed from c the observation. c c The corrections are computed using T/P altimetry, consequently c they are valid within the +/- 66 parallel. c By specifying a period (month,year ) and select terms these c will be estimated and added to the mean sea surface. c the KMS03 mean sea surface have computed using no-ib corrected c data because we. c An optional correction for the mean value of the inverse barometer c correction (as computed using the NASA Pathfinder data) c can also be applied c OA - 08-2003 c-------------------------- c The 27 grids are stored in one common gridfile. c Grid 1 = trend (cm/year) c Grid 2 = 9y mean value of Inverse barometer correction (<10 cm) c Grid 3 = 1993 annual mean deviation in ssh (<20 cm) c Grid 4 = 1994 annual mean deviation in ssh (<30 cm) c Grid 5 = 1995 annual mean deviation in ssh (<30 cm) c Grid 6 = 1996 annual mean deviation in ssh (<30 cm) c Grid 7 = 1997 annual mean deviation in ssh (<30 cm) c Grid 8 = 1998 annual mean deviation in ssh (<30 cm) c Grid 9 = 1999 annual mean deviation in ssh (<20 cm) c Grid 10 = 2000 annual mean deviation in ssh (<20 cm) c Grid 11 = 2001 annual mean deviation in ssh (<20 cm) c Grid 12 = 2002 annual mean deviation in ssh (<20 cm) c Grid 13 = 2003 annual mean deviation in ssh (<20 cm) c Grid 14 = 2004 annual mean deviation in ssh (<20 cm) c Grid 15-26 = Monthly mean deviation in ssh (<20 cm) c Grid 27 = Bathymetry to screen of land (>50 meters) c c--------------------------- implicit real*8 (a-h,o-z) real*8 acor(141,361,27) real*4 fac(27) character*30 ofile1,ofile2 character*1 ch logical lcasp open (20,file = 'msshcorr.grd',form = 'formatted',err=9004) open (21,file = 'msshspec.grd',form = 'formatted',err=9004) ip = 141 jp = 361 kp = 27 c number of years used for mssh computation iyear = 12 do k = 1,kp fac(k) = 0.0 read(20,*) alatmin,alatmax,alonmin,alonmax,dlat,dlon read(20,*) ((acor(i,j,k),j=1,jp),i = 1,ip) end do write(21,*) alatmin,alatmax,alonmin,alonmax,dlat,dlon write(*,*) 'Program computes the average of the sea level' write(*,*) 'variability relative to the DNSC08MSSH ' write(*,*) 'for a user-defined time epoch' write(*,*) '(must me subperiod to 1993-2004 used for DNSC08)' write(*,*) 'Optionally including the following signals:' write(*,*) ' Sea level trend' write(*,*) ' Intra-annual variation (<1y)' write(*,*) ' Inter-annual variation (>1.5 years)' Write(*,*) ' Inv. barometer effect (to create NoIB MSS)' Write(*,*) '-----------------------------------------------' write(*,*) 'Computation is only long wavelengt' write(*,*) 'using data from all 12 years' write(*,*) 'Some interpolation residuals is therefor likely' Write(*,*) '-----------------------------------------------' Write(*,*) 'Initially specify period to adjust variability to' 9001 write(*,*) 'Enter start time (month, year)' write(*,*) 'Valid years are 1993 -> 2004) ' write(*,*) 'Valud month are 01 -> 13 ' write(*,*) '(01 is 1.Jan, 12 is 1.Dec, 13 is 31.Dec) ' read (*,*) sm,sy if (sm.gt.13) sm = 13 if (sm.lt.1) sm = 1 if (sy.gt.2004) sy = 2004 if (sy.lt.1993) sy = 1993 write(*,*) 'Enter end time (month, year)' write(*,*) 'Valid years are 1993 -> 2004) ' write(*,*) 'Valud month are 01-13 ' write(*,*) '(01 is 1. Jan, 12 is 1. Dec, 13 is 31.Dec) ' read (*,*) em,ey if (em.gt.13) em = 13 if (em.lt.1) em = 1 if (ey.gt.2004) ey = 2004 if (ey.lt.1993) ey = 1993 write(*,*) 'Start time ', sm, ' ', sy write(*,*) 'End time ', em, ' ', ey write(*,*) 'Correct (y/n)' read(*,'(A)') ch if (ch.eq.'n'.or.ch.eq.'N') goto 9001 Write(*,*) '------------------------------------------' Write(*,*) 'Specify signals to include in computation' write(*,*) 'Effect of 12y linear sea level trend signal' write(*,*) 'over specified period (y/n)' read(*,'(A)') ch if (ch.eq.'y'.or.ch.eq.'Y') then sp = sy+(sm-1)/12.0 ep = ey+(em-1)/12.0 fac(1) = (sp+ep)/2.0-1999.0 endif write(*,*) fac(1) write(*,*) 'Inverse Barometer effect(for No-ib.MSSH(y/n)' read(*,'(A)') ch if (ch.eq.'y'.or.ch.eq.'Y') fac(2) = 100.0 write(*,*) 'Inter-annual signal' write(*,*) '(Deviation from annual mean)' write(*,*) 'over specified period (y/n)' read(*,'(A)') ch if (ch.eq.'y'.or.ch.eq.'Y') then write(*,*) sy,ey do k = sy,ey fac(k-1990) = 1.0 enddo fac(sy-1990) = (13-sm)/12.0 fac(ey-1990) = (em-1)/12.0 endif write(*,*) 'Intra-annual signal (mean monthly deviations)' write(*,*) '(over specified period)(y/n)' write(*,*) '(only if period includes submonthly signal)' read(*,'(A)') ch if (ch.eq.'y'.or.ch.eq.'Y') then write(*,*) sm,em if (em.gt.sm) then do k = sm,em-1 fac(iyear+2+k) = 1.0 enddo endif if (em.lt.sm) then do k = em,12 fac(iyear+2+k) = 1.0 end do do k = 1,sm-1 fac(iyear+2+k) = 1.0 end do endif endif write(*,*) (fac(k),k=1,2) write(*,*) (fac(k),k=3,14) write(*,*) (fac(k),k=15,26) cmean = 0.0 cstd = 0.0 cmax = -1000000.0 cmin = 1000000.0 do i = 1,ip do j = 1,jp lcasp = .false. if (i.lt.36.and.j.gt.30.and.j.lt.70) lcasp=.true. cout = 0.0 c Here a 5 meter depth mask is enforced for the computation if (acor(i,j,27).lt.-5.and..not.lcasp) then inum = inum + 1 do k = 1,kp cout = cout + fac(k)*acor(i,j,k) end do if (cout.gt.cmax) cmax = cout if (cout.lt.cmin) cmin = cout cmean = cmean + cout cstd = cstd + cout*cout endif write(21,*) cout enddo enddo cmean = cmean/inum cstd = sqrt((cstd/(inum-1))-cmean*cmean) write(*,*) 'Correction: Minumum, maximum, mean, std (cm)' write(*,*) cmin,cmax,cmean,cstd close(20) close(21) goto 9005 9004 write(*,*) ' Error opening input file' 9005 continue end