module presolver
use const
use comm
use var
implicit real*8 (a-h,o-z)
# include "pnetcdf.inc"
! preprocessor for 1/2 deg res global model july, 2007
!!cyclic longitudinal b.c.'s (i.e., global or periodic)
! ----------------------------------------------------------------------
! block data nishal
! ----------------------------------------------------------------------
real*8 hh,d1
common/contro/r0
common/lpi/mvfreq
common/boyanc/d1(i0,j0,k1)
common/bathy/hh(i0,j0),rmets(i0,j0)
common/bcs/tsurf(i2,j2,12),ssurf(i2,j2,12),snsp(i2,10,k2,12),sssp(i2,10,k2,12),tnsp(i2,10,k2,12),tssp(i2,10,k2,12)
! lsolver solver control
! lsolver:0 no solver setting
! lsolver:1 evp
! lsolver:2 bi-cgstab
! ie(nb0)=j1
! we initialize ie in do loop 250 of subroutine initfs
! data ie/6,11,16,...,j1/
! ----------------------------
! explanation of parameters
! ----------------------------
! i0,j0: horizontal dimension parameters (includes ghost zones around
! entire model domain; climatological boundary values are stored in
! ghost zones)
! k0: vertical dimension parameter (number of layer interfaces, equal
! the number of layers or pressure levels plus 1, does not include ghost
! zone)
! nb0: number of evp ellipti! solver restart vectors, used to control
! solver roundoff error (decreases with increasing nb0)
! case=5-character case name
! dscrib=63-character case descriptor
! r0=earth radius
! daodt=number of time steps per day
! b=thermal expansion coefficient when linear equation of state is used
! g=gravity
! dm0,de0 are horizontal heat and momentum diffusivities (cm-cm/sec)
! rzmx=vertical cell reynolds number limit
! drag=drag coefficient (standard value is .002)
! fltw=time filter coefficient on time filtered leapfrog method
! ktrm=thermocline level (used only for animation graphics)
! lrstrt=flag to initialize from restart file; if lrstrt=0, initial vel
! is zero (but not necessarily boundary values)
! mxit=last time step number of current run (last time step of previous
! run is saved on restart file)
! lopen=flag for open lateral boundaries (derive geostrophi! inflows)
! lmovi=flag to save data for animation of the results
! lsolver=flag for this program to initialize ellipti! solver data for the
! lprecon=flag for precodition of bi-cgstab(lsolver=2)
! given geometry and also to generate initial temperature and salinity
! ie=array for evp ellipti! solver (controls roundoff error); ie(nb0)=j1
contains
! ----------------------------------------------------------------------
subroutine solvergen
! ----------------------------------------------------------------------
! input data: etopo5, levitus (t and s)
! the model now reads input wind data directly
! output data:
! input data translated to model variable latitudinal resolution grid
integer*2 intmp(i2,j2),ing(i0t,j0t),inbuf(i2t*j2t)
integer*2 isw(j0*k1),ise(j0*k1),isn(i0*k1),iss(i0*k1),irw(j0*k1),ire(j0*k1),irn(i0*k1),irs(i0*k1)
real*8 scr_tmp(i0,j0)
daodt=576
r0=6.4e8
b=2.05e-4
g=980.
dm0=5.e6
de0=25.e5
dmz0=.01
rzmx=50
drag=.001
fltw=.1
rawf=.53
! this may be unstable
!data ktrm,lrstrt,mxit,isav,mxsav/5,0,0,5400,5400/
ktrm=5
lrstrt=0
mxit=0
isav=5400
mxsav=5400
!data lopen,lmovi,mvfreq,lfsrf,lsolver,lprecon/0,1,5,1,2,1/
lopen=0
lmovi=1
mvfreq=5
lfsrf=1
lsolver=2
lprecon=1
case="GLOBO"
dscrib="global model (.5 deg longitudinal resolution) "
if (myid .eq. 0) then
open(79,file='namelist.txt',form='formatted')
read(79,*)lrstrt
read(79,*)daodt
read(79,*)irun
close(79)
mxit=irun*daodt*360
endif
call mpi_bcast(lrstrt,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(daodt,1,mpi_real8,0,mpi_comm_world,ierr)
call mpi_bcast(mxit,1,mpi_integer,0,mpi_comm_world,ierr)
! output data
! open(10,file='data/yzgrid'//grd,form='unformatted')
! open(11,file='data/depth'//grd,form='unformatted')
! open(70,file='data/annualevitus'//grd,form='unformatted')
! output data
! open(7,file='data/vbk_ascii'//grd)
! open(8,file='data/windmix'//grd,form='unformatted')
! open(9,file='data/rundata'//grd,form='unformatted')
! open(17,file='data/kbview'//grd)
! open(60,file='data/zkb'//grd,form='unformatted')
! open(62,file='data/boundaries'//grd,form='unformatted')
! check solver
if(lsolver.eq.1 .or. lsolver.eq.2) then
! open(99,file='data/evp'//grd,form='unformatted')
! rewind 99
elseif(lsolver.eq.0) then
if(myid.eq.0) print*,'no solver setting'
else
if(myid.eq.0) print*,'lsolver=',lsolver,' err'
stop 0
endif
! rewind 10
! rewind 11
! rewind 70
! rewind 7
! rewind 8
! rewind 9
! rewind 17
! rewind 60
! rewind 62
call rdnc_yzgrid2(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 i0,k0 to values from metgen.f.')
stop 0
endif
! continue
!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
!k01=k0+k1
!write(21,10) (.01*z(k),k=1,k01,2)
!10 format('z-levels(m)'/(10f8.2))
!coriolis parameter array
omega2=3.141592654/21600.
! degrees to radians
pi_180=3.141592654/180.
do 12 j=2,j1
! spherical curvature parameter
tanphi(j-1)=tan(pi_180*ydeg(j))/r0
12 f(j-1)=omega2*sin(pi_180*ydeg(j))
! write(21,14) i0,j0
!14 format('longitude(x) grid dimension:',i4/'latitude(y) grid dimension:',i4/)
! rmets data
! rmets(i,1) corresponds to southern most internal grid point
! stop
call rdnc_depth(rmets)
! read(11) rmets
write(*,*) '==============',maxval(rmets),minval(rmets)
write(*,*) 'rmets filtered'
write(*,*) '=============='
call mpi_type_vector(j0,1,i0,mpi_real8,m_vlon_2d,ierr) !dog
call mpi_type_vector(1,i0,j0,mpi_real8,m_vlat_2d,ierr)
call mpi_type_commit(m_vlon_2d,ierr)
call mpi_type_commit(m_vlat_2d,ierr)
call mpi_sendrecv(rmets(2,1),1,m_vlon_2d,m_w,1,rmets(i0,1),1,m_vlon_2d,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(i1,1),1,m_vlon_2d,m_e,1,rmets(1,1),1,m_vlon_2d,m_w,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(1,2),1,m_vlat_2d,m_s,1,rmets(1,j0),1,m_vlat_2d,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(1,j1),1,m_vlat_2d,m_n,1,rmets(1,1),1,m_vlat_2d,m_s,1,mpi_comm_world,istat,ierr)
!cts change from mxpas=2 to 4
mxpas=4
do n=1,mxpas
do 15 j=2,j1
do 15 i=2,i1
! don't filter izu ridge and se honshu coastal region
if(xdeg(i).gt.134.5.and.xdeg(i).lt.147..and.ydeg(j).gt.-62.5.and.ydeg(j).lt.-51.5) goto 15
! don't filter gs window region
if(xdeg(i).gt.284.5.and.xdeg(i).lt.320..and.ydeg(j).gt.-54.5.and.ydeg(j).lt.-29.) goto 15
! don't filter regions near land (they are adequately known)
! filter only points farther than nwet grid intervals from land
if (rmets(i,j).lt. 300.) goto 15
scr_tmp(i,j)=.25*(rmets(i,j-1)+rmets(i,j+1)+rmets(i-1,j)+rmets(i+1,j))
15 continue
do 16 j=2,j1
do 16 i=2,i1
! don't filter izu ridge and se honshu coastal region
if(xdeg(i).gt.134.5.and.xdeg(i).lt.147..and.ydeg(j).gt.-62.5.and.ydeg(j).lt.-51.5) goto 16
! don't filter gs window region
if(xdeg(i).gt.284.5.and.xdeg(i).lt.320..and.ydeg(j).gt.-54.5.and.ydeg(j).lt.-29.) goto 16
if (rmets(i,j).lt. 300.) goto 16
rmets(i,j)=scr_tmp(i,j)
16 continue
call mpi_sendrecv(rmets(2,1),1,m_vlon_2d,m_w,1,rmets(i0,1),1,m_vlon_2d,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(i1,1),1,m_vlon_2d,m_e,1,rmets(1,1),1,m_vlon_2d,m_w,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(1,2),1,m_vlat_2d,m_s,1,rmets(1,j0),1,m_vlat_2d,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(1,j1),1,m_vlat_2d,m_n,1,rmets(1,1),1,m_vlat_2d,m_s,1,mpi_comm_world,istat,ierr)
enddo
!cts20150921
if(.false.)then
! 555 510 64 53
! make izu ridge into a narrow ridge at sill level
!cw do 18 j=547,527,-1
!cw jloc=j-mylat*j2
do 18 j=2,j1
do 18 i=3,i2
! rmets is in meters
! 15 write(*,16) j,(int(.1*rmets(i,j)),i=555,565)
! 16 format(i3,1x,11i4)
! find shallowest sill level in longitudinal range i=555-565
kmn=10000
!cw do 17 i=556,564
!cw iloc=i-mylon*i2
! if(iloc.le.2.or.iloc.ge.i1.or.jloc.le.2.and.jloc.ge.j1)
!cw if(iloc.lt.2.or.iloc.gt.i1.or.jloc.lt.2.or.jloc.gt.j1)
!cw 1 goto 18
if(xdeg(i).gt.138..and.xdeg(i).lt.142..and.ydeg(j).gt.30..and.ydeg(j).lt.35.) then
! k=kmn
! kmn=min(kmn,int(.1*rmets(i,j)))
! 17 if (kmn.ne.k) imn=i
! make narrow ridge
! tmp=max(rmets(imn-1,j),rmets(imn-2,j))
! rmets(imn-2,j)=tmp
! rmets(imn-1,j)=tmp
! tmp=max(rmets(imn+1,j),rmets(imn+2,j))
! rmets(imn+2,j)=tmp
! rmets(imn+1,j)=tmp
! write(*,*) imn,kmn
k=kmn
!cw kmn=min(kmn,int(.1*rmets(iloc,jloc)))
kmn=min(kmn,int(.1*rmets(i,j)))
!cw 17 if (kmn.ne.k) imn=iloc
17 if (kmn.ne.k) imn=i
! make narrow ridge
tmp=max(rmets(imn-1,j),rmets(imn-2,j))
rmets(imn-2,j)=tmp
rmets(imn-1,j)=tmp
tmp=max(rmets(imn+1,j),rmets(imn+2,j))
rmets(imn+2,j)=tmp
rmets(imn+1,j)=tmp
!cw tmp=max(rmets(imn-1,jloc),rmets(imn-2,jloc))
!cw rmets(imn-2,jloc)=tmp
!cw rmets(imn-1,jloc)=tmp
!cw tmp=max(rmets(imn+1,jloc),rmets(imn+2,jloc))
!cw rmets(imn+2,jloc)=tmp
!cw rmets(imn+1,jloc)=tmp
! write(*,*) imn,kmn
endif
18 continue
do i=2,i1
do j=2,j1
! smoother coastline near cozumel
if(xdeg(i).gt.273.2.and.xdeg(i).lt.273.5.and.ydeg(j).gt.21.1.and.ydeg(j).lt.21.5) then
rmets(i,j)=0.
endif
! eliminate shallow water northwest of easternmost florida point
! rmets(1118:1120,516:530)=0.
if(xdeg(i).gt.279.1.and.xdeg(i).lt.279.9.and.ydeg(j).gt.28.1.and.ydeg(j).lt.31.6) then
rmets(i,j)=0.
endif
! rmets(1121,516:530)=28.
if(xdeg(i).gt.280.0.and.xdeg(i).lt.280.3.and.ydeg(j).gt.28.1.and.ydeg(j).lt.31.6) then
rmets(i,j)=28.
endif
! rmets(1121,521:525)=45.
if(xdeg(i).gt.280.0.and.xdeg(i).lt.280.3.and.ydeg(j).gt.29.2.and.ydeg(j).lt.30.5) then
rmets(i,j)=45.
endif
! rmets(1120,522:523)=28.
if(xdeg(i).gt.279.5.and.xdeg(i).lt.280.0.and.ydeg(j).gt.29.5.and.ydeg(j).lt.30.0) then
rmets(i,j)=28.
endif
if(ydeg(j).gt.85.0)then
rmets(i,j)=0.
endif
enddo
enddo
! stop
! artificial high latitude shelves (shortcircuits arcti! ocean,
! truncates antarcti! seas and bays, and avoids artificial excessive bt
! mode vortex stretching at latitudinal boundaries)
!cw new shallowing scheme
jtt=min(j2t/60,k2)
!cw jtt=min(j2t/30,k2)
kinv=int(dble(k2)/dble(jtt)+.5)
kzl=int(70.*dble(k1)/60.)
if(myid.eq.0) print*,'kzl=',kzl
do 23 j=2,j1
if(mylat*j2+j-1.le.jtt) then
jin=mylat*j2+j-2
kin=1+kinv*jin
!cw kk=min(k0+k1,2*kin+3)
kk=min(k0+k1,2*kin+kzl+1)
tmp=.01*z(kk)
do 21 i=1,i0
21 rmets(i,j)=min(rmets(i,j),tmp)
endif
if(mylat*j2+j.ge.j0t-jtt) then
jin=j0t-mylat*j2-j-1
kin=1+kinv*jin
!cw kk=min(k0+k1,2*kin+3)
kk=min(k0+k1,2*kin+kzl+1)
tmp=.01*z(kk)
do 22 i=1,i0
22 rmets(i,j)=min(rmets(i,j),tmp)
endif
23 continue
!cts20150921 if(.false.)then
endif
! if (mylat .eq. 0) then
! n=0
! do 22 k=1,k1
! n=n+1
! kk=min(k0+k1,2*k+3)
! tmp=.01*z(kk)
! do 22 i=1,i0
! 22 rmets(i,n+1)=min(rmets(i,n+1),tmp)
! endif
!
! if (mylat .eq. ngy1) then
! n=0
! do 23 k=1,k1
! n=n+1
! kk=min(k0+k1,2*k+3)
! tmp=.01*z(kk)
! do 23 i=1,i0
! 23 rmets(i,j0-n)=min(rmets(i,j0-n),tmp)
! endif
! bad point is (39,401) starting at step 91100
! write(*,24) (j,(int(rmets(i,j)),i=31,47),j=j1,j1-50,-1)
! 24 format(i3,4x,17i4)
! write(*,24) (j,(int(rmets(i,j)),i=261,300),j=206,167,-1)
!c24 format(/(i3,1x,40i1))
! make all water least 2 layers deep
! to allow baroclini! ventillation all the way to coast
tmp=.01*z(5)
!convert 1-layer deep water to 2-layer-deep water
! alternatively, !convert water less than 2 layers deep to land
!cw20150910
!cw tmp1=.01*(z(3)+z(5))/2.
tmp1=.01*(z(5)+z(7))/2.
do 25 j=2,j1
do 25 i=2,i1
!convert land to water
! 25 if (rmets(i,j).gt.tmp1) rmets(i,j)=max(rmets(i,j),tmp)
!convert water to land
25 if (rmets(i,j).lt.tmp1) rmets(i,j)=0.
! periodic
call mpi_sendrecv(rmets(2,1),1,m_vlon_2d,m_w,1,rmets(i0,1),1,m_vlon_2d,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(i1,1),1,m_vlon_2d,m_e,1,rmets(1,1),1,m_vlon_2d,m_w,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(1,2),1,m_vlat_2d,m_s,1,rmets(1,j0),1,m_vlat_2d,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(rmets(1,j1),1,m_vlat_2d,m_n,1,rmets(1,1),1,m_vlat_2d,m_s,1,mpi_comm_world,istat,ierr)
do 33 j=1,j0
do 33 i=1,i0
33 rmets(i,j)=min(tlz,100.*rmets(i,j))
! determine "logical grid" rmets array kb(i,j)
! kb(i,j)=number of layers at pressure point (i,j)
do 34 k=1,k1
zz=.5*(z(2*k-1)+z(2*k+1))
do 34 j=1,j0
do 34 i=1,i0
! if rmets is deeper than mid-layer (zz), increment kb by 1
l=rmets(i,j)/zz
34 kb(i,j)=kb(i,j)+min(1,l)
! widen northern izu gap to allow more flow through to entrain more
! oyashio water. this was added at day 90 of year 6. before that the
! kuroshio extension was good, but a small southward loop downstream from
! the gap is decreased by this change.
! deepen widened gap at start of year 14
! n=1
! further deepen gap at start of year 16
n=2
do i=2,i1
do j=2,j1
if(.false.)then
!cts 20150911 remove new york bight singularity
if(xdeg(i).gt.289.5.and.xdeg(i).lt.293.5.and.ydeg(j).gt.42.2.and.ydeg(j).lt.45.2) then
kb(i,j)=min(kb(i,j),3)
endif
if(xdeg(i).gt.289.5.and.xdeg(i).lt.293.5.and.ydeg(j).gt.41.5.and.ydeg(j).lt.45.2) then
kb(i,j)=min(kb(i,j),7)
endif
endif
!cts seamounts around 44.7n, 321e
if(xdeg(i).gt.319.and.xdeg(i).lt.322.and.ydeg(j).gt.44.8.and.ydeg(j).lt.45.2) then
kb(i,j)=min(kb(i,j),32)
endif
if(.false.)then
!cts kuroshio seperation-japan coast should be sahllower
if(xdeg(i).gt.138.8.and.xdeg(i).lt.141.4.and.ydeg(j).gt.34.5.and.ydeg(j).lt.35.) then
kb(i,j)=min(kb(i,j),18)
endif
if(xdeg(i).gt.140.9.and.xdeg(i).lt.141.8.and.ydeg(j).gt.35..and.ydeg(j).lt.36.) then
kb(i,j)=min(kb(i,j),15)
endif
if(xdeg(i).gt.140.9.and.xdeg(i).lt.141.8.and.ydeg(j).gt.36.and.ydeg(j).lt.37.) then
kb(i,j)=min(kb(i,j),11)
endif
if(xdeg(i).gt.141.8.and.xdeg(i).lt.142.8.and.ydeg(j).gt.35.1.and.ydeg(j).lt.39.4) then
kb(i,j)=kb(i,j)-5
endif
endif
!cw 20150910
if(.false.) then
!cw do j=539,541 i=561,563
!cw kb(iloc,jloc)=13+n
if(xdeg(i).gt.140.0.and.xdeg(i).lt.140.8.and.ydeg(j).gt.33.2.and.ydeg(j).lt.34.0) then
kb(i,j)=13+n
endif
!cw iloc=561-mylon*i2 jloc=539-mylat*j2
!cw kb(iloc,jloc)=13+n
if(xdeg(i).gt.140.0.and.xdeg(i).lt.140.3.and.ydeg(j).gt.33.2.and.ydeg(j).lt.33.6) then
kb(i,j)=13+n
endif
!cw iloc=562-mylon*i2 jloc=539-mylat*j2
!cw kb(iloc,jloc)=13+n
if(xdeg(i).gt.140.3.and.xdeg(i).lt.140.5.and.ydeg(j).gt.33.2.and.ydeg(j).lt.33.6) then
kb(i,j)=13+n
endif
!cw iloc=561-mylon*i2 jloc=540-mylat*j2
!cw kb(iloc,jloc)=13+n
if(xdeg(i).gt.140.0.and.xdeg(i).lt.140.3.and.ydeg(j).gt.33.5.and.ydeg(j).lt.33.7) then
kb(i,j)=13+n
endif
!cw iloc=562-mylon*i2 jloc=540-mylat*j2
!cw kb(iloc,jloc)=13+n
if(xdeg(i).gt.140.3.and.xdeg(i).lt.140.5.and.ydeg(j).gt.33.5.and.ydeg(j).lt.33.7) then
kb(i,j)=13+n
endif
!cw iloc=562-mylon*i2 jloc=541-mylat*j2
!cw kb(iloc,jloc)=13+n
if(xdeg(i).gt.140.3.and.xdeg(i).lt.140.5.and.ydeg(j).gt.33.7.and.ydeg(j).lt.33.9) then
kb(i,j)=13+n
endif
! eliminate anomalous ny bight island
!cw iloc=1142-mylon*i2 jloc=557-mylat*j2
!cw kb(iloc,jloc)=2
if(xdeg(i).gt.285.3.and.xdeg(i).lt.285.5.and.ydeg(j).gt.36.9.and.ydeg(j).lt.37.1) then
kb(i,j)=2
endif
!cw iloc=1142-mylon*i2 jloc=558-mylat*j2
!cw kb(iloc,jloc)=2
if(xdeg(i).gt.285.3.and.xdeg(i).lt.285.5.and.ydeg(j).gt.37.1.and.ydeg(j).lt.37.3) then
kb(i,j)=2
endif
!cw iloc=1142-mylon*i2 jloc=559-mylat*j2
!cw kb(iloc,jloc)=2
if(xdeg(i).gt.285.3.and.xdeg(i).lt.285.5.and.ydeg(j).gt.37.3.and.ydeg(j).lt.37.5) then
kb(i,j)=2
endif
!cw iloc=1144-mylon*i2 jloc=561-mylat*j2
!cw kb(iloc,jloc)=2
if(xdeg(i).gt.285.8.and.xdeg(i).lt.290.1.and.ydeg(j).gt.37.7.and.ydeg(j).lt.37.9) then
kb(i,j)=2
endif
!cts no above modification
endif
!cts new check and topographical modification
if(xdeg(i).gt.140.5.and.xdeg(i).lt.141.5.and.ydeg(j).gt.41.0.and.ydeg(j).lt.42.) then
kb(i,j)=10
endif
!cw20150902
if(xdeg(i).gt.120.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.24.5.and.ydeg(j).lt.25.0) then
kb(i,j)=9
endif
if(xdeg(i).gt.121.5.and.xdeg(i).lt.122.5.and.ydeg(j).gt.25.5.and.ydeg(j).lt.26.0) then
kb(i,j)=13
endif
!cw indonesia throughflow
if(xdeg(i).gt.119.5.and.xdeg(i).lt.120.5.and.ydeg(j).gt.-6..and.ydeg(j).lt.0.25) then
kb(i,j)=23
endif
if(xdeg(i).gt.119.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.0.5.and.ydeg(j).lt.0.75) then
kb(i,j)=20
endif
if(xdeg(i).gt.119.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.0.75.and.ydeg(j).lt.1.) then
kb(i,j)=24
endif
if(xdeg(i).gt.119.5.and.xdeg(i).lt.121.5.and.ydeg(j).gt.1..and.ydeg(j).lt.1.5) then
kb(i,j)=26
endif
!cw allow lombok
if(xdeg(i).gt.118.5.and.xdeg(i).lt.119.5.and.ydeg(j).gt.-12.5.and.ydeg(j).lt.-11.75) then
kb(i,j)=26
endif
if(xdeg(i).gt.118.5.and.xdeg(i).lt.119.5.and.ydeg(j).gt.-10.75.and.ydeg(j).lt.-10.25) then
kb(i,j)=28
endif
!cw allow ombai
if(xdeg(i).gt.122.5.and.xdeg(i).lt.123.5.and.ydeg(j).gt.-8.5.and.ydeg(j).lt.-8.) then
kb(i,j)=26
endif
if(xdeg(i).gt.123.5.and.xdeg(i).lt.124.5.and.ydeg(j).gt.-8.5.and.ydeg(j).lt.-8.) then
kb(i,j)=25
endif
if(xdeg(i).gt.126.5.and.xdeg(i).lt.127.5.and.ydeg(j).gt.-9.and.ydeg(j).lt.-8.5) then
kb(i,j)=24
endif
if(xdeg(i).gt.129.5.and.xdeg(i).lt.131.5.and.ydeg(j).gt.-3.75.and.ydeg(j).lt.0.75) then
kb(i,j)=24
endif
enddo
enddo
do j=2,j1
do i=2,i1
if(ydeg(j).ge.34..and.ydeg(j).le.36.5.and.xdeg(i).ge.352.5.and.xdeg(i).le.356)then
kb(i,j)=max(kb(i,j),19)
endif
enddo
enddo
do j=2,j1
do i=2,i1
if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
write(*,*)'topo1',i,j
do jj=j+6,j-8,-1
write(*,'(a,11i4)')'topo1',kb(i-5:i+5,jj)
enddo
endif
enddo
enddo
!cw check
! do i=2,i1
! do j=2,j1
! it=i+mylon*i2-1
! jt=j+mylat*j2-1
! if(it.eq.107.and.jt.eq.144) print*,'1 kb 107 144=',kb(i,j)
! 1 ,xdeg(i),ydeg(j)
! if(it.eq.116.and.jt.eq.144) print*,'1 kb 116 144=',kb(i,j)
! 1 ,xdeg(i),ydeg(j)
! if(it.eq.117.and.jt.eq.144) print*,'1 kb 117 144=',kb(i,j)
! enddo
! enddo
! derive land/sea mask array
do 52 j=1,j0
do 52 i=1,i0
it=kb(i,j)
if (it.ne.0) go to 50
kb(i,j)=1
go to 52
50 do 51 k=1,it
51 in(i,j,k)=1
52 continue
do 53 j=2,j1
do 53 i=2,i1
53 intmp(i-1,j-1)=in(i,j,1)
do j=2,j1
do i=2,i1
if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
write(*,*)'topo1',i,j
do jj=j+6,j-8,-1
write(*,'(a,11i4)')'topo1',in(i-5:i+5,jj,1)
enddo
endif
enddo
enddo
nsend=i2*j2
call mpi_gather(intmp,nsend,mpi_integer2,inbuf,nsend,mpi_integer2,0,mpi_comm_world,ierr)
if (myid .eq. 0) then
do 57 nx=1,ngx0
ngxt=(nx-1)*ngy0
do 57 ny=1,ngy0
ngt=((ny-1)+ngxt)*i2*j2
do 57 j=1,j2
nt=ngt+(j-1)*i2
jt=j+(ny-1)*j2
do 57 i=1,i2
! it=mod(i+(nx-1)*i2+i2t/2,i2t)
it=i+(nx-1)*i2
57 ing(it+1,jt+1)=inbuf(nt+i)
! in and ing is the same with in in global_serial1
!remark that boundary condition is computed in value of ghost cell.
! ------------------------------------------------------
! start with single point in central equatorial pacific
! check pt
n=0
! ing(i2t/6,j0t/2)=2
ing(2*i0t/3,j0t/2)=2
! ------------------------------------------------------
! due to periodi! logical grid boundary at 0 deg longitude
! we need a second point in the arcti! ocean
! rmk: i=2 alway lies in processor 0.
!cw20140604 ing(2,j1t)=2
153 n=n+1
nsum=0
do 155 j=2,j1t
do 154 i=2,i1t
nn=ing(i+1,j)
if (ing(i+1,j).eq.1.and.ing(i,j).eq.2) ing(i+1,j)=2
nsum=nsum+ing(i+1,j)-nn
nn=ing(i-1,j)
if (ing(i-1,j).eq.1.and.ing(i,j).eq.2) ing(i-1,j)=2
nsum=nsum+ing(i-1,j)-nn
nn=ing(i,j+1)
if (ing(i,j+1).eq.1.and.ing(i,j).eq.2) ing(i,j+1)=2
nsum=nsum+ing(i,j+1)-nn
nn=ing(i,j-1)
if (ing(i,j-1).eq.1.and.ing(i,j).eq.2) ing(i,j-1)=2
154 nsum=nsum+ing(i,j-1)-nn
do 155 i=i1t,2,-1
nn=ing(i+1,j)
if (ing(i+1,j).eq.1.and.ing(i,j).eq.2) ing(i+1,j)=2
nsum=nsum+ing(i+1,j)-nn
nn=ing(i-1,j)
if (ing(i-1,j).eq.1.and.ing(i,j).eq.2) ing(i-1,j)=2
nsum=nsum+ing(i-1,j)-nn
nn=ing(i,j+1)
if (ing(i,j+1).eq.1.and.ing(i,j).eq.2) ing(i,j+1)=2
nsum=nsum+ing(i,j+1)-nn
nn=ing(i,j-1)
if (ing(i,j-1).eq.1.and.ing(i,j).eq.2) ing(i,j-1)=2
155 nsum=nsum+ing(i,j-1)-nn
write(*,156) nsum,n
156 format(i8,'changes during iteration',i4)
if (nsum.ne.0) go to 153
do 157 j=1,j0t
ing(i0t,j)=ing(2,j)
157 ing(1,j)=ing(i1t,j)
endif
call mpi_bcast(ing,i0t*j0t,mpi_integer2,0,mpi_comm_world,ierr)
! widen northern izu gap to allow more flow through to entrain more
! oyashio water. this was added at day 90 of year 6. before that the
! kuroshio extension was good, but a small southward loop downstream from
! the gap is decreased by this change.
!cw iloc=561-mylon*i2
!cw jloc=539-mylat*j2
!cw if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw 1 kb(iloc,jloc)=13
!cw iloc=562-mylon*i2
!cw jloc=539-mylat*j2
!cw if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw 1 kb(iloc,jloc)=13
!cw iloc=561-mylon*i2
!cw jloc=540-mylat*j2
!cw if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw 1 kb(iloc,jloc)=13
!cw iloc=562-mylon*i2
!cw jloc=540-mylat*j2
!cw if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw 1 kb(iloc,jloc)=13
!cw iloc=562-mylon*i2
!cw jloc=541-mylat*j2
!cw if(iloc.ge.2.and.iloc.le.i1.and.jloc.ge.2.and.jloc.le.j1)
!cw 1 kb(iloc,jloc)=13
it=mylon*i2
jt=mylat*j2
do 158 j=1,j0
do 158 i=1,i0
in(i,j,1)=ing(it+i,jt+j)
n=in(i,j,1)
158 if (in(i,j,1).ne.2) kb(i,j)=0
do j=2,j1
do i=2,i1
if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
write(*,*)'topo2',i,j
do jj=j+6,j-8,-1
write(*,'(a,11i4)')'toin2',in(i-5:i+5,jj,1)
enddo
endif
enddo
enddo
do 159 k=1,k1
do 159 j=1,j0
do 159 i=1,i0
159 in(i,j,k)=0
do j=2,j1
do i=2,i1
if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
write(*,*)'topo2',i,j
do jj=j+6,j-8,-1
write(*,'(a,11i4)')'topo2',kb(i-5:i+5,jj)
enddo
endif
enddo
enddo
! do 8886 j=j1,2,-1
!8886 write(*,8887) j,(kb(i,j),i=1,122)
!8887 format(i3,1x,122i1)
! do 8888 j=j1,2,-1
!8888 write(*,8887) j,(kb(i,j),i=121,242)
! do 8889 j=j1,2,-1
!8889 write(*,8887) j,(kb(i,j),i=241,362)
! redefine in with exterraneous pools eliminated
do 162 j=2,j1
do 162 i=1,i0
it=kb(i,j)
if (it.ne.0) go to 160
kb(i,j)=1
go to 162
160 do 161 k=1,it
161 in(i,j,k)=1
162 continue
do j=2,j1
do i=2,i1
if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
write(*,*)'topo3',i,j
do jj=j+6,j-8,-1
write(*,'(a,11i4)')'topo3',kb(i-5:i+5,jj)
enddo
endif
enddo
enddo
! periodic
call mpi_type_vector(1,i0,j0,mpi_integer2,m_vlat_2di,ierr)
call mpi_type_commit(m_vlat_2di,ierr)
call mpi_sendrecv(kb(1,2),1,m_vlat_2di,m_s,1,kb(1,j0),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(kb(1,j1),1,m_vlat_2di,m_n,1,kb(1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)
call mpi_type_vector(k1,i0,i0*j0,mpi_integer2, m_vlat_2di,ierr)
call mpi_type_commit(m_vlat_2di,ierr)
call mpi_sendrecv(in(1,2,1),1,m_vlat_2di,m_s,1,in(1,j0,1),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(in(1,j1,1),1,m_vlat_2di,m_n,1,in(1,1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)
!cw check
! do i=2,i1
! do j=2,j1
! it=i+(mylon)*i2-1
! jt=j+(mylat)*j2-1
! if(it.eq.107.and.jt.eq.144) print*,'2 kb 107 144=',kb(i,j)
! if(it.eq.116.and.jt.eq.144) print*,'2 kb 116 144=',kb(i,j)
! if(it.eq.117.and.jt.eq.144) print*,'2 kb 117 144=',kb(i,j)
! enddo
! enddo
! mask out points having no lateral ventillation
69 n=0
do 70 k=2,k1
do 70 j=2,j1
do 70 i=2,i1
if (in(i-1,j,k)+in(i+1,j,k)+in(i,j-1,k)+in(i,j+1,k).ne.0) go to 70
if (in(i,j,k).ne.0) n=n+1
in(i,j,k)=0
70 continue
write(*,71) n
71 format('mask out',i4,' points having no lateral ventillation')
if (n.ne.0) go to 69
! redefine kb based on modified in
do 72 k=2,k1
do 72 j=2,j1
do 72 i=2,i1
if (in(i,j,k).eq.1) go to 72
if (in(i,j,k-1).eq.1) kb(i,j)=k-1
! this may be redundant
kb(i,j)=max(1,kb(i,j))
72 continue
do j=2,j1
do i=2,i1
if(ydeg(j).ge.36..and.ydeg(j).le.36.5.and.xdeg(i).ge.354.and.xdeg(i).le.356)then
write(*,*)'topo4',i,j
do jj=j+6,j-8,-1
write(*,'(a,11i4)')'topo4',kb(i-5:i+5,jj)
enddo
endif
enddo
enddo
! periodic
call mpi_barrier(mpi_comm_world,ierr)
call mpi_type_vector(1,i0,j0,mpi_integer2,m_vlat_2di,ierr)
call mpi_type_commit(m_vlat_2di,ierr)
call mpi_type_vector(j0,1,i0,mpi_integer2,m_vlon_2di,ierr)
call mpi_type_commit(m_vlon_2di,ierr)
call mpi_sendrecv(kb(1,2),1,m_vlat_2di,m_s,1,kb(1,j0),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(kb(1,j1),1,m_vlat_2di,m_n,1,kb(1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(kb(2,1),1,m_vlon_2di,m_w,1,kb(i0,1),1,m_vlon_2di,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(kb(i1,1),1,m_vlon_2di,m_e,1,kb(1,1),1,m_vlon_2di,m_w,1,mpi_comm_world,istat,ierr)
call mpi_barrier(mpi_comm_world,ierr)
call mpi_type_vector(k1,i0,i0*j0,mpi_integer2,m_vlat_2di,ierr)
call mpi_type_commit(m_vlat_2di,ierr)
call mpi_type_vector(k1*j0,1,i0,mpi_integer2,m_vlon_2di,ierr)
call mpi_type_commit(m_vlon_2di,ierr)
call mpi_sendrecv(in(1,2,1),1,m_vlat_2di,m_s,1,in(1,j0,1),1,m_vlat_2di,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(in(1,j1,1),1,m_vlat_2di,m_n,1,in(1,1,1),1,m_vlat_2di,m_s,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(in(2,1,1),1,m_vlon_2di,m_w,1,in(i0,1,1),1,m_vlon_2di,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(in(i1,1,1),1,m_vlon_2di,m_e,1,in(1,1,1),1,m_vlon_2di,m_w,1,mpi_comm_world,istat,ierr)
call mpi_barrier(mpi_comm_world,ierr)
!cw check
! do i=2,i1
! do j=2,j1
! it=i+(mylon)*i2-1
! jt=j+(mylat)*j2-1
! if(it.eq.107.and.jt.eq.144) print*,'3 kb 107 144=',kb(i,j)
! if(it.eq.116.and.jt.eq.144) print*,'3 kb 116 144=',kb(i,j)
! if(it.eq.117.and.jt.eq.144) print*,'3 kb 117 144=',kb(i,j)
! enddo
! enddo
! other mask arrays
do 64 k=1,k1
do 64 j=1,j0
do 64 i=1,i1
64 iu(i,j,k)=in(i,j,k)*in(i+1,j,k)
do 65 k=1,k1
do 65 j=1,j1
do 65 i=1,i0
65 iv(i,j,k)=in(i,j,k)*in(i,j+1,k)
do 66 k=1,k2
do 66 j=1,j0
do 66 i=1,i0
66 iw(i,j,k+1)=in(i,j,k)*in(i,j,k+1)
do k=1,k1
jt=(k-1)*j0
do j=1,j0
ise(jt+j)=iu(i1,j,k)
isw(jt+j)=iu(2,j,k)
enddo
enddo
do k=1,k1
it=(k-1)*i0
do i=1,i0
isn(it+i)=iv(i,j2,k)
iss(it+i)=iv(i,2,k)
enddo
enddo
nlon=j0*k1
nlat=i0*k1
call mpi_barrier(mpi_comm_world,ierr)
call mpi_sendrecv(isw,nlon,mpi_integer2,m_w,1,ire,nlon,mpi_integer2,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(ise,nlon,mpi_integer2,m_e,1,irw,nlon,mpi_integer2,m_w,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(iss,nlat,mpi_integer2,m_s,1,irn,nlat,mpi_integer2,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(isn,nlat,mpi_integer2,m_n,1,irs,nlat,mpi_integer2,m_s,1,mpi_comm_world,istat,ierr)
do k=1,k1
jt=(k-1)*j0
do j=1,j0
ise(jt+j)=iu(i2,j,k)
isw(jt+j)=iu(3,j,k)
iu(i0,j,k)=ire(jt+j)
iu(1,j,k)=irw(jt+j)
enddo
enddo
if (mylat .ne.ngy1) then
do k=1,k1
it=(k-1)*i0
do i=1,i0
iv(i,j0,k)=irn(it+i)
enddo
enddo
endif
if (mylat .ne. 0) then
do k=1,k1
it=(k-1)*i0
do i=1,i0
iv(i,0,k)=irs(it+i)
enddo
enddo
endif
call mpi_barrier(mpi_comm_world,ierr)
call mpi_sendrecv(isw,nlon,mpi_integer2,m_w,1,ire,nlon,mpi_integer2,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(ise,nlon,mpi_integer2,m_e,1,irw,nlon,mpi_integer2,m_w,1,mpi_comm_world,istat,ierr)
do k=1,k1
jt=(k-1)*j0
do j=1,j0
iu(i0+1,j,k)=ire(jt+j)
iu(0,j,k)=irw(jt+j)
enddo
enddo
do j=2,j1
do i=0,i0+1
iu0(i,j)=iu(i,j,1)
iu(i,j,1)=1
enddo
enddo
do j=0,j0
do i=2,i1
iv0(i,j)=iv(i,j,1)
iv(i,j,1)=1
enddo
enddo
if(mylat.eq.0)then
do i=1,i0
iv(i,0,1)=0
iv(i,1,1)=0
enddo
endif
if(mylat.eq.ngy1)then
do i=1,i0
iv(i,j0,1)=0
iv(i,j1,1)=0
enddo
endif
! write(60) tlz,z,odz,odzw,kb,in,iu,iv,iw,iu0,iv0
! write(60) rmets
call wr8nc_zkb !dog
! ===================================================================
! coupling to 1/8 deg resolution tai and gom grids; one may still run
! entire medina model with coupling to global model, because it uses
! no equatorial shelf.
! this global model, the medina model and norpa! tai grid all fit in
! the memory of a single "fast" processor, but may run less than a
! day per hour of cpu. with one way coupling and separate processors
! for each model, the speed will be almost as fast as the full norpac
! duo grid model.
! ===================================================================
! write rmets data for tai grid (at 150 deg e)
! this covers coupling region with tai grid (j0=488 in tai grid)
! from equator northward
! open(80,file='tai_kbe')
! write(80,80) (kb(601,j),kb(601,j),j=399,642)
! 80 format(40i3)
! open(81,file='gom_kbe')
! write rmets data for gom grid (at 300 deg e)
! this covers coupling region with gom grid (j0=336 in gom grid)
! from 10 deg latitude to 46 deg n
! write(81,80) (kb(1201,j),kb(1201,j),j=439,606)
! write(*,99) (j,(in(i,j,1),i=1,140),j=j0,1,-1)
! 99 format(i3,1x,140i1)
call initfs
end subroutine
! ****************
! initialization *
! ****************
! ----------------------------------------------------------------------
subroutine pre(ax,ay,bb,cx,cy,rinv,rinv1,hh,ie)
! ----------------------------------------------------------------------
real*8 rinv,rinv1,hh,run,rtp,rtemp,rp
dimension ax(i2,j2),ay(i2,j2),bb(i2,j2),cx(i2,j2),cy(i2,0:j2),rinv(ibir,i0,nbg0),rinv1(ibir,i0,nbg1),hh(i0,j0),ie(nb0),rp(i0t*i0t),rtemp(i2t,i0t),rtp(i2t,i2t),ipvt(i2t)
dimension axt(i2,j2),ayt(i2,j2),bbt(i2,j2),cxt(i2,j2),cyt(i2,j2),ct(i2)
do 100 ng=1,ngy0
jl=1
if (mylat .eq. ng-1) then
do 99 j=1,j2
do 99 i=1,i2
axt(i,j)=ax(i,j)
ayt(i,j)=ay(i,j)
bbt(i,j)=bb(i,j)
cxt(i,j)=cx(i,j)
99 cyt(i,j)=cy(i,j)
endif
if (mylat .eq. ng-2) then
do 98 i=1,i2
98 ct(i)=cy(i,j2)
endif
nsd=i2*j2
call mpi_bcast(axt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(ayt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(bbt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(cxt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(cyt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
if (ng .ne. 1)call mpi_bcast(ct,i2,mpi_real8,ng-2,m_comm_lat,ierr)
if (mylat .eq. ng-1) then
do 97 i=1,i2
97 cy(i,0)=ct(i)
endif
do 100 nb=1,nb0
nbg=(ng-1)*nb0+nb
write(*,111) nbg,nbg0
111 format ('processing block #',i3,' out of ',i3,' total evp solver')
jh=ie(nb)
jhp=jh+1
jhm=jh-2
jg=jl+1
do 250 ii=1,ibir
ig=mylat*ibir+ii
idxy=(ig-1)/i2
do 210 j=jl,jhp
do 210 i=1,i0
210 hh(i,j)=0.d0
if (idxy .eq. mylon) then
if (mod(ig,i2) .eq. 0) then
hh(i1,jg)=1.d0
else
hh(mod(ig,i2)+1,jg)=1.d0
endif
endif
call mpi_sendrecv(hh(2,jg),1,mpi_real8,m_w,1,hh(i0,jg),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(hh(i1,jg),1,mpi_real8,m_e,1,hh(1,jg),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
if (nbg .eq. 1) go to 220
if (idxy .eq. mylon) then
if (nb .ne. 1) then
if (mod(ig,i2) .eq. 0) then
ctemp=cyt(i2,jg-2)
else
ctemp=cyt(mod(ig,i2),jg-2)
endif
else
if (mod(ig,i2) .eq. 0) then
ctemp=ct(i2)
else
ctemp=ct(mod(ig,i2))
endif
endif
endif
call mpi_bcast(ctemp,1,mpi_real8,idxy,m_comm_lon,ierr)
do 218 n=1,i0
218 hh(n,jl)=rinv1(ii,n,nbg-1)*ctemp
220 continue
do 226 j=jl,jhm
do 225 i=1,i2
225 hh(i+1,j+2)=-(axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1))/cyt(i,j)
call mpi_sendrecv(hh(2,j+2),1,mpi_real8,m_w,1,hh(i0,j+2),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
226 call mpi_sendrecv(hh(i1,j+2),1,mpi_real8,m_e,1,hh(1,j+2),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
j=jh-1
do 230 i=1,i2
230 rinv(ii,i+1,nbg)=axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1)
if (nbg.eq.nbg0) go to 250
j=ie(nb)
do 240 n=1,i0
240 rinv(ii,n,nbg0)=hh(n,j)
if (idxy .eq. mylon) then
if (mod(ig,i2) .eq. 0) then
hh(i1,jg)=0.d0
else
hh(mod(ig,i2)+1,jg)=0.d0
endif
endif
250 continue
nsend=i2*ibir
call mpi_gather(rinv(1,2,nbg),nsend,mpi_real8,rp(1),nsend,mpi_real8,0,mpi_comm_world,ierr)
if (myid .eq. 0) then
do 497 nx=1,ngx0
ngxt=(nx-1)*ngy0
do 497 ny=1,ngy0
ngt=((ny-1)+ngxt)*ibir*i2
do 497 j=1,i2
nt=ngt+(j-1)*ibir
jt=j+(nx-1)*i2
do 497 i=1,ibir
it=i+(ny-1)*ibir
497 rtemp(it,jt)=rp(nt+i)
do 498 i=1,i2t
do 498 j=1,i2t
498 rtp(j,i)=0.d0
do 499 i=1,i2t
499 rtp(i,i)=1.d0
call dgesv(i2t,i2t,rtemp,i2t,ipvt,rtp,i2t,info)
do 500 i=1,i2t
do 500 j=1,i2t
500 rtemp(j,i)=rtp(i,j)
endif
call mpi_bcast(info,1,mpi_integer,0,mpi_comm_world,ierr)
if (info .ne. 0) then
print*,'preporcessor failed because rinv is singular at nbg =',nbg
stop
endif
nsend=i2t*i2t
call mpi_bcast(rtp,nsend,mpi_real8,0,mpi_comm_world,ierr)
it=mylon*i2
jt=mylat*ibir
jss=1
jff=i0
if (mylon .eq.0) then
do 501 i=1,ibir
501 rinv(i,1,nbg)=rtp(jt+i,i2t)
jss=2
endif
if (mylon .eq. ngx1) then
do 502 i=1,ibir
502 rinv(i,i0,nbg)=rtp(jt+i,1)
jff=i1
endif
do 503 j=jss,jff
do 503 i=1,ibir
503 rinv(i,j,nbg)=rtp(jt+i,it+j-1)
if (nbg.eq.nbg0) return
nsend=i0*ibir
call mpi_allgather(rinv(1,1,nbg0),nsend,mpi_real8,rp,nsend,mpi_real8,m_comm_lat,ierr)
do 504 ny=1,ngy0
ngt=(ny-1)*ibir*i0
do 504 j=1,i0
nt=ngt+(j-1)*ibir
do 504 i=1,ibir
it=i+(ny-1)*ibir
504 rtemp(it,j)=rp(nt+i)
it=ibir*mylat
do 300 i=1,ibir
do 300 j=1,i0
rinv1(i,j,nbg)=0.d0
do 300 k=1,i2t
300 rinv1(i,j,nbg)=rinv1(i,j,nbg)-rtp(it+i,k)*rtemp(k,j)
jl=jh
100 call mpi_barrier(mpi_comm_world,ierr)
end subroutine
! ----------------------------------------------------------------------
subroutine pre2(ax,ay,bb,cx,cy,rinv,rinv1,hh,ie)
! ----------------------------------------------------------------------
real*8 rinv,rinv1,h,run,rtp,rtemp,rp
dimension ax(i2,j2),ay(i2,j2),bb(i2,j2),cx(i2,j2),cy(i2,0:j2),rinv(ibir,i0,nbg0),rinv1(ibir,i0,nbg1),hh(i0,j0),ie(nb0),rp(ig0*ig0),rtemp(ig2,ig0),rtp(ig2,ig2),ipvt(ig2)
dimension axt(i2,j2),ayt(i2,j2),bbt(i2,j2),cxt(i2,j2),cyt(i2,j2),ct(i2)
! if (myid .eq. 0)
do 100 ng=1,ngy0
jl=1
if (mylat .eq. ng-1) then
do 99 j=1,j2
do 99 i=1,i2
axt(i,j)=ax(i,j)
ayt(i,j)=ay(i,j)
bbt(i,j)=bb(i,j)
cxt(i,j)=cx(i,j)
99 cyt(i,j)=cy(i,j)
endif
if (mylat .eq. ng-2) then
do 98 i=1,i2
98 ct(i)=cy(i,j2)
endif
nsd=i2*j2
call mpi_bcast(axt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(ayt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(bbt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(cxt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
call mpi_bcast(cyt,nsd,mpi_real8,ng-1,m_comm_lat,ierr)
if (ng .ne. 1) call mpi_bcast(ct,i2,mpi_real8,ng-2,m_comm_lat,ierr)
if (mylat .eq. ng-1) then
do 97 i=1,i2
97 cy(i,0)=ct(i)
endif
do 100 nb=1,nb0
nbg=(ng-1)*nb0+nb
write(*,111) nbg,nbg0
111 format ('processing block #',i3,' out of ',i3,' total evp solver')
jh=ie(nb)
jhp=jh+1
jhm=jh-2
jg=jl+1
do 250 ii=1,ibir
ig=mylat*ibir+ii
idxy=(ig-1)/i2
do 210 j=jl,jhp
do 210 i=1,i0
210 hh(i,j)=0.d0
if (idxy .eq. mylon) then
if (mod(ig,i2) .eq. 0) then
hh(i1,jg)=1.d0
else
hh(mod(ig,i2)+1,jg)=1.d0
endif
endif
call mpi_sendrecv(hh(2,jg),1,mpi_real8,m_w,1,hh(i0,jg),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(hh(i1,jg),1,mpi_real8,m_e,1,hh(1,jg),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
if (nbg .eq. 1) go to 220
if (idxy .eq. mylon) then
if (nb .ne. 1) then
if (mod(ig,i2) .eq. 0) then
ctemp=cyt(i2,jg-2)
else
ctemp=cyt(mod(ig,i2),jg-2)
endif
else
if (mod(ig,i2) .eq. 0) then
ctemp=ct(i2)
else
ctemp=ct(mod(ig,i2))
endif
endif
endif
call mpi_bcast(ctemp,1,mpi_real8,idxy,m_comm_lon,ierr)
do 218 n=1,i0
218 hh(n,jl)=rinv1(ii,n,nbg-1)*ctemp
220 continue
do 226 j=jl,jhm
do 225 i=1,i2
225 hh(i+1,j+2)=-(axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1))/cyt(i,j)
call mpi_sendrecv(hh(2,j+2),1,mpi_real8,m_w,1,hh(i0,j+2),1,mpi_real8,m_e,1,mpi_comm_world,istat,ierr)
226 call mpi_sendrecv(hh(i1,j+2),1,mpi_real8,m_e,1,hh(1,j+2),1,mpi_real8,m_w,1,mpi_comm_world,istat,ierr)
j=jh-1
do 230 i=1,i2
230 rinv(ii,i+1,nbg)=axt(i,j)*hh(i,j+1)+ayt(i,j)*hh(i+1,j)+bbt(i,j)*hh(i+1,j+1)+cxt(i,j)*hh(i+2,j+1)
if (nbg.eq.nbg0) go to 250
j=ie(nb)
do 240 n=1,i0
240 rinv(ii,n,nbg0)=hh(n,j)
if (idxy .eq. mylon) then
if (mod(ig,i2) .eq. 0) then
hh(i1,jg)=0.d0
else
hh(mod(ig,i2)+1,jg)=0.d0
endif
endif
250 continue
nsend=i2*ibir
call mpi_gather(rinv(1,2,nbg),nsend,mpi_real8,rp(1),nsend,mpi_real8,0,mpi_comm_world,ierr)
if (myid .eq. 0) then
do 497 nx=1,ngx0
ngxt=(nx-1)*ngy0
do 497 ny=1,ngy0
ngt=((ny-1)+ngxt)*ibir*i2
do 497 j=1,i2
nt=ngt+(j-1)*ibir
jt=j+(nx-1)*i2
do 497 i=1,ibir
it=i+(ny-1)*ibir
497 rtemp(it,jt)=rp(nt+i)
do 498 i=1,ig2
do 498 j=1,ig2
498 rtp(j,i)=0.d0
do 499 i=1,ig2
499 rtp(i,i)=1.d0
! write(1233)rtemp
call dgesv(ig2,ig2,rtemp,ig2,ipvt,rtp,ig2,info)
do 500 i=1,ig2
do 500 j=1,ig2
500 rtemp(j,i)=rtp(i,j)
endif
call mpi_bcast(info,1,mpi_integer,0,mpi_comm_world,ierr)
if (info .ne. 0) then
print*,'preporcessor failed because rinv is singular at nbg =',nbg
stop
endif
nsend=ig2*ig2
call mpi_bcast(rtp,nsend,mpi_real8,0,mpi_comm_world,ierr)
it=mylon*i2
jt=mylat*ibir
jss=1
jff=i0
if (mylon .eq.0) then
do 501 i=1,ibir
501 rinv(i,1,nbg)=rtp(jt+i,ig2)
jss=2
endif
if (mylon .eq. ngx1) then
do 502 i=1,ibir
502 rinv(i,i0,nbg)=rtp(jt+i,1)
jff=i1
endif
do 503 j=jss,jff
do 503 i=1,ibir
503 rinv(i,j,nbg)=rtp(jt+i,it+j-1)
if (nbg.eq.nbg0) return
nsend=i0*ibir
call mpi_allgather(rinv(1,1,nbg0),nsend,mpi_real8,rp,nsend,mpi_real8,m_comm_lat,ierr)
do 504 ny=1,ngy0
ngt=(ny-1)*ibir*i0
do 504 j=1,i0
nt=ngt+(j-1)*ibir
do 504 i=1,ibir
it=i+(ny-1)*ibir
504 rtemp(it,j)=rp(nt+i)
it=ibir*mylat
do 300 i=1,ibir
do 300 j=1,i0
rinv1(i,j,nbg)=0.d0
do 300 k=1,ig2
300 rinv1(i,j,nbg)=rinv1(i,j,nbg)-rtp(it+i,k)*rtemp(k,j)
jl=jh
100 call mpi_barrier(mpi_comm_world,ierr)
end subroutine
! ----------------------------------------------------------------------
subroutine sipinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct,alpha)
! ----------------------------------------------------------------------
dimension al(i2,j2),ar(i2,j2),ab(i2,j2),ac(i2,j2),at(i2,0:j2),cl(i2,j2),cr(i2,j2),cb(i2,j2),cc(i2,j2),ct(i2,j2),act(i2,j2)
do 50 j=1,j2
do 50 i=1,i2
50 act(i,j)=ac(i,j)
cl(1,1)=al(1,1)
cb(1,1)=ab(1,1)
cc(1,1)=1/ac(1,1)
ct(1,1)=at(1,1)*cc(1,1)
cr(1,1)=ar(1,1)*cc(1,1)
do 100 i=2,i2
cl(i,1)=al(i,1)
cb(i,1)=ab(i,1)/(1+alpha*cr(i-1,1))
cc(i,1)=1/(act(i,1)+alpha*(cb(i,1)*cr(i-1,1))-cb(i,1)*ct(i-1,1))
ct(i,1)=at(i,1)*cc(i,1)
100 cr(i,1)=(ar(i,1)-alpha*cb(i,1)*ar(i-1,1))*cc(i,1)
do 200 j=2,j2
cl(1,j)=al(1,j)/(1+alpha*ct(1,j-1))
cb(1,j)=ab(1,j)
cc(1,j)=1/(act(1,j)+alpha*(cl(1,j)*ct(1,j-1))-cl(1,j)*cr(1,j-1))
ct(1,j)=(at(1,j)-alpha*cl(1,j)*ct(1,j-1))*cc(1,j)
cr(1,j)=ar(1,j)*cc(1,j)
do 200 i=2,i2
cl(i,j)=al(i,j)/(1+alpha*ct(i,j-1))
cb(i,j)=ab(i,j)/(1+alpha*cr(i-1,j))
cc(i,j)=1/(act(i,j)+alpha*(cl(i,j)*ct(i,j-1)+cb(i,j)*cr(i-1,j))-cl(i,j)*cr(i,j-1)-cb(i,j)*ct(i-1,j))
ct(i,j)=(at(i,j)-alpha*cl(i,j)*ct(i,j-1))*cc(i,j)
200 cr(i,j)=(ar(i,j)-alpha*cb(i,j)*cr(i-1,j))*cc(i,j)
end subroutine
! ----------------------------------------------------------------------
subroutine cholinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
! ----------------------------------------------------------------------
dimension al(i2,j2),ar(i2,j2),ab(i2,j2),ac(i2,j2),at(i2,0:j2),cl(i2,j2),cr(i2,j2),cb(i2,j2),cc(i2,j2),ct(i2,j2),act(i2,j2)
do 50 j=1,j2
do 50 i=1,i2
50 act(i,j)=ac(i,j)
do 100 j=1,j0-3
do 150 i=1,i0-3
cc(i,j)=1/sqrt(act(i,j))
cr(i,j)=ar(i,j)*cc(i,j)
cl(i,j)=cr(i,j)
ct(i,j)=at(i,j)*cc(i,j)
cb(i,j)=ct(i,j)
act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)*cc(i,j)
150 act(i,j+1)=act(i,j+1)-at(i,j)*at(i,j)*cc(i,j)*cc(i,j)
cc(i2,j)=1/sqrt(act(i2,j))
ct(i2,j)=at(i2,j)*cc(i2,j)
cb(i2,j)=ct(i2,j)
100 act(i2,j+1)=act(i2,j+1)-at(i2,j)*at(i2,j)*cc(i2,j)*cc(i2,j)
j=j2
do 200 i=1,i3
cc(i,j)=1/sqrt(act(i,j))
cr(i,j)=ar(i,j)*cc(i,j)
cl(i,j)=cr(i,j)
200 act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)*cc(i,j)
cc(i2,j2)=1/sqrt(ac(i2,j2))
end subroutine
! ----------------------------------------------------------------------
subroutine luinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
! ----------------------------------------------------------------------
dimension al(i2,j2),ar(i2,j2),ab(i2,j2),ac(i2,j2),at(i2,0:j2),cl(i2,j2),cr(i2,j2),cb(i2,j2),cc(i2,j2),ct(i2,j2),act(i2,j2)
do 50 j=1,j2
do 50 i=1,i2
50 act(i,j)=ac(i,j)
do 100 j=1,j0-3
do 150 i=1,i0-3
cc(i,j)=1/act(i,j)
cl(i,j)=al(i,j)*cc(i,j)
cr(i,j)=ar(i,j)
cb(i,j)=ab(i,j)*cc(i,j)
ct(i,j)=at(i,j)
act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)
150 act(i,j+1)=act(i,j+1)-at(i,j)*at(i,j)*cc(i,j)
cc(i2,j)=1/act(i2,j)
cb(i2,j)=ab(i2,j)*cc(i2,j)
ct(i2,j)=at(i2,j)
100 act(i2,j+1)=act(i2,j+1)-at(i2,j)*at(i2,j)*cc(i2,j)
j=j2
do 200 i=1,i3
cc(i,j)=1/act(i,j)
cl(i,j)=al(i,j)*cc(i,j)
cr(i,j)=ar(i,j)
200 act(i+1,j)=act(i+1,j)-ar(i,j)*ar(i,j)*cc(i,j)
cc(i2,j2)=1/ac(i2,j2)
end subroutine
! ---------------------------------------------------------------------
subroutine bounds
! ----------------------------------------------------------------------
real*8 mpd,sltd,pisd,p01,p02,sb8,rb8,dg
dimension ttemp(i0,j0,k1),stemp(i0,j0,k1)
call mpi_type_vector(k1*j0,1,i0,mpi_real8,m_vlon,ierr) !dog
call mpi_type_vector(k1,i0,i0*j0,mpi_real8,m_vlat,ierr)
call mpi_type_commit(m_vlon,ierr)
call mpi_type_commit(m_vlat,ierr)
call rdnc_levitus(1) !dog
! read(70) s1,t1
do 32 n=1,50
write(*,20) n
20 format('iteration #',i3)
do 31 k=1,k1
smin=100.
smax=0.
tmin=100.
tmax=0.
do 30 j=2,j1
do 30 i=2,i1
m=in(i-1,j,k)+in(i+1,j,k)+in(i,j-1,k)+in(i,j+1,k)
if (m*in(i,j,k).eq.0) go to 30
tmp=1./m
ttemp(i,j,k)=tmp*(in(i+1,j,k)*t1(i+1,j,k)+in(i-1,j,k)*t1(i-1,j,k)+in(i,j+1,k)*t1(i,j+1,k)+in(i,j-1,k)*t1(i,j-1,k))
stemp(i,j,k)=tmp*(in(i+1,j,k)*s1(i+1,j,k)+in(i-1,j,k)*s1(i-1,j,k)+in(i,j+1,k)*s1(i,j+1,k)+in(i,j-1,k)*s1(i,j-1,k))
smin=min(smin,stemp(i,j,k))
smax=max(smax,stemp(i,j,k))
tmin=min(tmin,ttemp(i,j,k))
tmax=max(tmax,ttemp(i,j,k))
30 continue
31 continue
do 33 k=1,k1
do 33 j=2,j1
do 33 i=2,i1
t1(i,j,k)=ttemp(i,j,k)
s1(i,j,k)=stemp(i,j,k)
33 continue
call mpi_sendrecv(s1(2,1,1),1,m_vlon,m_w,1,s1(i0,1,1),1,m_vlon,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(s1(i1,1,1),1,m_vlon,m_e,1,s1(1,1,1),1,m_vlon,m_w,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(s1(1,2,1),1,m_vlat,m_s,1,s1(1,j0,1),1,m_vlat,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(s1(1,j1,1),1,m_vlat,m_n,1,s1(1,1,1),1,m_vlat,m_s,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(t1(2,1,1),1,m_vlon,m_w,1,t1(i0,1,1),1,m_vlon,m_e,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(t1(i1,1,1),1,m_vlon,m_e,1,t1(1,1,1),1,m_vlon,m_w,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(t1(1,2,1),1,m_vlat,m_s,1,t1(1,j0,1),1,m_vlat,m_n,1,mpi_comm_world,istat,ierr)
call mpi_sendrecv(t1(1,j1,1),1,m_vlat,m_n,1,t1(1,1,1),1,m_vlat,m_s,1,mpi_comm_world,istat,ierr)
call mpi_barrier(mpi_comm_world,ierr)
32 continue
! levitus initial conditions have been modified
! now, we check for statically unstable points that may exist
50 do 55 k=1,k1
do 55 j=2,j1
do 55 i=2,i1
tmpd=t1(i,j,k)
sltd=s1(i,j,k)
pisd=.001*z(2*k)
! dan wright's full e.o.s.
! rho(i,j,k)=r(tmp,slt,pis)
! function r(t,s,p)
p02=1747.4508988+tmpd*(11.51588-0.046331033*tmpd)-sltd*(3.85429655+0.01353985*tmpd)
p01=pisd+5884.81703666+tmpd*(39.803732+tmpd*(-0.3191477+tmpd*0.0004291133))+2.6126277*sltd
55 d1(i,j,k)=p01/(p02+0.7028423*p01)
do 57 k=2,k1
n=0
do 56 j=2,j1
do 56 i=2,i1
if (in(i,j,k).eq.0.or.d1(i,j,k).ge.d1(i,j,k-1)) go to 56
n=n+1
t1(i,j,k)=t1(i,j,k-1)
s1(i,j,k)=s1(i,j,k-1)
tmpd=t1(i,j,k)
sltd=s1(i,j,k)
pisd=.001*z(2*k)
p02=1747.4508988+tmpd*(11.51588-0.046331033*tmpd)-sltd*(3.85429655+0.01353985*tmpd)
p01=pisd+5884.81703666+tmpd*(39.803732+tmpd*(-0.3191477+tmpd*0.0004291133))+2.6126277*sltd
d1(i,j,k)=p01/(p02+0.7028423*p01)
56 continue
57 write(*,58) n,k
58 format(i6,' unstable points at level',i3)
! periodic
! remark that there doesn't not exchange the b.c. the b! will be exchange in fsglo
! write climatology based initial conditions
! climatology is modified in prep program
! wenien
! open(69,file='data/initial'//grd,form='unformatted')
! write(69) s1,t1
! close(69)
call wr8nc_preini !dog
! 199 rewind 70
! if (.false.) then
! seasonal boundary conditions
do 7000 nm=1,12
write(*,238) nm
238 format('month',i3,' in bounds')
! read in monthly t,s boundary climatologies
call rdnc_levitus(nm) !dog
! read(70) s1,t1
smin=100.
smax=0.
tmin=100.
tmax=0.
do 2400 k=1,k1
do 2400 j=2,j1
do 2400 i=2,i1
s1(i,j,k)=min(max(s1(i,j,k),smax),smin)
2400 t1(i,j,k)=min(max(t1(i,j,k),tmax),tmin)
do 240 j=2,j1
do 240 i=2,i1
tsurf(i-1,j-1,nm)=t1(i,j,1)
240 ssurf(i-1,j-1,nm)=s1(i,j,1)
do 260 k=1,k1
avg1=0.
avg2=0.
areaxy=0.
smax=0.
smin=100.
do 250 j=2,j1
tmp=dx(j)*dy(j)
do 250 i=2,i1
temp=s1(i,j,k)
s1(i,j,k)=max(s1(i,j,k),25.)
if (in(i,j,k).eq.0) go to 250
smax=max(smax,s1(i,j,k))
smin=min(smin,s1(i,j,k))
if (s1(i,j,k).eq.25.) write(*,29) i,j,k,temp,s1(i,j,k)
29 format('fixed bad point i,j,k,s1bad,s1fix=',3i5,2f5.1)
areaxy=areaxy+tmp
avg1=avg1+tmp*t1(i,j,k)
avg2=avg2+tmp*s1(i,j,k)
250 continue
tp=avg1
call mpi_allreduce(tp,avg1,1,mpi_real8,mpi_sum,mpi_comm_world,ierr)
tp=avg2
call mpi_allreduce(tp,avg2,1,mpi_real8,mpi_sum,mpi_comm_world,ierr)
tp=areaxy
call mpi_allreduce(tp,areaxy,1,mpi_real8,mpi_sum,mpi_comm_world,ierr)
write(*,251) nm,k,smin,smax,avg2/areaxy
251 format('nm,k,smin,smax,savg=',2i3,1x,3f5.1)
! if(areaxy.eq.0.) print*,mylon,'_',mylat,'k=',k,'areaxy=',areaxy
txycli(k,nm)=avg1/areaxy
260 sxycli(k,nm)=avg2/areaxy
write(*,261) nm,(sxycli(k,nm),k=1,k1)
261 format(i2/(21f5.1))
! sponge layer climatology
! k=1 is included in ssurf, tssurf
! le
! generally, j0 will not less than 10, so there requires to fix tssp etc.,
! if j0 is less than 10.
if (mylat .eq. 0) then
! do 270 k=1,k1
! do 270 j=1,j2t/18
do 270 k=2,k1
do 270 j=1,10
do 270 i=2,i1
! tssp(i-1,j,k,nm)=t1(i,j+1,k)
! 270 sssp(i-1,j,k,nm)=s1(i,j+1,k)
tssp(i-1,j,k-1,nm)=t1(i,j+1,k)
270 sssp(i-1,j,k-1,nm)=s1(i,j+1,k)
endif
if (mylat .eq. ngy0-1) then
! do 280 k=1,k1
! do 280 j=1,j2t/18
do 280 k=2,k1
do 280 j=1,10
do 280 i=2,i1
! tnsp(i-1,j,k,nm)=t1(i,j0-j,k)
! 280 snsp(i-1,j,k,nm)=s1(i,j0-j,k)
tnsp(i-1,j,k-1,nm)=t1(i,j0-j,k)
280 snsp(i-1,j,k-1,nm)=s1(i,j0-j,k)
endif
7000 continue
! although the following is superfluous when using seasonal climatology,
! (because the model uses linear interpolations anyway),
! we will need 12-month arrays when we use the monthly climatology
! interpolate tsurf, ssurf in time
tmp=1./3.
tem=2./3.
do 600 n=1,10,3
nm=n+3
if (nm.eq.13) nm=1
do 600 j=1,j2
do 600 i=1,i2
tsurf(i,j,n+1)=tsurf(i,j,n)+tmp*(tsurf(i,j,nm)-tsurf(i,j,n))
tsurf(i,j,n+2)=tsurf(i,j,n)+tem*(tsurf(i,j,nm)-tsurf(i,j,n))
ssurf(i,j,n+1)=ssurf(i,j,n)+tmp*(ssurf(i,j,nm)-ssurf(i,j,n))
600 ssurf(i,j,n+2)=ssurf(i,j,n)+tem*(ssurf(i,j,nm)-ssurf(i,j,n))
! interpolate seasonal sponge climatologies
do 650 n=1,10,3
nm=n+3
if (nm.eq.13) nm=1
do 648 k=1,k1
txycli(k,n+1)=txycli(k,n)+tmp*(txycli(k,nm)-txycli(k,n))
txycli(k,n+2)=txycli(k,n)+tem*(txycli(k,nm)-txycli(k,n))
sxycli(k,n+1)=sxycli(k,n)+tmp*(sxycli(k,nm)-sxycli(k,n))
648 sxycli(k,n+2)=sxycli(k,n)+tem*(sxycli(k,nm)-sxycli(k,n))
! do 650 k=1,k1
! do 650 j=1,j2t/18
do 650 k=1,k2
do 650 j=1,10
do 650 i=1,i2
tssp(i,j,k,n+1)=tssp(i,j,k,n)+tmp*(tssp(i,j,k,nm)-tssp(i,j,k,n))
tssp(i,j,k,n+2)=tssp(i,j,k,n)+tem*(tssp(i,j,k,nm)-tssp(i,j,k,n))
tnsp(i,j,k,n+1)=tnsp(i,j,k,n)+tmp*(tnsp(i,j,k,nm)-tnsp(i,j,k,n))
tnsp(i,j,k,n+2)=tnsp(i,j,k,n)+tem*(tnsp(i,j,k,nm)-tnsp(i,j,k,n))
sssp(i,j,k,n+1)=sssp(i,j,k,n)+tmp*(sssp(i,j,k,nm)-sssp(i,j,k,n))
sssp(i,j,k,n+2)=sssp(i,j,k,n)+tem*(sssp(i,j,k,nm)-sssp(i,j,k,n))
snsp(i,j,k,n+1)=snsp(i,j,k,n)+tmp*(snsp(i,j,k,nm)-snsp(i,j,k,n))
650 snsp(i,j,k,n+2)=snsp(i,j,k,n)+tem*(snsp(i,j,k,nm)-snsp(i,j,k,n))
! write surface climatology t,s and specified inflow velocities
! write(62) ssurf,tsurf,sssp,snsp,tssp,tnsp
! write(62) sxycli,txycli
call wr8nc_boundaries !dog
end subroutine
! ----------------------------------------------------------------------
subroutine initfs
! initfs initializes all controlling arrays and derived scalars
!(scalar control parameters are set in block data at start of this file)
! ----------------------------------------------------------------------
dt=172800./daodt
odt=1./dt
orzmx=1./rzmx
ofltw=1.-2.*fltw
mvi=mvfreq*daodt
prn=dm0/de0
! ====================================================
! set background vertical eddy viscosity & diffusivity
! ====================================================
do 90 k=1,k2
! augment molecular viscosity & diffusivity by parameterized synopti! wind
! events, and by breaking (u,v,t,s, near surface only) and non-breaking
! (u,v only, at all rmetss) internal waves.
! synopti! wind forced mixing having 20m e-folding scale.
! bigger momentum mixing due to pressure-xfers
! associated with internal waves that have no counterpart in t,s xfers.
! after yr 7
! 10m scale height for t,s mixing due to breaking internal waves
tmpts=exp(-.001*(z(2*k+1)-z(3)))
!
! deeper tmpuv enhancement is motivated by findings of hurlburt & hogan
! that deep bathymetry and deep transient eddies influences gs path.
! tmpts is not enhanced (layer model has none across layer interface!)
! 1000m scale height for u,v mixing (due to pressure-xfers even without
! internal wave breaking). this enhances deep bathymetry effe!cts on
! upper level flow (e.g., gs separation, kuroshio meanders)
tmpuv=exp(-.00001*(z(2*k+1)-z(3)))
vbk(k)=.01+.05*tmpuv
90 hbk(k)=.002+.02*tmpts
! write(8) vbk,hbk
! close(8)
call wr8nc_winmx !dog
if(myid.eq.0)then
write(*,91) (k,int(.01*z(2*k+1)),vbk(k),hbk(k),k=1,k2)
! write(7,91) (k,int(.01*z(2*k+1)),vbk(k),hbk(k),k=1,k2)
91 format(' K Z VBK HBK'/(i2,i5,2f7.3))
endif
! wenien
! write(9) case,dscrib,dt,tlz,b,g,dm0,de0,dmz0,f,tanphi,rzmx,drag,fltw,daodt,odt,prn,orzmx,ofltw,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv,ktrm,lrstrt,mxit,isav,mxsav,lopen,lmovi,mvi,lfsrf,lsolver,xdeg,lprecon,rawf
! close(9)
call wr8nc_rundata !dog
call mpi_barrier(mpi_comm_world,ierr)
call bounds
if (lsolver.eq.0) then
print*,'preprocessor running without solver setting has finished.'
stop 0
endif
do 250 n=1,nb0
250 ie(n)=min(1+n0*n,j1)
! write(*,251) ie
! 251 format('ie'/(20i4))
if (ie(nb0).gt.ie(nb1).and.j1.lt.ie(nb1)+9) go to 260
! write(17,252) ie(nb1),j1
write(*,252) ie(nb1),j1
252 format('ie(nb1),j1=',2i5,' not properly defined. Fix above do 250 loop and restart.')
260 ie(nb0)=j1
! preprocessor for ellipti! solver used to maintain
! incompressibility for rigid-lid hydrostati! flows
! first, do top layer (all water to regularize evp domain)
tmp=1./odz(1)
do 280 j=1,j2
dzx=odx(j+1)**2*tmp
dzyb=csv(j)*ocs(j+1)*odyv(j)*ody(j+1)*tmp
dzyt=csv(j+1)*ocs(j+1)*odyv(j+1)*ody(j+1)*tmp
do 280 i=1,i2
al(i,j)=dzx
ar(i,j)=dzx
ab(i,j)=dzyb
280 at(i,j)=dzyt
! remaining layers
do 282 k=2,k1
tmp=1./odz(k)
do 282 j=1,j2
dzx=odx(j+1)**2*tmp
dzyb=csv(j)*ocs(j+1)*odyv(j)*ody(j+1)*tmp
dzyt=csv(j+1)*ocs(j+1)*odyv(j+1)*ody(j+1)*tmp
do 282 i=1,i2
al(i,j)=al(i,j)+dzx*iu(i,j+1,k)
ar(i,j)=ar(i,j)+dzx*iu(i+1,j+1,k)
ab(i,j)=ab(i,j)+dzyb*iv(i+1,j,k)
282 at(i,j)=at(i,j)+dzyt*iv(i+1,j+1,k)
if (mylat .eq. 0) then
do 285 i=1,i2
285 ab(i,1)=0.
endif
if (mylat .eq. ngy1) then
do 286 i=1,i2
286 at(i,j2)=0.
endif
gamma_tmp=0.
if (lfsrf .eq. 1) gamma_tmp=1./(0.5*dt**2)
do 287 j=1,j2
do 287 i=1,i2
al(i,j)=-al(i,j)
ar(i,j)=-ar(i,j)
ab(i,j)=-ab(i,j)
at(i,j)=-at(i,j)
287 ac(i,j)=(-al(i,j)-ar(i,j)-ab(i,j)-at(i,j))+gamma_tmp
call mpi_barrier(mpi_comm_world,ierr)
if(lsolver .eq.1) then
if(myid.eq.0) print*,'enter evp preprocessor',myid
call pre(al,ab,ac,ar,at,rinv,rinv1,hh,ie)
! write(99) al,ar,ab,at,ac,rinv,rinv1,ie
! close(99)
call wr8nc_evp1 !dog
if(myid.eq.0) print*,'finished evp preprocessor',myid
elseif(lsolver .eq.2) then
if(myid.eq.0) print*,'enter preconditioner',myid
if(lprecon .le. 1) then
if(myid.eq.0) print*,'sip precondition'
alpha=.5d0
call sipinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct,alpha)
elseif(lprecon .eq. 2) then
if(myid.eq.0) print*,'lu precondition'
call luinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
elseif(lprecon .eq. 3) then
if(myid.eq.0) print*,'chol precondition'
call cholinc(al,ab,ac,ar,at,cl,cb,cc,cr,ct)
else
if(myid.eq.0) print*,'lprecon=',lprecon,' err'
stop 0
endif
! write(99) al,ar,ab,at,ac,cl,cb,cc,cr,ct
! close(99)
call wr8nc_evp2 !dog
endif
print*,'finished preconditioner',myid
end subroutine
subroutine wr8nc_evp1
use const, only: i2t,j2t,i0t,ibir,ngy0,nbg0,nbg1,nb0
use comm, only: myid,mylon,mylat,m_cart
implicit none
!declare pnetcdf var
character :: fn*11,desp*30
integer :: ierr,ncid,nt(3),nd(2),ns(1)
integer :: i2id,j2id,i0id,ibid,ibg0,ibg1,j2idd,inb0
integer :: ial,iar,iab,iac,iat,irin,irin1,ieid
integer(kind=mpi_offset_kind) :: dim_nn(7),leng
integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3),ns_ij(2),nc_ij(2),ns_i,nc_i
!derive the size of each mpiid task,320,288
!output file name
fn = 'data/evp.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) = i0*ngx0
dim_nn(4) = ibir*ngy0
dim_nn(5) = nbg0
dim_nn(6) = nbg1
dim_nn(7) = nb0
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,"i0_global",dim_nn(3),i0id)
ierr=nfmpi_def_dim(ncid,"j2+1_global",dim_nn(2)+ngy0,j2idd)
ierr=nfmpi_def_dim(ncid,"ibir_global",dim_nn(4),ibid)
ierr=nfmpi_def_dim(ncid,"nbg0",dim_nn(5),ibg0)
ierr=nfmpi_def_dim(ncid,"nbg1",dim_nn(6),ibg1)
ierr=nfmpi_def_dim(ncid,"nb0",dim_nn(7),inb0)
!define variables and assign varID
nd(1) = j2id
nd(2) = i2id
ierr=nfmpi_def_var(ncid,"al",nf_double,2,nd,ial)
ierr=nfmpi_def_var(ncid,"ar",nf_double,2,nd,iar)
ierr=nfmpi_def_var(ncid,"ab",nf_double,2,nd,iab)
ierr=nfmpi_def_var(ncid,"ac",nf_double,2,nd,iac)
nd(1) = i2id
nd(2) = j2idd
ierr=nfmpi_def_var(ncid,"at",nf_double,2,nd,iat)
nt(1) = ibid
nt(2) = i0id
nt(3) = ibg0
ierr=nfmpi_def_var(ncid,"rinv",nf_double,3,nt,irin)
nt(1) = ibid
nt(2) = i0id
nt(3) = ibg1
ierr=nfmpi_def_var(ncid,"rinv1",nf_double,3,nt,irin1)
ns(1) = inb0
ierr=nfmpi_def_var(ncid,"ie",nf_int,1,ns,ieid)
desp='input files for lsolver = 1'
leng=len(trim(desp))
ierr=nfmpi_put_att_text(ncid,NF_GLOBAL,'desp',leng,trim(desp))
ierr=nfmpi_enddef(ncid)
!put variables
ns_ij(1) = mylon*i2+1
ns_ij(2) = mylat*j2+1
nc_ij(1) = i2
nc_ij(2) = j2
ierr=nfmpi_put_vara_double_all(ncid,ial,ns_ij,nc_ij,al)
ierr=nfmpi_put_vara_double_all(ncid,iar,ns_ij,nc_ij,ar)
ierr=nfmpi_put_vara_double_all(ncid,iab,ns_ij,nc_ij,ab)
ierr=nfmpi_put_vara_double_all(ncid,iac,ns_ij,nc_ij,ac)
ns_ij(1) = mylon*j2+1
ns_ij(2) = mylat*(j2+1)+1
nc_ij(1) = i2
nc_ij(2) = j2+1
ierr=nfmpi_put_vara_double_all(ncid,iat,ns_ij,nc_ij,at)
ns_ijk(1) = mylat*ibir+1
ns_ijk(2) = mylon*i0+1
ns_ijk(3) = 1
nc_ijk(1) = ibir
nc_ijk(2) = i0
nc_ijk(3) = nbg0
ierr=nfmpi_put_vara_double_all(ncid,irin,ns_ijk,nc_ijk,rinv)
ns_ijk(1) = mylat*ibir+1
ns_ijk(2) = mylon*i0+1
ns_ijk(3) = 1
nc_ijk(1) = ibir
nc_ijk(2) = i0
nc_ijk(3) = nbg1
ierr=nfmpi_put_vara_double_all(ncid,irin1,ns_ijk,nc_ijk,rinv1)
ns_i = 1
nc_i = nb0
ierr=nfmpi_put_vara_int_all(ncid,ieid,ns_i,nc_i,ie)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_evp1
subroutine wr8nc_evp2
use const, only: i2,j2,ngy0,ngx0
use comm, only: myid,mylon,mylat,m_cart
use var, only: grd
implicit none
!declare pnetcdf var
character :: fn*11,desp*120
integer :: ierr,ncid,nd(2)
integer :: i2id,j2id,j2idd
integer :: ial,iar,iab,iac,iat,icl,icb,icc,icr,ict
integer(kind=mpi_offset_kind) :: dim_nn(3),leng
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
!derive the size of each mpiid task
!output file name
fn = 'data/evp.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) = (j2+1)*ngy0
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,"j2+1_global",dim_nn(3),j2idd)
!define variables and assign varID
nd(1) = i2id
nd(2) = j2id
ierr=nfmpi_def_var(ncid,"al",nf_double,2,nd,ial)
ierr=nfmpi_def_var(ncid,"ar",nf_double,2,nd,iar)
ierr=nfmpi_def_var(ncid,"ab",nf_double,2,nd,iab)
ierr=nfmpi_def_var(ncid,"ac",nf_double,2,nd,iac)
nd(2) = j2idd
ierr=nfmpi_def_var(ncid,"at",nf_double,2,nd,iat)
nd(2) = j2id
ierr=nfmpi_def_var(ncid,"cl",nf_double,2,nd,icl)
ierr=nfmpi_def_var(ncid,"cb",nf_double,2,nd,icb)
ierr=nfmpi_def_var(ncid,"cc",nf_double,2,nd,icc)
ierr=nfmpi_def_var(ncid,"cr",nf_double,2,nd,icr)
ierr=nfmpi_def_var(ncid,"ct",nf_double,2,nd,ict)
desp='input files for lsolver = 2, note that the global domain &
is decomposed into ngx0*ngy0 subdomains'
leng=len(trim(desp))
ierr=nfmpi_put_att_text(ncid,NF_GLOBAL,'desp',leng,trim(desp))
ierr=nfmpi_enddef(ncid)
!put variables
ns_ij(1) = mylon*i2+1
ns_ij(2) = mylat*j2+1
nc_ij(1) = i2
nc_ij(2) = j2
ierr=nfmpi_put_vara_double_all(ncid,ial,ns_ij,nc_ij,al)
ierr=nfmpi_put_vara_double_all(ncid,iar,ns_ij,nc_ij,ar)
ierr=nfmpi_put_vara_double_all(ncid,iab,ns_ij,nc_ij,ab)
ierr=nfmpi_put_vara_double_all(ncid,iac,ns_ij,nc_ij,ac)
ns_ij(1) = mylon*i2+1
ns_ij(2) = mylat*(j2+1)+1
nc_ij(1) = i2
nc_ij(2) = j2+1
ierr=nfmpi_put_vara_double_all(ncid,iat,ns_ij,nc_ij,at)
ns_ij(1) = mylon*i2+1
ns_ij(2) = mylat*j2+1
nc_ij(1) = i2
nc_ij(2) = j2
ierr=nfmpi_put_vara_double_all(ncid,icl,ns_ij,nc_ij,cl)
ierr=nfmpi_put_vara_double_all(ncid,icb,ns_ij,nc_ij,cb)
ierr=nfmpi_put_vara_double_all(ncid,icc,ns_ij,nc_ij,cc)
ierr=nfmpi_put_vara_double_all(ncid,icr,ns_ij,nc_ij,cr)
ierr=nfmpi_put_vara_double_all(ncid,ict,ns_ij,nc_ij,ct)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_evp2
subroutine wr8nc_winmx
use const, only: k2
use var, only: z
use comm, only: m_cart
implicit none
!declare pnetcdf var
character :: fn*15,desp*100
integer :: ierr,ncid,ns
integer :: k2id,i
integer :: ivbk,ihbk,idep
integer(kind=mpi_offset_kind) :: dim_nn(1),leng
integer(kind=mpi_offset_kind) :: sns,cns
integer :: dep(k2)
!derive the size of each mpiid task
sns = 1
cns = k2
do i = 1,k2
dep(i) = int(.01*z(2*i+1))
end do
!output file name
fn = 'data/windmix.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) = k2
ierr=nfmpi_def_dim(ncid,"k2",dim_nn(1),k2id)
!define variables and assign varID
ns = k2id
ierr=nfmpi_def_var(ncid,"depth",nf_int,1,ns,idep)
ierr=nfmpi_def_var(ncid,"vbk",nf_double,1,ns,ivbk)
ierr=nfmpi_def_var(ncid,"hbk",nf_double,1,ns,ihbk)
desp='input files for vertical eddy viscosity and diffusivity '
leng=len(trim(desp))
ierr=nfmpi_put_att_text(ncid,NF_GLOBAL,'desp',leng,trim(desp))
ierr=nfmpi_enddef(ncid)
!put variables
ierr=nfmpi_put_vara_int_all(ncid,idep,sns,cns,dep)
ierr=nfmpi_put_vara_double_all(ncid,ivbk,sns,cns,vbk)
ierr=nfmpi_put_vara_double_all(ncid,ihbk,sns,cns,hbk)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_winmx
subroutine wr8nc_rundata
use const, only: j2,j1,j0,i0,j2t,j1t,j0t,i0t,ngx0,ngy0
use comm, only: myid,mylat,mylon,m_cart
implicit none
!declare pnetcdf var
character :: fn*15
integer :: ncid,ierr,ns
integer :: siid,i0id,j0id,j1id,j2id,cid,dsid
integer :: icase,idsp,idt,itlz,ib,ig,idm0,ide0,idmz0
integer :: irzmx,idrag,ifltw,idadt,iodt,iprn,iorz,ioflt,irawf
integer :: iktrm,irst,imxit,isv,imxsv,ilop,ilmov,imvi,ilfs,ilsol,ilpre
integer :: ixdeg,ifcol,iphi,icsv,iocs,idxv,iodx,idyv,iody,iydeg
integer :: iy,iyv,iyvd,ics,iocs2,idx,iodxx,idy,iodyy
integer(kind=mpi_offset_kind) :: dim_nn(5),leng
integer(kind=mpi_offset_kind) :: ns_i,nc_i
!output file name
fn = "data/rundata.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) = 1
dim_nn(2) = i0*ngx0
dim_nn(3) = j0*ngy0
dim_nn(4) = j1*ngy0
dim_nn(5) = j2*ngy0
ierr=nfmpi_def_dim(ncid,"char_case",dim_nn(1)*5,cid)
ierr=nfmpi_def_dim(ncid,"char_desp",dim_nn(1)*63,dsid)
ierr=nfmpi_def_dim(ncid,"pt_val",dim_nn(1),siid)
ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(2),i0id)
ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(3),j0id)
ierr=nfmpi_def_dim(ncid,"j1_global",dim_nn(4),j1id)
ierr=nfmpi_def_dim(ncid,"j2_global",dim_nn(5),j2id)
!define variables and assign varID
ns = cid
ierr=nfmpi_def_var(ncid,"case",nf_char,1,ns,icase)
ns = dsid
ierr=nfmpi_def_var(ncid,"describ",nf_char,1,ns,idsp)
ns = siid
ierr=nfmpi_def_var(ncid,"dt",nf_double,1,ns,idt)
ierr=nfmpi_def_var(ncid,"tlz",nf_double,1,ns,itlz)
ierr=nfmpi_def_var(ncid,"b",nf_double,1,ns,ib)
ierr=nfmpi_def_var(ncid,"g",nf_double,1,ns,ig)
ierr=nfmpi_def_var(ncid,"dm0",nf_double,1,ns,idm0)
ierr=nfmpi_def_var(ncid,"de0",nf_double,1,ns,ide0)
ierr=nfmpi_def_var(ncid,"dmz0",nf_double,1,ns,idmz0)
ierr=nfmpi_def_var(ncid,"rzmx",nf_double,1,ns,irzmx)
ierr=nfmpi_def_var(ncid,"drag",nf_double,1,ns,idrag)
ierr=nfmpi_def_var(ncid,"fltw",nf_double,1,ns,ifltw)
ierr=nfmpi_def_var(ncid,"daodt",nf_double,1,ns,idadt)
ierr=nfmpi_def_var(ncid,"odt",nf_double,1,ns,iodt)
ierr=nfmpi_def_var(ncid,"prn",nf_double,1,ns,iprn)
ierr=nfmpi_def_var(ncid,"orzmx",nf_double,1,ns,iorz)
ierr=nfmpi_def_var(ncid,"ofltw",nf_double,1,ns,ioflt)
ierr=nfmpi_def_var(ncid,"rawf",nf_double,1,ns,irawf)
ierr=nfmpi_def_var(ncid,"ktrm",nf_int,1,ns,iktrm)
ierr=nfmpi_def_var(ncid,"lrstrt",nf_int,1,ns,irst)
ierr=nfmpi_def_var(ncid,"mxit",nf_int,1,ns,imxit)
ierr=nfmpi_def_var(ncid,"isav",nf_int,1,ns,isv)
ierr=nfmpi_def_var(ncid,"mxsav",nf_int,1,ns,imxsv)
ierr=nfmpi_def_var(ncid,"lopen",nf_int,1,ns,ilop)
ierr=nfmpi_def_var(ncid,"lmovi",nf_int,1,ns,ilmov)
ierr=nfmpi_def_var(ncid,"mvi",nf_int,1,ns,imvi)
ierr=nfmpi_def_var(ncid,"lfsrf",nf_int,1,ns,ilfs)
ierr=nfmpi_def_var(ncid,"lsolver",nf_int,1,ns,ilsol)
ierr=nfmpi_def_var(ncid,"lprecon",nf_int,1,ns,ilpre)
ns = j2id
ierr=nfmpi_def_var(ncid,"f",nf_double,1,ns,ifcol)
ierr=nfmpi_def_var(ncid,"tanphi",nf_double,1,ns,iphi)
ns = j1id ! csv,ocsv,dxv,odxv,dyv,odyv
ierr=nfmpi_def_var(ncid,"csv",nf_double,1,ns,icsv)
ierr=nfmpi_def_var(ncid,"ocsv",nf_double,1,ns,iocs)
ierr=nfmpi_def_var(ncid,"dxv",nf_double,1,ns,idxv)
ierr=nfmpi_def_var(ncid,"odxv",nf_double,1,ns,iodx)
ierr=nfmpi_def_var(ncid,"dyv",nf_double,1,ns,idyv)
ierr=nfmpi_def_var(ncid,"odyv",nf_double,1,ns,iody)
ns = j0id ! y,yv,yvdeg,cs,ocs,dx,odx,dy,ody
ierr=nfmpi_def_var(ncid,"y",nf_double,1,ns,iy)
ierr=nfmpi_def_var(ncid,"yv",nf_double,1,ns,iyv)
ierr=nfmpi_def_var(ncid,"yvdeg",nf_double,1,ns,iyvd)
ierr=nfmpi_def_var(ncid,"ydeg",nf_double,1,ns,iydeg)
ierr=nfmpi_def_var(ncid,"cs",nf_double,1,ns,ics)
ierr=nfmpi_def_var(ncid,"ocs",nf_double,1,ns,iocs2)
ierr=nfmpi_def_var(ncid,"dx",nf_double,1,ns,idx)
ierr=nfmpi_def_var(ncid,"odx",nf_double,1,ns,iodxx)
ierr=nfmpi_def_var(ncid,"dy",nf_double,1,ns,idy)
ierr=nfmpi_def_var(ncid,"ody",nf_double,1,ns,iodyy)
ns = i0id
ierr=nfmpi_def_var(ncid,"xdeg",nf_double,1,ns,ixdeg)
ierr=nfmpi_enddef(ncid)
!put variables in collective data mode
ns_i = mylat*j2+1
nc_i = j2
ierr=nfmpi_put_vara_double_all(ncid,ifcol,ns_i,nc_i,f)
ierr=nfmpi_put_vara_double_all(ncid,iphi,ns_i,nc_i,tanphi)
ns_i = mylat*j1+1
nc_i = j1
ierr=nfmpi_put_vara_double_all(ncid,icsv,ns_i,nc_i,csv)
ierr=nfmpi_put_vara_double_all(ncid,iocs,ns_i,nc_i,ocsv)
ierr=nfmpi_put_vara_double_all(ncid,idxv,ns_i,nc_i,dxv)
ierr=nfmpi_put_vara_double_all(ncid,iodx,ns_i,nc_i,odxv)
ierr=nfmpi_put_vara_double_all(ncid,idyv,ns_i,nc_i,dyv)
ierr=nfmpi_put_vara_double_all(ncid,iody,ns_i,nc_i,odyv)
ns_i = mylat*j0+1
nc_i = j0
ierr=nfmpi_put_vara_double_all(ncid,iy,ns_i,nc_i,y)
ierr=nfmpi_put_vara_double_all(ncid,iyv,ns_i,nc_i,yv)
ierr=nfmpi_put_vara_double_all(ncid,iydeg,ns_i,nc_i,ydeg)
ierr=nfmpi_put_vara_double_all(ncid,iyvd,ns_i,nc_i,yvdeg)
ierr=nfmpi_put_vara_double_all(ncid,ics,ns_i,nc_i,cs)
ierr=nfmpi_put_vara_double_all(ncid,iocs2,ns_i,nc_i,ocs)
ierr=nfmpi_put_vara_double_all(ncid,idx,ns_i,nc_i,dx)
ierr=nfmpi_put_vara_double_all(ncid,iodxx,ns_i,nc_i,odx)
ierr=nfmpi_put_vara_double_all(ncid,idy,ns_i,nc_i,dy)
ierr=nfmpi_put_vara_double_all(ncid,iodyy,ns_i,nc_i,ody)
ns_i = mylon*i0+1
nc_i = i0
ierr=nfmpi_put_vara_double_all(ncid,ixdeg,ns_i,nc_i,xdeg)
!put variables in independent data mode
ierr=nfmpi_begin_indep_data(ncid)
ierr=nfmpi_put_var_text(ncid,icase,case)
ierr=nfmpi_put_var_text(ncid,idsp,dscrib)
ierr=nfmpi_put_var_double(ncid,idt,dt)
ierr=nfmpi_put_var_double(ncid,itlz,tlz)
ierr=nfmpi_put_var_double(ncid,ib,b)
ierr=nfmpi_put_var_double(ncid,ig,g)
ierr=nfmpi_put_var_double(ncid,idm0,dm0)
ierr=nfmpi_put_var_double(ncid,ide0,de0)
ierr=nfmpi_put_var_double(ncid,idmz0,dmz0)
ierr=nfmpi_put_var_double(ncid,irzmx,rzmx)
ierr=nfmpi_put_var_double(ncid,idrag,drag)
ierr=nfmpi_put_var_double(ncid,ifltw,fltw)
ierr=nfmpi_put_var_double(ncid,idadt,daodt)
ierr=nfmpi_put_var_double(ncid,iodt,odt)
ierr=nfmpi_put_var_double(ncid,iprn,prn)
ierr=nfmpi_put_var_double(ncid,iorz,orzmx)
ierr=nfmpi_put_var_double(ncid,ioflt,ofltw)
ierr=nfmpi_put_var_double(ncid,irawf,rawf)
ierr=nfmpi_put_var_int(ncid,iktrm,ktrm)
ierr=nfmpi_put_var_int(ncid,irst,lrstrt)
ierr=nfmpi_put_var_int(ncid,imxit,mxit)
ierr=nfmpi_put_var_int(ncid,isv,isav)
ierr=nfmpi_put_var_int(ncid,imxsv,mxsav)
ierr=nfmpi_put_var_int(ncid,ilop,lopen)
ierr=nfmpi_put_var_int(ncid,ilmov,lmovi)
ierr=nfmpi_put_var_int(ncid,imvi,mvi)
ierr=nfmpi_put_var_int(ncid,ilfs,lfsrf)
ierr=nfmpi_put_var_int(ncid,ilsol,lsolver)
ierr=nfmpi_put_var_int(ncid,ilpre,lprecon)
ierr=nfmpi_end_indep_data(ncid)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_rundata
subroutine wr8nc_zkb
use comm, only: myid,mylat,mylon,m_cart
use const, only: i0,j0,k0,i01,k0,k1,ngx0,ngy0
implicit none
!declare pnetcdf var
character :: fn*11
integer :: ncid,ierr,ns,nd(2),nt(3)
integer :: siid,k0id,k1id,kkid,i01id,i0id,j0id,i11id,j11id
integer :: itlz,iz,iodz,iodzw,ikb,irmt,iiv0,iiu0
integer :: iiu,iiv,iiw,iin
integer(kind=mpi_offset_kind) :: dim_nn(9),ns_ijk(3),nc_ijk(3),ns_ij(2),nc_ij(2)
integer(kind=mpi_offset_kind) :: ns_i,nc_i
!output file name
fn = 'data/zkb.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) = 1
dim_nn(2) = k0
dim_nn(3) = k1
dim_nn(4) = k0+k1
dim_nn(5) = i01*ngx0
dim_nn(6) = i0*ngx0
dim_nn(7) = j0*ngy0
dim_nn(8) = (i01+1)*ngx0
dim_nn(9) = (j0+1)*ngy0
ierr=nfmpi_def_dim(ncid,"pt_val",dim_nn(1),siid)
ierr=nfmpi_def_dim(ncid,"k0",dim_nn(2),k0id)
ierr=nfmpi_def_dim(ncid,"k1",dim_nn(3),k1id)
ierr=nfmpi_def_dim(ncid,"k0+k1",dim_nn(4),kkid)
ierr=nfmpi_def_dim(ncid,"i01_global",dim_nn(5),i01id)
ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(6),i0id)
ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(7),j0id)
ierr=nfmpi_def_dim(ncid,"i01+1_global",dim_nn(8),i11id)
ierr=nfmpi_def_dim(ncid,"j0+1_global",dim_nn(9),j11id)
!define variables and assign varID
ns = siid
ierr=nfmpi_def_var(ncid,"tlz",nf_double,1,ns,itlz)
ns = kkid
ierr=nfmpi_def_var(ncid,"z",nf_double,1,ns,iz)
ns = k1id
ierr=nfmpi_def_var(ncid,"odz",nf_double,1,ns,iodz)
ns = k0id
ierr=nfmpi_def_var(ncid,"odzw",nf_double,1,ns,iodzw)
nd(1) = i0id
nd(2) = j0id
ierr=nfmpi_def_var(ncid,"kb",nf_short,2,nd,ikb)
ierr=nfmpi_def_var(ncid,"rmets",nf_double,2,nd,irmt)
nd(1) = i0id
nd(2) = j11id
ierr=nfmpi_def_var(ncid,"iv0",nf_short,2,nd,iiv0)
nd(1) = i11id
nd(2) = j0id
ierr=nfmpi_def_var(ncid,"iu0",nf_short,2,nd,iiu0)
nt(1) = i0id
nt(2) = j0id
nt(3) = k0id
ierr=nfmpi_def_var(ncid,"iw",nf_short,3,nt,iiw)! (i0,j0,k0)
nt(3) = k1id
ierr=nfmpi_def_var(ncid,"in",nf_short,3,nt,iin)! (i0,j0,k1)
nt(2) = j11id
ierr=nfmpi_def_var(ncid,"iv",nf_short,3,nt,iiv)! (i0,0:j0,k1)
nt(1) = i11id
nt(2) = j0id
ierr=nfmpi_def_var(ncid,"iu",nf_short,3,nt,iiu)
ierr=nfmpi_enddef(ncid)
!put variables
ns_i = 1
nc_i = 1
ierr=nfmpi_put_vara_double_all(ncid,itlz,ns_i,nc_i,tlz)
nc_i = k0+k1
ierr=nfmpi_put_vara_double_all(ncid,iz,ns_i,nc_i,z)
nc_i = k1
ierr=nfmpi_put_vara_double_all(ncid,iodz,ns_i,nc_i,odz)
nc_i = k0
ierr=nfmpi_put_vara_double_all(ncid,iodzw,ns_i,nc_i,odzw)
ns_ij(1) = mylon*i0+1
ns_ij(2) = mylat*j0+1
nc_ij(1) = i0
nc_ij(2) = j0
ierr=nfmpi_put_vara_int2_all(ncid,ikb,ns_ij,nc_ij,kb)
ierr=nfmpi_put_vara_double_all(ncid,irmt,ns_ij,nc_ij,rmets)
ns_ij(2) = mylat*(j0+1)+1
nc_ij(2) = j0+1
ierr=nfmpi_put_vara_int2_all(ncid,iiv0,ns_ij,nc_ij,iv0)
ns_ij(1) = mylon*(i01+1)+1
ns_ij(2) = mylat*j0+1
nc_ij(1) = i01+1
nc_ij(2) = j0
ierr=nfmpi_put_vara_int2_all(ncid,iiu0,ns_ij,nc_ij,iu0)
ns_ijk(1) = mylon*i0+1
ns_ijk(2) = mylat*j0+1
ns_ijk(3) = 1
nc_ijk(1) = i0
nc_ijk(2) = j0
nc_ijk(3) = k0
ierr=nfmpi_put_vara_int2_all(ncid,iiw,ns_ijk,nc_ijk,iw)
nc_ijk(3) = k1
ierr=nfmpi_put_vara_int2_all(ncid,iin,ns_ijk,nc_ijk,in)
ns_ijk(2) = mylat*(j0+1)+1
nc_ijk(2) = j0+1
ierr=nfmpi_put_vara_int2_all(ncid,iiv,ns_ijk,nc_ijk,iv)
ns_ijk(1) = mylon*(i01+1)+1
ns_ijk(2) = mylat*j0+1
nc_ijk(1) = i01+1
nc_ijk(2) = j0
ierr=nfmpi_put_vara_int2_all(ncid,iiu,ns_ijk,nc_ijk,iu)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_zkb
subroutine wr8nc_boundaries
use const, only: i2,j2,ngx0,ngy0,k1,k2
use comm, only: m_cart,mylat,mylon
implicit none
!declare pnetcdf var
character :: fn*18
integer :: ierr,ncid,nd(2),nt(3),nf(4)
integer :: i2id,j2id,spid,moid,k1id,k2id
integer :: isxy,itxy,itsu,issu,isnp,issp,itnp,itsp
integer(kind=mpi_offset_kind) :: dim_nn(6)
integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4),ns_ijk(3),nc_ijk(3)
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
!output file name
fn = 'data/boundaries.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) = 10*ngy0
dim_nn(4) = 12
dim_nn(5) = k1
dim_nn(6) = k2
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,"j0-sponge_global",dim_nn(3),spid)
ierr=nfmpi_def_dim(ncid,"month",dim_nn(4),moid)
ierr=nfmpi_def_dim(ncid,"k1",dim_nn(5),k1id)
ierr=nfmpi_def_dim(ncid,"k2",dim_nn(6),k2id)
!define variables and assign varID
nd(1) = k1id
nd(2) = moid
ierr=nfmpi_def_var(ncid,"sxycli",nf_double,2,nd,isxy)
ierr=nfmpi_def_var(ncid,"txycli",nf_double,2,nd,itxy)
nt(1) = i2id
nt(2) = j2id
nt(3) = moid
ierr=nfmpi_def_var(ncid,"tsurf",nf_double,3,nt,itsu)
ierr=nfmpi_def_var(ncid,"ssurf",nf_double,3,nt,issu)
nf(1) = i2id
nf(2) = spid
nf(3) = k2id
nf(4) = moid
ierr=nfmpi_def_var(ncid,"snsp",nf_double,4,nf,isnp)
ierr=nfmpi_def_var(ncid,"sssp",nf_double,4,nf,issp)
ierr=nfmpi_def_var(ncid,"tnsp",nf_double,4,nf,itnp)
ierr=nfmpi_def_var(ncid,"tssp",nf_double,4,nf,itsp)
ierr=nfmpi_enddef(ncid)
!put variables
ns_ij(1) = 1
ns_ij(2) = 1
nc_ij(1) = k1
nc_ij(2) = 12
ierr=nfmpi_put_vara_double_all(ncid,isxy,ns_ij,nc_ij,sxycli)
ierr=nfmpi_put_vara_double_all(ncid,itxy,ns_ij,nc_ij,txycli)
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,issu,ns_ijk,nc_ijk,ssurf)
ierr=nfmpi_put_vara_double_all(ncid,itsu,ns_ijk,nc_ijk,tsurf)
ns_ijkt(1) = mylon*i2+1
ns_ijkt(2) = mylat*10+1
ns_ijkt(3) = 1
ns_ijkt(4) = 1
nc_ijkt(1) = i2
nc_ijkt(2) = 10
nc_ijkt(3) = k2
nc_ijkt(4) = 12
ierr=nfmpi_put_vara_double_all(ncid,isnp,ns_ijkt,nc_ijkt,snsp)
ierr=nfmpi_put_vara_double_all(ncid,issp,ns_ijkt,nc_ijkt,sssp)
ierr=nfmpi_put_vara_double_all(ncid,itnp,ns_ijkt,nc_ijkt,tnsp)
ierr=nfmpi_put_vara_double_all(ncid,itnp,ns_ijkt,nc_ijkt,tssp)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_boundaries
subroutine wr8nc_preini
use const, only: i0,j0,k1,ngx0,ngy0
use comm, only: m_cart,mylat,mylon
implicit none
!declare pnetcdf var
character :: fn*14
integer :: ierr,ncid,nt(3)
integer :: i0id,j0id,k1id
integer :: isd,itd
integer(kind=mpi_offset_kind) :: dim_nn(3)
integer(kind=mpi_offset_kind) :: ns_ijk(3),nc_ijk(3)
!output file name
fn = 'data/preini.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
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)
!define variables and assign varID
nt(1) = i0id
nt(2) = j0id
nt(3) = k1id
ierr=nfmpi_def_var(ncid,"s1",nf_double,3,nt,isd)
ierr=nfmpi_def_var(ncid,"t1",nf_double,3,nt,itd)
ierr=nfmpi_enddef(ncid)
!put variables
ns_ijk(1) = mylon*i0+1
ns_ijk(2) = mylat*j0+1
ns_ijk(3) = 1
nc_ijk(1) = i0
nc_ijk(2) = j0
nc_ijk(3) = k1
ierr=nfmpi_put_vara_double_all(ncid,isd,ns_ijk,nc_ijk,s1)
ierr=nfmpi_put_vara_double_all(ncid,itd,ns_ijk,nc_ijk,t1)
ierr=nfmpi_close(ncid)
end subroutine wr8nc_preini
subroutine rdnc_levitus(nm)
use comm, only: m_cart,myid,mylat,mylon
use const, only: i0,j0,k1,ngy0,ngx0
implicit none
character :: fn*20
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ijkt(4),nc_ijkt(4)
integer,intent(in) :: nm
fn = 'data/annualevitus.nc'
ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
ns_ijkt(1) = mylon*i0+1
ns_ijkt(2) = mylat*j0+1
ns_ijkt(3) = 1
ns_ijkt(4) = nm
nc_ijkt(1) = i0
nc_ijkt(2) = j0
nc_ijkt(3) = k1
nc_ijkt(4) = 1
ierr=nfmpi_inq_varid(ncid,'s1_data',itmp)
ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijkt,nc_ijkt,s1)
ierr=nfmpi_inq_varid(ncid,'t1_data',itmp)
ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ijkt,nc_ijkt,t1)
ierr=nfmpi_close(ncid)
endsubroutine rdnc_levitus
subroutine rdnc_depth(rmets)
use comm, only: m_cart,myid,mylat,mylon
use const, only: i0,j0
implicit none
character :: fn*14
integer :: ncid,ierr,itmp
integer(kind=mpi_offset_kind) :: ns_ij(2),nc_ij(2)
real*8,intent(inout) :: rmets(i0,j0)
fn = 'data/depth.nc'
ierr=nfmpi_open(m_cart,fn,nf_nowrite,mpi_info_null,ncid)
ierr=nfmpi_inq_varid(ncid,'rmets',itmp)
ns_ij(1) = mylon*i0+1
ns_ij(2) = mylat*j0+1
nc_ij(1) = i0
nc_ij(2) = j0
ierr=nfmpi_get_vara_double_all(ncid,itmp,ns_ij,nc_ij,rmets)
ierr=nfmpi_close(ncid)
endsubroutine rdnc_depth
subroutine rdnc_yzgrid2(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_yzgrid2
end module