module interp use const use comm use var !TS use io implicit real*8 (a-h,o-z) # include "pnetcdf.inc" !dog contains subroutine inmets ! ====================================================================== ! extract regional subset of etopo5 bathymetry ! user specified lat-long window and resolution ! by david dietrich, may 1, 1995 ! ! reads one latitude at a time in a simple do loop. ! reads full etopo5 data sets into memory, then extracts desired subset ! simple and versatile. ! most workstations have the required storage capacity. ! etopo5 has 360 x 12 integer*2 data elements per (latitudinal) record. ! this gives 8640 bytes or 2160 sgi (or vax) "words". thus, under sun ! fortran, recl=8640, while under sgi or vax fortran recl=2160. ! ====================================================================== !nres is the etopo resolution =5 (5min) parameter(nres=2) parameter(nx=21600/nres,ny=10800/nres,nxpd=nx/360,nypd=ny/180) integer*2 meters,in ! ======================================================== ! tlev, slev is levitus climatology data with landfills ! t,s is tlev, slev interpolated to model horizontal grid ! t1, s1 is t,s vertically interpolated to model z-levels ! ======================================================== dimension rmets(i0,j0),meters(nx,ny),in(i0,j0) dimension depth(nx,ny) dimension istart(2),icount(2),dt(nx,ny) ! read nc data, which is single precision real,dimension(nx,ny) :: dep_tmp ! ====================================================================== ! i0= number of grid points in the x-direction including ghost zones ! j0= number of grid points in the y-direction including ghost zones ! westdeg is western most longitude degrees (west cell face of i=2 zone) ! longitude coordinate increases eastward ! 0 < westdeg < 360 ! longitudinal increment in minutes,dxmnut, is read from yzgrid ! latitudes are read from yzgrid ! ====================================================================== ! open(9,file='data/westdeg'//grd) ! rewind 9 ! read(9,'(f9.4)') westdeg ! close(9) call rdnc_westdeg !dog ! ================ ! input data files ! ================ ! open(10,file='data/yzgrid'//grd,form='unformatted') ! new, single record integer*2 data ! no record length bullshit! if (myid .eq. 0) then if(nres==1)then istatus=nf_open('../worldata/etopo/etopo1/etopo1_bed_c_gmt4.grd',nf_nowrite,ncid) print*,'istatus1',istatus,ncid,nx,ny elseif(nres==2) then istatus=nf_open('../worldata/etopo/etopo2/etopo2v2g_f4.nc',nf_nowrite,ncid) print*,'istatus2',istatus,ncid,nx,ny,NF_STRERROR(istatus) elseif(nres==5) then open(12,file='../worldata/etopo/etopo5/itopo5',form='unformatted') endif endif ! rewind 10 rewind 12 rewind 35 ! ================= ! output data files ! ================= ! open(8,file='data/depth'//grd,form='unformatted') ! open(14,file='data/inview'//grd) ! rewind 8 ! rewind 14 ! ============================================== ! read metrics data file yzgrid (from metgen.f) ! ============================================== call rdnc_yzgrid(j,k,dxmnut) !dog ! read(10) j,k ! if (j0.eq.j.and.k0.eq.k) go to 9 if (j0.ne.j.and.k0.ne.k) then write(*,6) 6 format('stop due to resolution mismatch. set j0,k0 to values from metgen.f.') stop endif ! 9 read(10) tlz,z,odz,odzw,dxmnut,xdeg,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv ! new single record input if (myid .eq. 0) then if(nres.le.2)then istart(1) = 1 istart(2) = 1 icount(1) = nx icount(2) = ny istatus=nf_get_vara_real(ncid,3,istart,icount,dep_tmp) depth = dep_tmp do j=1,ny do i=1,nx/2 depth(i,j)=depth(i,j)+depth(i+nx/2,j) depth(i+nx/2,j)=depth(i,j)-depth(i+nx/2,j) depth(i,j)=depth(i,j)-depth(i+nx/2,j) enddo enddo print*,'icount,',icount do j=1,ny do i=1,nx dt(i,j)=depth(i,ny-j+1) enddo enddo do j=1,ny do i=1,nx depth(i,j)=dt(i,j) enddo enddo !ts will add this later else if(nres.eq.5) then read(12) meters do j=1,ny do i=1,nx depth(i,j)=meters(i,j) enddo enddo endif endif call mpi_bcast(depth,nx*ny,mpi_real8,0,mpi_comm_world,ierr) ! ============================================================== ! determine depths at model grid points and store in rmets array ! ============================================================== ! first element is at 0.0 degrees east (or 360 degrees east). ! if(nres.eq.1)then xinc=dxmnut/1. else if(nres.eq.2)then xinc=dxmnut/2. else if(nres.eq.5)then xinc=dxmnut/5. endif if(myid.eq.0) print*,'xinc',xinc ! xin! is the grid distance of etop2. therefore, for 15' grid size ! each grid has 15/2 = 7.5 grid point in etop2 data set. do 100 jj=1,j0 ! ylat is measured from north pole. first point is at pole. ! xlon is measured units 1/12 deg east, from grenwich, with first point ! at 1/12 deg east (meters(1,j) is at 1/12 deg east; meters(i,1) is ! at north pole, or ydeg=90) ylat=nypd*(90-ydeg(jj))+1 ! average depths in lat-long window j=ylat jp=j-1 xlon=nxpd*westdeg-.5*xinc ! print*,'xinc',xinc,xlon,ydeg(524) do 100 ii=2,i1 xlon=xlon+xinc i=xlon tmp1=xlon-i tmp2=1.-tmp1 ip=mod(i+1,nx)+1 i=mod(i,nx)+1 !choke point depths are critical in modeling straits. when they are not ! adequately resolved they are always too shallow when evaluating them from ! world data base. thus, at a given control volume location, it is ! physically better to assign the deepest adjacent point in the control ! volume region than to use the actual local value. ! special for 1/4 deg global grid in critical gom and gibraltar regions if (xdeg(ii).lt.269.625.or.xdeg(ii).gt.359.125) go to 80 ! avoid cape hatteras and northern north atlanti! region if (xdeg(ii).lt.299.625.and.ydeg(jj).gt.29.91505) go to 80 ! this may make some single-point islands wet, but that's ok! ! note: negative meters values mean positive ocean depths rmets(ii,jj)=min(depth(i,j),depth(ip,j),depth(i,jp),depth(ip,jp)) rmets(ii,jj)=-rmets(ii,jj) rmets(ii,jj)=max(0.,rmets(ii,jj)) go to 100 80 dsw=depth(i,jp) dse=depth(ip,jp) dnw=depth(i,j) dne=depth(ip,j) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw rmets(ii,jj)=(ylat-j)*ds+(1-ylat+j)*dn rmets(ii,jj)=-rmets(ii,jj) rmets(ii,jj)=max(0.,rmets(ii,jj)) 100 continue ! open(333,file='data/rmets'//grd,form='unformatted') ! write(333) rmets ! close(333) ! write(8) rmets call wr8nc_depth(rmets) !dog m=0 n=0 tmp=0. temp=.01*tlz do j=2,j1 do i=2,i1 m=m+in(i,j) tmp=max(tmp,rmets(i,j)) if (rmets(i,j).gt.temp) n=n+1 enddo enddo write(*,103) m,n,temp,tmp 103 format(i8,' total wet points,',i6,' deeper than ',f5.0,' m, deepest= ',f6.0,' m') end subroutine subroutine indata ! ====================================================================== ! read bathymetry file from inmets.f ! read salinity, temperature and wind data sets ! interpolate data to model grid ! ====================================================================== ! lscowf: .true. scow wind climatological forcing ! lscowf: .false. hr wind climatological forcing ! lclim: 1 levitus 94 climatology ! lclim: 2 phc3 climatology parameter(lscowf=.true.,lclim=2) parameter(nxscow=1440,nyscow=560,nxhel=360,nyhel=180) integer*2 nzdat ! real*4 slon,slat,ilon,ilat,tlevs,slevs ! real*4 tmp variables for file reading real*4,dimension(360,180,33) :: swoa_tmp,twoa_tmp real*4,dimension(360,180,33) :: tlevs_tmp,slevs_tmp real*4,dimension(nxhel,nyhel) :: tx_tmp,ty_tmp real*4,dimension(nxscow,nyscow) :: sctx_tmp,scty_tmp ! ======================================================== ! tlev, slev is levitus climatology data with landfills ! t,s is tlev, slev interpolated to model horizontal grid ! t1_data, s1_data is t_data,s_data vertically interpolated to model z-levels ! ======================================================== dimension tlev(360,180,33),slev(360,180,33),rho_data(360,180,33), & t_data(i0,j0,33),s_data(i0,j0,33),t1_data(i0,j0,k1),s1_data(i0,j0,k1),taux_data(i2,j2), & tauy_data(i2,j2),tauxs(i2,j2),tauys(i2,j2),nzdat(33) dimension tlevs(360,180,33),slevs(360,180,33) dimension twoa(360,180,33),swoa(360,180,33) dimension twoal(0:361,0:181),swoal(0:361,0:181) dimension tx_data(nxhel,nyhel),ty_data(nxhel,nyhel) dimension scowtx_data(nxscow,nyscow),scowty_data(nxscow,nyscow),istart(2),icount(2) character*59 txc(12) character*59 tyc(12) character*59 season(4) character*60 twoaf(12),swoaf(12),seasonswoa(4),seasontwoa(4) character*80 scowtaux,scowtauy integer*2 maxj(i2,j2) dimension scr_data(i0,j0,2) real*8,dimension(i0,j0,k1,12) :: clidata_t1,clidata_s1 real*8,dimension(i2,j2,12) :: cli_tauxs,cli_tauys real*8,dimension(i2,j2,12) :: clidata_tauy,clidata_taux data txc/ & '../worldata/hellerman/tx-jan.dat', & '../worldata/hellerman/tx-feb.dat', & '../worldata/hellerman/tx-mar.dat', & '../worldata/hellerman/tx-apr.dat', & '../worldata/hellerman/tx-may.dat', & '../worldata/hellerman/tx-jun.dat', & '../worldata/hellerman/tx-jul.dat', & '../worldata/hellerman/tx-aug.dat', & '../worldata/hellerman/tx-sep.dat', & '../worldata/hellerman/tx-oct.dat', & '../worldata/hellerman/tx-nov.dat', & '../worldata/hellerman/tx-dec.dat'/ data tyc/ & '../worldata/hellerman/ty-jan.dat', & '../worldata/hellerman/ty-feb.dat', & '../worldata/hellerman/ty-mar.dat', & '../worldata/hellerman/ty-apr.dat', & '../worldata/hellerman/ty-may.dat', & '../worldata/hellerman/ty-jun.dat', & '../worldata/hellerman/ty-jul.dat', & '../worldata/hellerman/ty-aug.dat', & '../worldata/hellerman/ty-sep.dat', & '../worldata/hellerman/ty-oct.dat', & '../worldata/hellerman/ty-nov.dat', & '../worldata/hellerman/ty-dec.dat'/ data season/ & '../worldata/levitus/modelev_win.dat', & '../worldata/levitus/modelev_spr.dat', & '../worldata/levitus/modelev_sum.dat', & '../worldata/levitus/modelev_aut.dat'/ data swoaf/ & '../worldata/phc3/salt01_p3', & '../worldata/phc3/salt02_p3', & '../worldata/phc3/salt03_p3', & '../worldata/phc3/salt04_p3', & '../worldata/phc3/salt05_p3', & '../worldata/phc3/salt06_p3', & '../worldata/phc3/salt07_p3', & '../worldata/phc3/salt08_p3', & '../worldata/phc3/salt09_p3', & '../worldata/phc3/salt10_p3', & '../worldata/phc3/salt11_p3', & '../worldata/phc3/salt12_p3'/ data twoaf/ & '../worldata/phc3/temp01_p3', & '../worldata/phc3/temp02_p3', & '../worldata/phc3/temp03_p3', & '../worldata/phc3/temp04_p3', & '../worldata/phc3/temp05_p3', & '../worldata/phc3/temp06_p3', & '../worldata/phc3/temp07_p3', & '../worldata/phc3/temp08_p3', & '../worldata/phc3/temp09_p3', & '../worldata/phc3/temp10_p3', & '../worldata/phc3/temp11_p3', & '../worldata/phc3/temp12_p3'/ data seasonswoa/ & '../worldata/phc3/salt13_p3', & '../worldata/phc3/salt15_p3', & '../worldata/phc3/salt15_p3', & '../worldata/phc3/salt13_p3'/ data seasontwoa/ & '../worldata/phc3/temp13_p3', & '../worldata/phc3/temp15_p3', & '../worldata/phc3/temp15_p3', & '../worldata/phc3/temp13_p3'/ data scowtauy/ & '../worldata/scow/ws_v.nc'/ data scowtaux/ & '../worldata/scow/ws_u.nc'/ ! ------------------------------------- ! standard levitus data depths (m) ! ------------------------------------- data nzdat/ 0, 10, 20, 30, 50, 75, 100, 125, 150, 200, 250, & 300, 400, 500, 600, 700, 800, 900,1000,1100,1200,1300,1400,1500, & 1750,2000,2500,3000,3500,4000,4500,5000,5500/ !lat: /lat (degree_north) ordered (89.5s) to (89.5n) by 1.0 n= 180 pts :grid !lon: /lon (degree_east) periodi! (0.5e) to (0.5w) by 1.0 n= 360 pts :grid ! ====================================================================== ! i0= number of grid points in the x-direction including ghost zones ! j0= number of grid points in the y-direction including ghost zones ! westdeg is western most longitude degrees (west cell face of i=2 zone) ! longitude coordinate increases eastward ! 0 < westdeg < 360 ! longitudinal increment in minutes,dxmnut, is read from yzgrid ! latitudes are read from yzgrid ! ====================================================================== ! ================ ! input data files ! ================ ! open(10,file='data/yzgrid'//grd,form='unformatted') ! open(14,file='data/invu'//grd) ! rewind 10 ! rewind 14 ! ================= ! output data files ! ================= if (myid .eq. 0) open(41,file='data/fixedlev',form='unformatted') ! annual cycle levitus climatology on model grid ! open(70,file='data/annualevitus'//grd,form='unformatted') ! s,t at surface and lateral boundaries ! open(72,file='data/winds'//grd,form='unformatted') rewind 41 ! rewind 70 ! rewind 72 ! ============================================== ! read metrics data file yzgrid (from metrics.f) ! ============================================== call rdnc_yzgrid(j,k,dxmnut) !dog ! read(10) j,k if (j0.ne.j.and.k0.ne.k) then write(*,6) 6 format('stop due to resolution mismatch. set j0,k0 to values from metgen.f.') stop endif ! 9 read(10) tlz,z,odz,odzw,dxmnut,xdeg,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv ! scow wind forcing if (myid.eq.0.and.lscowf.eq..true.) then ist=nf_open(scowtaux,nf_nowrite,nctx) ist=nf_open(scowtauy,nf_nowrite,ncty) istart(1)=1 istart(2)=1 icount(1)=nxscow icount(2)=nyscow endif ! ================================================== ! prepare monthly winds and seasonal s,t climatology ! ================================================== do m=1,12 if (myid .eq. 0) then rewind 41 write(*,111) m 111 format(i3) ! hellerman open(15,file=txc(m)) open(16,file=tyc(m)) rewind 15 rewind 16 if (mod(m,3).eq.1) then nm=(m+2)/3 ! =================================================== ! read modelfied (with landfills) levitus climatology ! and fix negative stratification areas ! =================================================== open(36,file=season(nm),form='unformatted') rewind 36 if(lclim.ne.1)then open(37,file=seasonswoa(nm)) open(38,file=seasontwoa(nm)) rewind 37 rewind 38 endif ! =================================================== ! read modelfied (with landfills) levitus climatology ! and fix negative stratification areas ! =================================================== ! read(40) ! 1 ((tlev(i,j,1),i=1,360),j=1,180),((slev(i,j,1),i=1,360),j=1,180) print*,'read seasonal data' do k=1,33 read(36) ((tlevs_tmp(i,j,k),i=1,360),j=1,180),((slevs_tmp(i,j,k),i=1,360),j=1,180) enddo tlevs = tlevs_tmp slevs = slevs_tmp if(lclim.ne.1)then do k=1,33 do j=1,180 do i=1,36 read(37,'(10f8.4)') (swoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10) read(38,'(10f8.4)') (twoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10) enddo enddo enddo swoa = swoa_tmp twoa = twoa_tmp endif endif ! read(40) ! 1 ((tlev(i,j,1),i=1,360),j=1,180),((slev(i,j,1),i=1,360),j=1,180) if(lclim.ne.1)then open(39,file=swoaf(m)) open(40,file=twoaf(m)) rewind 39 rewind 40 ! =================================================== ! read modelfied (with landfills) levitus climatology ! and fix negative stratification areas ! =================================================== do 124 k=1,24 do 124 j=1,180 do 124 i=1,36 read(40,'(10f8.4)') (twoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10) 124 read(39,'(10f8.4)') (swoa_tmp(ii,j,k),ii=(i-1)*10+1,i*10) swoa = swoa_tmp twoa = twoa_tmp endif ! print*,'month twoa',(twoa(180,j,1),j=1,180) 127 continue if(lclim.eq.1)then do k=1,33 do j=1,180 do i=1,360 swoa(i,j,k)=slevs(i,j,k) twoa(i,j,k)=tlevs(i,j,k) enddo enddo enddo else do 128 k=1,33 do 128 j=1,180 do 128 i=1,360 if (swoa(i,j,k) .le. -10.) swoa(i,j,k)=slevs(i,j,k) 128 if (twoa(i,j,k) .le. -10.) twoa(i,j,k)=tlevs(i,j,k) endif write(41) ((twoa(i,j,1),i=1,360),j=1,180),((swoa(i,j,1),i=1,360),j=1,180) ! nzdat are levitus levels zlev = dble(nzdat(1)) do 142 j=1,180 do 142 i=1,360 if ((twoa(i,j,1) .lt. -10.0)) twoa(i,j,1) = 0.0 if ((swoa(i,j,1) .lt. -10.0)) swoa(i,j,1) = 0.0 twoa(i,j,1)=theta(swoa(i,j,1),twoa(i,j,1),zlev,dble(0.)) 142 rho_data(i,j,1)=r(twoa(i,j,1),swoa(i,j,1),zlev) ! eliminate unstable stratification data ssum=0. n=0 lb=1 lt=2 do 150 k=2,33 ! read(40) ! 1 ((tlev(i,j,lt),i=1,360),j=1,180),((slev(i,j,lt),i=1,360),j=1,180) ! for density referenced to in situ pressure zlev = dble(nzdat(k)) ! for density referenced to mean pressure level ! zlev = float(nzdat(k-1)+nzdat(k))/2.0 ! write(*,*) zlev do 145 j=1,180 do 145 i=1,360 if ((twoa(i,j,k) .lt. -10.0)) twoa(i,j,k) = 0.0 if ((swoa(i,j,k) .lt. -10.0)) swoa(i,j,k) = 0.0 twoa(i,j,k)=theta(swoa(i,j,k),twoa(i,j,k),zlev,dble(0.)) 145 rho_data(i,j,k)=r(twoa(i,j,k),swoa(i,j,k),zlev) do 148 j=1,180 do 148 i=1,360 rhok=rho_data(i,j,k) rhokm=rho_data(i,j,k-1) if (rhok.ge.rhokm) go to 148 n=n+1 twoa(i,j,k)=twoa(i,j,k-1) swoa(i,j,k)=swoa(i,j,k-1) rho_data(i,j,k)=rhokm 148 continue 150 write(41) ((twoa(i,j,k),i=1,360),j=1,180),((swoa(i,j,k),i=1,360),j=1,180) insum=32*180*360 ! ================================================= ! interpolate levitus data to model horizontal grid ! ================================================= close(40) endif rewind 41 do 230 k=1,33 if (myid .eq. 0) then read(41) ((tlev(i,j,1),i=1,360),j=1,180),((slev(i,j,1),i=1,360),j=1,180) endif call mpi_bcast(tlev,360*180,mpi_real8,0,mpi_comm_world,ierr) call mpi_bcast(slev,360*180,mpi_real8,0,mpi_comm_world,ierr) ! if(k.eq.1.or.k.eq.2) then ! write(*,'(i2,a6,10f10.7)'),myid,' slev=',slev(340,161:170,1) ! write(*,'(i2,a6,10f10.7)'),myid,' slev=',slev(320,101:110,1) ! endif slon=.5 slat=-89.5 ilon=1. ilat=1. do 230 jj=1,j0 ! xlon is model data point, measured in degrees east of grenwich ! ylat is model data point, measured in degrees north of south pole ! ydeg is measured from -90 at south pole to 90 at north pole ! tlev(1,1) is at 0.5 deg e, -89.5 deg lat ! tlev(360,180,) is at 359.5 deg e, 89.5 deg lat ! allowed model latitudes: -89.5.lt.ydeg.lt.89.5 ! allowed model longitudes: 0.le.xlon.le.720 ! ylat=90.+ydeg(jj) ! j=ylat+0.5 ylat=(ydeg(jj)-slat)*ilat+1. j=ylat jp=j+1 do 230 ii=1,i0 ! xlon=xdeg(ii) xlon=(xdeg(ii)-slon)*ilon+1. if (xlon.gt.360) xlon=xlon-360. ! xplus is model data point distance east of 0.5 deg west(!) ! note that xplus is always positive due to xlon restrictions ! xplus=xlon+0.5 xplus=xlon i=xplus ip=i+1 ! tmp1 is weighting of east point tmp1=xplus-dble(i) tmp2=1.-tmp1 ! special treatments when 359.5 < xlon < 360 and 0 < xlon < 0.5 if (ip.eq.361) ip=1 ! if (i.eq.0) i=359 if (i.eq.0) i=360 dnw=tlev(i,jp,1) dne=tlev(ip,jp,1) dsw=tlev(i,j,1) dse=tlev(ip,j,1) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw t_data(ii,jj,k)=(ylat-j)*dn+(jp-ylat)*ds dnw=slev(i,jp,1) dne=slev(ip,jp,1) dsw=slev(i,j,1) dse=slev(ip,j,1) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw 230 s_data(ii,jj,k)=(ylat-j)*dn+(jp-ylat)*ds tmax=-100. tmin=100. do 240 j=2,j1 do 240 i=2,i1 tmax=max(tmax,t_data(i,j,1)) 240 tmin=min(tmin,t_data(i,j,1)) ! write(*,245) m,tmin,tmax ! 245 format('month',i3,' tmin,tmax=',2f6.2) ! =============================================== ! interpolate from levitus levels to model levels ! =============================================== z2=100.*nzdat(1) k=1 zc=z(2) n=1 320 continue if(n.lt.33) then n=n+1 z1=z2 z2=100.*nzdat(n) endif if (z2.lt.zc.and.n.lt.33) go to 320 322 cc1=(z2-zc)/(z2-z1) cc2=(zc-z1)/(z2-z1) ! write(14,323) k,n ! 323 format('for diecast model level',i3,', lower input level=',i3) do 325 j=1,j0 do 325 i=1,i0 t1_data(i,j,k)=cc1*t_data(i,j,n-1)+cc2*t_data(i,j,n) 325 s1_data(i,j,k)=cc1*s_data(i,j,n-1)+cc2*s_data(i,j,n) k=k+1 if (k.eq.k0) go to 390 zc=z(2*k) if (zc.gt.z2) go to 320 go to 322 ! write all seasons or months so we can diagnose annual cycle ! 390 write(70) s1_data,t1_data 390 clidata_t1(:,:,:,m) = t1_data clidata_s1(:,:,:,m) = s1_data ! wenien ! open(777,file='/home/opt/s1t1'//grd,form='unformatted') ! write(777) s1,t1,s,t,swoa,twoa ! close(777) ! call mpi_finalize(ierr) ! stop 390 enddo ! increment westdeg by one grid interval, and start at second latitude ! (no ghost zones in winds) do m=1,12 if (myid .eq. 0) then open(15,file=txc(m)) open(16,file=tyc(m)) read(15,400) i,j,x_tmp,x_tmp,x_tmp,x_tmp,fmt read(16,400) i,j,x_tmp,x_tmp,x_tmp,x_tmp,fmt 400 format(2i4/4(3x,f10.7),a10) ! ================ ! read hellerman winds ! ================ read(15,415) tx_tmp read(16,415) ty_tmp tx_data = tx_tmp ty_data = ty_tmp 415 format(5e15.7) write(*,416) 416 format('read wind stress') endif call mpi_bcast(tx_data,nxhel*nyhel,mpi_real8,0,mpi_comm_world,ierr) call mpi_bcast(ty_data,nxhel*nyhel,mpi_real8,0,mpi_comm_world,ierr) ! xlon is model data point, measured in degrees east of grenwich ! ylat is model data point, measured in degrees north of south pole ! ydeg is measured from -90 at south pole to 90 at north pole ! tlev(1,1) is at 0.5 deg e, -89.5 deg lat ! tlev(360,180,) is at 359.5 deg e, 89.5 deg lat ! allowed model latitudes: -89.5.lt.ydeg.lt.89.5 ! allowed model longitudes: 0.le.xlon.le.720 slon=.5 slat=-89.5 ilon=1. ilat=1. do jj=2,j1 ! ylat=90.+ydeg(jj) ! j=ylat+0.5 ylat=(ydeg(jj)-slat)*ilat+1. j=ylat jp=j+1 do ii=2,i1 ! xlon=xdeg(ii+1) xlon=(xdeg(ii)-slon)*ilon+1. if (xlon.gt.dble(nxhel)) xlon=xlon-dble(nxhel) ! xplus=xlon+0.5 xplus=xlon i=xplus ip=i+1 tmp1=xplus-dble(i) tmp2=1.-tmp1 if (ip.eq.(nxhel+1)) ip=1 ! if (i.eq.0) i=nxhel-1 if (i.eq.0) i=nxhel dnw=tx_data(i,jp) dne=tx_data(ip,jp) dsw=tx_data(i,j) dse=tx_data(ip,j) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw ! taux(ii-1,jj-1)=(.5+ylat-j)*dn+(.5-ylat+j)*ds taux_data(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds dnw=ty_data(i,jp) dne=ty_data(ip,jp) dsw=ty_data(i,j) dse=ty_data(ip,j) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw tauy_data(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds enddo enddo ! if(lscowf)then if (myid .eq. 0) then ! hellerman ierr=nf_get_vara_real(nctx,m+2,istart,icount,sctx_tmp) ierr=nf_get_vara_real(ncty,m+2,istart,icount,scty_tmp) scowtx_data = sctx_tmp scowty_data = scty_tmp !dog ! print*,nf_strerror(ierr) ! write(1233)tx,ty do 420 j=1,nyscow do 420 i=1,nxscow if(abs(scowtx_data(i,j)).lt.10..and.abs(scowty_data(i,j)).lt.10.)go to 420 scowtx_data(i,j)=0. scowty_data(i,j)=0. 420 continue endif ! stop call mpi_bcast(scowtx_data,nxscow*nyscow,mpi_real8,0,mpi_comm_world,ierr) call mpi_bcast(scowty_data,nxscow*nyscow,mpi_real8,0,mpi_comm_world,ierr) ! =================================================== ! interpolate scow data to model horizontal grid ! =================================================== slon=.125 slat=-69.875 ilon=4. ilat=4. do jj=2,j1 ! ============================================================ ! tx_data(i,1) is at latitude -69.875; tx_data(i,560) is at latitude 69.875 ! tx_data(1,j) is at longitude 0.125; tx_data(1440,j) is at longitude 359.875 ! xlon is measured in degrees east from 0. (xlon=0) ! ydeg is measured in degrees north from -69.875 (ydeg=0) ! ============================================================ ! xlon is measured in degrees east from 0 at grenwich ! ydeg is measured from -90 at south pole to 90 at north pole ylat=(min(ydeg(jj),69.875))*ilat+1. ! ylat=(ydeg(jj)-slat)*ilat+1. ! ylat=(70.125+ydeg(jj))*4. j=ylat if(j.le.0) j=1 jp=j+1 do ii=2,i1 xlon=(xdeg(ii)-slon)*ilon+1. ! xlon=xdeg(ii)*4.+0.5 ! xlon=xdeg(ii)*4. if (xlon.gt.dble(nxscow)) xlon=xlon-dble(nxscow) ! xplus is model data point distance east of 0.5 deg west(!) ! note that xplus is always positive due to xlon restrictions xplus=xlon i=xplus ip=i+1 ! tmp1 is weighting of east point tmp1=xplus-i tmp2=1.-tmp1 ! special treatments when 359.5 < xlon < 360 and 0 < xlon < 0.5 if (ip.eq.(nxscow+1)) ip=1 ! if (i.eq.0) i=nxscow-1 if (i.eq.0) i=nxscow dnw=scowtx_data(i,jp) dne=scowtx_data(ip,jp) dsw=scowtx_data(i,j) dse=scowtx_data(ip,j) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw tauxs(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds dnw=scowty_data(i,jp) dne=scowty_data(ip,jp) dsw=scowty_data(i,j) dse=scowty_data(ip,j) dn=tmp1*dne+tmp2*dnw ds=tmp1*dse+tmp2*dsw tauys(ii-1,jj-1)=(ylat-j)*dn+(jp-ylat)*ds enddo enddo endif if(lscowf)then do i=1,i2 do j=1,j2 maxj(i,j)=0 if(abs(tauxs(i,j)*tauys(i,j)).le.1.e-4)then iss=i-1 iff=i+1 if(iss.eq.0) iss=1 if(iff.eq.i1) iff=i2 jss=j-1 jff=j+1 if(jss.eq.0) jss=1 if(jff.eq.j1) jff=j2 do ii=iss,iff do jj=jss,jff maxj(ii,jj)=1 enddo enddo tauxs(i,j)=taux_data(i,j) tauys(i,j)=tauy_data(i,j) endif enddo enddo ! lf smoothing !dog ??? call mpi_type_vector(2*j0,1,i0,mpi_real8,m_vlon3d,ierr) call mpi_type_vector(2,i0,i0*j0,mpi_real8,m_vlat3d,ierr) call mpi_type_commit(m_vlon3d,ierr) call mpi_type_commit(m_vlat3d,ierr) js=2 jf=j1 if(mylat.eq.0) js=3 if(mylat.eq.ngy1) jf=j2 mxpas=10 do n=1,mxpas do i=2,i1 do j=2,j1 scr_data(i,j,1)=tauxs(i-1,j-1) scr_data(i,j,2)=tauys(i-1,j-1) enddo enddo call mpi_sendrecv(scr_data(2,1,1),1,m_vlon3d,m_w,1,scr_data(i0,1,1),1,m_vlon3d,m_e,1,mpi_comm_world,istat,ierr) call mpi_sendrecv(scr_data(i1,1,1),1,m_vlon3d,m_e,1,scr_data(1,1,1),1,m_vlon3d,m_w,1,mpi_comm_world,istat,ierr) call mpi_sendrecv(scr_data(1,2,1),1,m_vlat3d,m_s,1,scr_data(1,j0,1),1,m_vlat3d,m_n,1,mpi_comm_world,istat,ierr) call mpi_sendrecv(scr_data(1,j1,1),1,m_vlat3d,m_n,1,scr_data(1,1,1),1,m_vlat3d,m_s,1,mpi_comm_world,istat,ierr) do i=2,i1 do j=js,jf if(maxj(i-1,j-1).eq.1) then tauxs(i-1,j-1)=.2*(scr_data(i,j,1)+scr_data(i,j-1,1)+scr_data(i,j+1,1)+scr_data(i-1,j,1)+scr_data(i+1,j,1)) tauys(i-1,j-1)=.2*(scr_data(i,j,2)+scr_data(i,j-1,2)+scr_data(i,j+1,2)+scr_data(i-1,j,2)+scr_data(i+1,j,2)) endif enddo enddo enddo ! write(72) tauxs,tauys if(myid.eq.0)print*,'scow windstress month:',m else ! write(72) taux_data,tauy_data if(myid.eq.0)print*,'hr windstress month:',m endif cli_tauys(:,:,m) = tauys cli_tauxs(:,:,m) = tauxs clidata_tauy(:,:,m) = tauy_data clidata_taux(:,:,m) = taux_data enddo call wr8nc_winds(cli_tauxs,cli_tauys,clidata_taux,clidata_tauy) !dog call wr8nc_levitus(clidata_s1,clidata_t1) !dog end subroutine function r(t,s,p) p02=1747.4508988+t*(11.51588-0.046331033*t)-s*(3.85429655+0.01353985*t) p01=p/10d0+5884.81703666+t*(39.803732+t*(-0.3191477+t*0.0004291133))+2.6126277*s r=p01/(p02+0.7028423*p01) return end function function theta(sal,ts,ps,pr) ! ...sal=salinity ! ...ts=in situ temp (deg cel) ! ...ps=in situ pressure (dbars) ! ...pr=reference pressure (dbars) delp=pr-ps hafp=ps+.5*delp delt1=delp*gamma0(sal,ts,ps) t1_tmp=ts+.5*delt1 delt2=delp*gamma0(sal,t1_tmp,hafp) t2_tmp=t1_tmp+.2928932*(delt2-delt1) delt3=delp*gamma0(sal,t2_tmp,hafp) t3_tmp=t2_tmp+1.707107*(delt3-0.5857864*delt2-0.1213203*delt1) delt4=delp*gamma0(sal,t3_tmp,pr) t4_tmp=t3_tmp+0.16666667*(delt4-6.828427*delt3+4.828427*delt2+delt1) theta=t4_tmp return end function function gamma0(ss,tt,p) ! ...adiabati! temperature gradient(deg c/dbar) ! ...according to bryden (1973),dsr,20,401-408 ! ...copied from bio auxilary library xx=ss-35 gamma0=0.35803e-4+0.18932e-5*xx+p*(0.18741e-7- & 0.11351e-9*xx-0.46206e-12*p) & +tt*(0.85258e-5-0.42393e-7*xx+p*(-0.67795e-9+ & 0.27759e-11*xx+0.18676e-13*p) & +tt*(-0.68360e-7+p*(0.87330e-11-0.21687e-15*p) & +tt*(0.66228e-9-0.54481e-13*p))) return end function subroutine wr8nc_depth(rmets) use const, only: i0,j0,ngx0,ngy0 use comm, only: m_cart,mylat,mylon implicit none !declare pnetcdf var character :: fn*13 integer :: ierr,ncid,nd(2) integer :: i0id,j0id integer :: irmt integer(kind=mpi_offset_kind) :: dim_nn(2) integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2) real*8,dimension(i0,j0) :: rmets !output file name fn = 'data/depth.nc' !create ncfile and id ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid) !define dimension and assign dimID dim_nn(1) = i0*ngx0 dim_nn(2) = j0*ngy0 ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(1),i0id) ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(2),j0id) !define variables and assign varID nd(1) = i0id nd(2) = j0id ierr=nfmpi_def_var(ncid,"rmets",nf_double,2,nd,irmt) ierr=nfmpi_enddef(ncid) !put variables ns_ij(1) = mylon*i0+1 ns_ij(2) = mylat*j0+1 nc_ij(1) = i0 nc_ij(2) = j0 ierr=nfmpi_put_vara_double_all(ncid,irmt,ns_ij,nc_ij,rmets) ierr=nfmpi_close(ncid) end subroutine wr8nc_depth subroutine wr8nc_levitus(clidata_s1,clidata_t1) use const, only: i0,j0,k1,ngx0,ngy0 use comm, only: m_cart,mylat,mylon,myid implicit none !declare pnetcdf var character :: fn*20 integer :: ierr,ncid,nq(4) integer :: i0id,j0id,k1id,moid integer :: is1d,it1d,itmp integer(kind=mpi_offset_kind) :: dim_nn(4) integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4) real*8,dimension(i0,j0,k1,12),intent(in) :: clidata_s1,clidata_t1 !output file name fn = 'data/annualevitus.nc' !create ncfile and id ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid) !define dimension and assign dimID dim_nn(1) = i0*ngx0 dim_nn(2) = j0*ngy0 dim_nn(3) = k1 dim_nn(4) = 12 ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(1),i0id) ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(2),j0id) ierr=nfmpi_def_dim(ncid,"k1",dim_nn(3),k1id) ierr=nfmpi_def_dim(ncid,"month",dim_nn(4),moid) !define variables and assign varID nq(1) = i0id nq(2) = j0id nq(3) = k1id nq(4) = moid ierr=nfmpi_def_var(ncid,"s1_data",nf_double,4,nq,is1d) ierr=nfmpi_def_var(ncid,"t1_data",nf_double,4,nq,it1d) ierr=nfmpi_enddef(ncid) !put variables ns_ijkt(1) = mylon*i0+1 ns_ijkt(2) = mylat*j0+1 ns_ijkt(3) = 1 ns_ijkt(4) = 1 nc_ijkt(1) = i0 nc_ijkt(2) = j0 nc_ijkt(3) = k1 nc_ijkt(4) = 12 ierr=nfmpi_put_vara_double_all(ncid,is1d,ns_ijkt,nc_ijkt,clidata_s1) ierr=nfmpi_put_vara_double_all(ncid,it1d,ns_ijkt,nc_ijkt,clidata_t1) ierr=nfmpi_close(ncid) end subroutine wr8nc_levitus subroutine wr8nc_winds(cli_tauxs,cli_tauys,clidata_taux,clidata_tauy) use const, only: i2,j2,ngy0,ngx0 use comm, only: m_cart implicit none !declare pnetcdf var character :: fn*13 integer :: ierr,ncid,nt(3) integer :: i2id,j2id,moid integer :: ixs,iys,ixd,iyd integer(kind=mpi_offset_kind) :: dim_nn(3) integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3) real*8,dimension(i2,j2,12),intent(in) :: cli_tauxs,cli_tauys real*8,dimension(i2,j2,12),intent(in) :: clidata_taux,clidata_tauy !output file name fn = 'data/winds.nc' !create ncfile and id ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid) !define dimension and assign dimID dim_nn(1) = i2*ngx0 dim_nn(2) = j2*ngy0 dim_nn(3) = 12 ierr=nfmpi_def_dim(ncid,"i2_global",dim_nn(1),i2id) ierr=nfmpi_def_dim(ncid,"j2_global",dim_nn(2),j2id) ierr=nfmpi_def_dim(ncid,"month",dim_nn(3),moid) !define variables and assign varID nt(1) = i2id nt(2) = j2id nt(3) = moid ierr=nfmpi_def_var(ncid,"tauxs",nf_double,3,nt,ixs) ierr=nfmpi_def_var(ncid,"tauys",nf_double,3,nt,iys) ierr=nfmpi_def_var(ncid,"taux_data",nf_double,3,nt,ixd) ierr=nfmpi_def_var(ncid,"tauy_data",nf_double,3,nt,iyd) ierr=nfmpi_enddef(ncid) !put variables ns_ijk(1) = mylon*i2+1 ns_ijk(2) = mylat*j2+1 ns_ijk(3) = 1 nc_ijk(1) = i2 nc_ijk(2) = j2 nc_ijk(3) = 12 ierr=nfmpi_put_vara_double_all(ncid,ixs,ns_ijk,nc_ijk,cli_tauxs) ierr=nfmpi_put_vara_double_all(ncid,iys,ns_ijk,nc_ijk,cli_tauys) ierr=nfmpi_put_vara_double_all(ncid,ixd,ns_ijk,nc_ijk,clidata_taux) ierr=nfmpi_put_vara_double_all(ncid,iyd,ns_ijk,nc_ijk,clidata_tauy) ierr=nfmpi_close(ncid) end subroutine wr8nc_winds subroutine rdnc_westdeg use comm, only: m_cart,myid implicit none character :: fn*15 integer :: ncid,ierr,itmp integer(kind=mpi_offset_kind) :: ns_i,nc_i real :: test fn = 'data/westdeg.nc' ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid) ierr=nfmpi_inq_varid(ncid,'westdeg',itmp) ns_i = myid+1 nc_i = 1 ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,westdeg) ierr=nfmpi_close(ncid) end subroutine rdnc_westdeg subroutine rdnc_yzgrid(j,k,dxmnut) use comm, only: m_cart,myid,mylat,mylon use const, only: j0,j1 implicit none character :: fn*14 integer :: ncid,ierr,itmp integer(kind=mpi_offset_kind) :: ns_i,nc_i integer,intent(inout) :: j,k real*8,intent(inout) :: dxmnut fn = 'data/yzgrid.nc' ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid) ierr=nfmpi_inq_varid(ncid,'j0',itmp) ierr=nfmpi_get_var_int_all(ncid,itmp,j) ierr=nfmpi_inq_varid(ncid,'k0',itmp) ierr=nfmpi_get_var_int_all(ncid,itmp,k) ierr=nfmpi_inq_varid(ncid,'dxmnut',itmp) ierr=nfmpi_get_var_double_all(ncid,itmp,dxmnut) ierr=nfmpi_inq_varid(ncid,'tlz_depth',itmp) ierr=nfmpi_get_var_double_all(ncid,itmp,tlz) ierr=nfmpi_inq_varid(ncid,'z',itmp) ierr=nfmpi_get_var_double_all(ncid,itmp,z) ierr=nfmpi_inq_varid(ncid,'odz',itmp) ierr=nfmpi_get_var_double_all(ncid,itmp,odz) ierr=nfmpi_inq_varid(ncid,'odzw',itmp) ierr=nfmpi_get_var_double_all(ncid,itmp,odzw) ns_i=mylon*i0+1 nc_i=i0 ierr=nfmpi_inq_varid(ncid,'xdeg',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,xdeg) ns_i=mylat*j0+1 nc_i=j0 ierr=nfmpi_inq_varid(ncid,'y',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,y) ierr=nfmpi_inq_varid(ncid,'yv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,yv) ierr=nfmpi_inq_varid(ncid,'yvdeg',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,yvdeg) ierr=nfmpi_inq_varid(ncid,'ydeg',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ydeg) ierr=nfmpi_inq_varid(ncid,'cs',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,cs) ierr=nfmpi_inq_varid(ncid,'ocs',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ocs) ierr=nfmpi_inq_varid(ncid,'dx',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dx) ierr=nfmpi_inq_varid(ncid,'odx',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odx) ierr=nfmpi_inq_varid(ncid,'dy',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dy) ierr=nfmpi_inq_varid(ncid,'ody',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ody) ns_i=mylat*j1+1 nc_i=j1 ierr=nfmpi_inq_varid(ncid,'csv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,csv) ierr=nfmpi_inq_varid(ncid,'ocsv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,ocsv) ierr=nfmpi_inq_varid(ncid,'dxv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dxv) ierr=nfmpi_inq_varid(ncid,'odxv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odxv) ierr=nfmpi_inq_varid(ncid,'dyv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,dyv) ierr=nfmpi_inq_varid(ncid,'odyv',itmp) ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_i,nc_i,odyv) ierr=nfmpi_close(ncid) end subroutine rdnc_yzgrid end module