PROGRAM HYCOM_PROFILE_NCODA_INC IMPLICIT NONE C C hycom_profile_ncoda_inc - Usage: hycom_profile_ncoda_inc blkdat.input C C apply ncoda_archv_inc to a HYCOM profile. C C blkdat.input is a ncoda_archv_inc parameter file, modified so C that there is a single ncoda file for t,s,u,v,p C C this version for "serial" Unix systems. C C Alan J. Wallcraft, Naval Research Laboratory, December 2015. C INTEGER IARGC INTEGER NARG CHARACTER*240 CARG C CHARACTER*240 flnm_b,CLINE real skip integer ki,ios c c --- from mod_plot c integer :: & ii,jj,kk,kkmax real, dimension (1,1,99) :: & u,v,ke,temp,saln,th3d,dw,dp,dpsd,p c real, dimension (1,1) :: & ubaro,vbaro,pbaro,kebaro, montg,srfht,steric,dpbl,dpmixl, & tmix,smix,thmix,umix,vmix,kemix, & depths,coast,plon,plat,qlon,qlat,ulon,ulat,vlon,vlat, & surflx,salflx,surtx,surty, ttrend,strend,emnp, & covice,thkice,temice, & scux,scvx,scpx,scuy,scvy,scpy,pang, & field c real, dimension (1) :: & theta c integer, dimension (1,1) :: & ip,iq,iu,iv c c --- from ncoda_archv_inc C real, parameter :: flag = 2.0**100 !data void marker c character label*81,text*18,flnm_i*240,flnm_o*240,flnm_d*240, & flnm_t*240,flnm_s*240,flnm_p*240,flnm_m*240, & flnm_u*240,flnm_v*240,flnm_c*240,flnm_z*240, & mark*3 logical initl,trcout,lsteric,icegln,larctic,vsigma logical lprsout,lpncz c integer artype,iexpt,iversn,yrflag,iweight,mntype integer i,ia,ibad,ih,j,ja,iit,jjt, & k,k2,kkin,kkout,kmld,l,newtop integer intflg,isoflg,nddflg,itest,jtest integer kncoda,kncodz integer nhybrd,nsigma real dp00,dp00x,dp00f,dp00i,ds00,ds00x,ds00f,dp0ij, & dp0k(99),dp0kf,dpm,dpms,isotop, & ds0k(99),ds0kf,dsm,dsms,dssk(99),dp0cum(99+1) real u1(99),v1(99),t1(99),s1(99),r1(99),p1(0:99), & uz(99),vz(99),tz(99),sz(99),rz(99),pz(0:99), & to(99),so(99),ro(99),po(0:99), & ti(99),si(99),ri(99), & zz(99),zi(0:99),rl(99), & pnc(99),pns(99) real sigma(99),thk(99),thbase,deniso,thnthk, & salmin,frzoff,frzinc,frcrng, & deqthu,deqthv,onem,qonem,thref,q, & hicemn,mldij,qk,pk,pnck,pnk,pnimin,pnmax, & vzero,uzero real hmina,hmaxa,hminn(99),hmaxn(99) real epsil,dpbot,dpmid,dpnew,dptop,dpinc double precision time3(3),time,year,mass_h,mass_n real prsout(1,1,99) real pout(1,1,99),theta3(1,1,99) real tncoda(1,1,99),sncoda(1,1,99),pncoda(1,1,99) real uinc(1,1,99),vinc(1,1,99) real uncoda(1,1,99),vncoda(1,1,99) real cncoda(1,1), mncoda(1,1) real tincmx(1,1) real pij(99),mldlay(1,1),work(1,1) integer incoda(1,1) c data trcout/.false./ ! must be .false. (no tracer remapping) data initl /.true. / c REAL*4 SIG_V,SIGLOC_V,SOFSIG_V,TOFSIG_V c integer sigver real g integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c integer thflag common/th/thflag c i0 = 0 j0 = 0 nbdy = 0 mnproc = 1 lp = 6 g = 9.806 uoff = 0 ii = 1 jj = 1 C C READ ARGUMENTS. C NARG = IARGC() C IF (NARG.EQ.1) THEN CALL GETARG(1,flnm_b) ELSE WRITE(6,*) + 'Usage: hycom_profile_ncoda_inc blkdat.input' CALL EXIT(1) ENDIF C OPEN(UNIT=uoff+99, FILE=flnm_b, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(flnm_b) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF c epsil = 1.0e-4 ! very small density increment thref = 1.0e-3 onem = 9806.0 ! g/thref qonem = 1.0/onem c c --- 'flnm_i' = name of original archive file c --- 'flnm_o' = name of target archive file c --- 'intflg' = vertical interpolation flag, must be 0=T&S. c --- 'isoflg' = Preserve isopycnal layer flag (0=n(def),1=y,2=y&layT,3=y&isoT) c --- 'nddflg' = ncoda density displacement flag (0=no,1=p,2=p+T&S,3=T&S) c --- 'iexpt ' = experiment number x10 (000=from archive file) c --- 'yrflag' = days in year flag (0=360J16,1=366J16,2=366J01,3=actual) c --- 'idm ' = longitudinal array size c --- 'jdm ' = latitudinal array size c --- 'itest ' = longitudinal test point (optional, default 0) c --- 'jtest ' = latitudinal test point (optional, default 0) c --- 'kdm ' = number of layers c read (uoff+99,'(a)') flnm_i if(mnproc.eq.1)then write (lp,'(2a)') ' input file: ',trim(flnm_i) call flush(lp) endif ! end if(mnproc.eq.1) block read (uoff+99,'(a)') flnm_o if(mnproc.eq.1)then write (lp,'(2a)') 'output file: ',trim(flnm_o) write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call blkini(intflg,'intflg') call blkini(isoflg,'isoflg') call blkini(nddflg,'nddflg') call blkini(iexpt ,'iexpt ') *** call blkini(yrflag,'yrflag') call blkini(iit, 'idm ') call blkini(jjt, 'jdm ') call blkini2(i,j, 'itest ','kdm ') !read itest or kdm if (j.eq.1) then itest = i call blkini(jtest, 'jtest ') call blkini(kkin, 'kdm ') else itest = 0 jtest = 0 kkin = i endif kkout = kkin ****** if (iit.ne.itdm .or. jjt.ne.jtdm) then ****** if(mnproc.eq.1)then ****** write(lp,*) ****** write(lp,*) 'error - wrong idm or jdm (should be:', ****** & itdm,jtdm,')' ****** write(lp,*) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call xcstop('ncoda_archv_inc') ****** stop ****** endif if (intflg.ne.0) then if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - must have intflg=0 for T&S increments' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (isoflg.eq.3) then if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - isoflg=3 not possible with T&S increments' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif c c --- 'nhybrd' = number of hybrid levels (0=all isopycnal) c --- 'nsigma' = number of sigma levels (nhybrd-nsigma z-levels) c --- 'dp00' = deep z-level spacing minimum thickness (m) c --- 'dp00x' = deep z-level spacing maximum thickness (m) c --- 'dp00f' = deep z-level spacing stretching factor (1.0=const.z) c --- 'ds00' = shallow z-level spacing minimum thickness (m) c --- 'ds00x' = shallow z-level spacing maximum thickness (m) c --- 'ds00f' = shallow z-level spacing stretching factor (1.0=const.z) c --- 'dp00i' = deep iso-pycnal spacing minimum thickness (m) c --- 'isotop' = shallowest depth for isopycnal layers (m) c c --- the above specifies a vertical coord. that is isopycnal or: c --- z in deep water, based on dp00,dp00x,dp00f c --- z in shallow water, based on ds00,ds00x,ds00f and nsigma c --- sigma between them, based on ds00,ds00x,ds00f and nsigma c --- for z-only set nsigma=0 (and ds00,ds00x,ds00f=dp00,dp00x,dp00f) c --- for sigma-z (shallow-deep) use a very small ds00 c --- (pure sigma-z also has ds00f=dp00f and ds00x=dp00x*ds00/dp00) c --- for z-sigma (shallow-deep) use a very large dp00 (not recommended) c --- for sigma-only set nsigma=kdm, dp00 large, and ds00 small c c --- or, in place of 'dp00','dp00x','dp00f','ds00','ds00x','ds00f' specify: c --- 'dp0k ' = layer k deep z-level spacing minimum thickness (m) c --- k=1,kdm; dp0k must be zero for k>nhybrd c --- 'ds0k ' = layer k shallow z-level spacing minimum thickness (m) c --- k=1,nsigma c c --- away from the surface, the minimum layer thickness is dp00i. c call blkini(nhybrd,'nhybrd') call blkini(nsigma,'nsigma') call blkinr2(dp00,k, 'dp00 ','(a6," =",f10.4," m")', & 'dp0k ','(a6," =",f10.4," m")' ) if (k.eq.1) then !dp00 call blkinr(dp00x, 'dp00x ','(a6," =",f10.4," m")') call blkinr(dp00f, 'dp00f ','(a6," =",f10.4," ")') call blkinr(ds00, 'ds00 ','(a6," =",f10.4," m")') call blkinr(ds00x, 'ds00x ','(a6," =",f10.4," m")') call blkinr(ds00f, 'ds00f ','(a6," =",f10.4," ")') else !dp0k dp0k(1) = dp00 dp00 = -1.0 !signal that dp00 is not input do k=2,kkout call blkinr(dp0k(k), 'dp0k ','(a6," =",f10.4," m")') c if (k.gt.nhybrd .and. dp0k(k).ne.0.0) then if(mnproc.eq.1)then write(lp,'(/ a,i3 /)') & 'error - dp0k must be zero for k>nhybrd' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif !k>nhybrd&dp0k(k)!=0 enddo !k do k=1,nsigma call blkinr(ds0k(k), 'ds0k ','(a6," =",f10.4," m")') enddo !k endif !dp00:dp0k c call blkinr(dp00i, 'dp00i ','(a6," =",f10.4," m")') call blkinr(isotop,'isotop','(a6," =",f10.4," m")') c if (nhybrd.gt.kkout) then if(mnproc.eq.1)then write(lp,'(/ a,i3 /)') & 'error - maximum nhybrd is kdmnew =',kkout call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (nsigma.gt.nhybrd) then if(mnproc.eq.1)then write(lp,'(/ a,i3 /)') & 'error - maximum nsigma is nhybrd =',nhybrd call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (dp00.ge.0.0) then if (dp00f.lt.1.0) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - must have dp00f>=1.0' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (dp00f.eq.1.0 .and. dp00.ne.dp00x) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - must have dp00x==dp00 for dp00f==1.0' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (dp00.gt.dp00x) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - dp00x must be at least dp00' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (ds00.gt.dp00 .or. ds00x.gt.dp00x .or. ds00f.gt.dp00f) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - must have ds00,ds00x,ds00f <= dp00,dp00x,dp00f' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (ds00.le.0.0) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - must have ds00>0.0' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (ds00f.lt.1.0) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - must have ds00f>=1.0' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (ds00f.eq.1.0 .and. ds00.ne.ds00x) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - must have ds00x==ds00 for ds00f==1.0' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (ds00.gt.ds00x) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - ds00x must be at least ds00' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif endif !dp00 used c c --- 'salmin' = minimum salinity (optional, default 0.0) c --- 'frzoff' = minimum temperature is freezing+frzoff (optional, default -9) c --- 'frcrng' = maximum absolute T increment is frcrng*vertical range c (optional, default 0.0 which turns off the limiter) c --- 'deniso' = isopycnal if layer is within deniso of target density c (for isoflg>0, large to recover previous behaviour) c --- 'thnthk' = minimum ratio of thin to thick isopycnal layers (0 to 1) c (for isoflg>0, zero to recover previous behaviour) c --- 'thbase' = new reference density (sigma units) c call blkinr9(q,k, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")', & 'salmin','("blkinr: ",a6," =",f11.4," psu")', & 'frzoff','("blkinr: ",a6," =",f11.4," degC")', & 'frcrng','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")') if (k.eq.1) then !deniso deniso = q salmin = 0.0 !default frzoff = -9.0 !default frcrng = 0.0 !default, turn off frcrng elseif (k.eq.2) then !salmin salmin = q call blkinr9(q,k, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")', & 'frzoff','("blkinr: ",a6," =",f11.4," degC")', & 'frcrng','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & 'XXXXXX','("blkinr: ",a6," =",f11.4," ")') if (k.eq.1) then !deniso frzoff = -9.0 !default frcrng = 0.0 !default, turn off frcrng deniso = q elseif (k.eq.2) then !frzoff frzoff = q call blkinr2(q,k, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")', & 'frcrng','("blkinr: ",a6," =",f11.4," ")') if (k.eq.1) then !deniso frcrng = 0.0 !default, turn off frcrng deniso = q elseif (k.eq.2) then frcrng = q call blkinr(deniso, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")') endif !deniso,frcrng else !frcrng frzoff = -9.0 !default frcrng = q call blkinr(deniso, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")') endif !deniso,frzoff,frcrng elseif (k.eq.3) then !frzoff salmin = 0.0 !default frzoff = q call blkinr2(q,k, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")', & 'frcrng','("blkinr: ",a6," =",f11.4," ")') if (k.eq.1) then !deniso frcrng = 0.0 !default, turn off frcrng deniso = q elseif (k.eq.2) then frcrng = q call blkinr(deniso, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")') endif !deniso,frcrng elseif (k.eq.4) then !frcrng salmin = 0.0 !default frzoff = -9.0 !default frcrng = q call blkinr(deniso, & 'deniso','("blkinr: ",a6," =",f11.4," kg/m^3")') else if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - unknown input (not salmin,frzoff,frcrng,deniso)' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif !salmin,frzoff,deniso call blkinr(thnthk, & 'thnthk','("blkinr: ",a6," =",f11.4," ")') call blkinr(thbase, & 'thbase','("blkinr: ",a6," =",f11.4," sig")') c c --- 'vsigma' = spacially varying isopycnal target densities (0=F,1=T) c call blkinl(vsigma,'vsigma') c c --- target layer densities (sigma units) c if(mnproc.eq.1)then write(lp,*) endif ! end if(mnproc.eq.1) block do k=1,kkout call blkinr(sigma(k), & 'sigma ','("blkinr: ",a6," =",f11.4," sig")') thk(k) = sigma(k) - thbase c if (k.gt.1) then if (sigma(k).le.sigma(k-1)) then if(mnproc.eq.1)then write(lp,'(/ a /)') & 'error - sigma is not stabally stratified' call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif endif enddo c c --- 'hicemn' = (ENLN) minimum ice thickness (m) c call blkinr(hicemn, & 'hicemn','("blkinr: ",a6," =",f11.4," m")') c dp00i =dp00i *onem isotop =isotop*onem c c --- calculate dp0k and ds0k? if (dp00.lt.0.0) then c --- dp0k and ds0k already input dpms = 0.0 do k=1,kkout dp0k(k) = dp0k(k)*onem dpm = dp0k(k) dpms = dpms + dpm if(mnproc.eq.1)then write(lp,135) k,dp0k(k)*qonem,dpm*qonem,dpms*qonem call flush(lp) endif ! end if(mnproc.eq.1) block enddo !k dsms = 0.0 do k=1,nsigma ds0k(k) = ds0k(k)*onem dsm = ds0k(k) dsms = dsms + dsm if(mnproc.eq.1)then write(lp,130) k,ds0k(k)*qonem,dsm*qonem,dsms*qonem call flush(lp) endif ! end if(mnproc.eq.1) block enddo !k if(mnproc.eq.1)then write(lp,*) endif ! end if(mnproc.eq.1) block else c --- calculate dp0k and ds0k c c --- logorithmic k-dependence of dp0 (deep z''s) dp00 =dp00 *onem dp00x =dp00x*onem dp0k(1)=dp00 dpm =dp0k(1) dpms =dpm if(mnproc.eq.1)then write(lp,*) write(lp,135) 1,dp0k(1)*qonem,dpm*qonem,dpms*qonem call flush(lp) endif ! end if(mnproc.eq.1) block 135 format('dp0k(',i2,') =',f7.2,' m', & ' thkns =',f7.2,' m', & ' depth =',f8.2,' m') c dp0kf=1.0 do k=2,kkout dp0kf=dp0kf*dp00f if (k.le.nhybrd) then dp0k(k)=min(dp00*dp0kf,dp00x) else dp0k(k)=0.0 endif dpm = dp0k(k) dpms = dpms + dpm if(mnproc.eq.1)then write(lp,135) k,dp0k(k)*qonem,dpm*qonem,dpms*qonem call flush(lp) endif ! end if(mnproc.eq.1) block enddo c c --- logorithmic k-dependence of ds0 (shallow z-s) ds00 =ds00 *onem ds00x =ds00x*onem ds0k(1)=ds00 dsm =ds0k(1) dsms =dsm if(mnproc.eq.1)then write(lp,*) write(lp,130) 1,ds0k(1)*qonem,dsm*qonem,dsms*qonem endif ! end if(mnproc.eq.1) block 130 format('ds0k(',i2,') =',f7.2,' m', & ' thkns =',f7.2,' m', & ' depth =',f8.2,' m') call flush(lp) c ds0kf=1.0 do k=2,nsigma ds0kf=ds0kf*ds00f ds0k(k)=min(ds00*ds0kf,ds00x) dsm = ds0k(k) dsms = dsms + dsm if(mnproc.eq.1)then write(lp,130) k,ds0k(k)*qonem,dsm*qonem,dsms*qonem call flush(lp) endif ! end if(mnproc.eq.1) block enddo if(mnproc.eq.1)then write(lp,*) endif ! end if(mnproc.eq.1) block endif !input:calculate dp0k,ds0k c c --- sigma-depth scale factors do k=1,nsigma dssk(k)=ds0k(k)/dsms ! fraction of depths in sigma layer k enddo do k= nsigma+1,kkout ds0k(k)=dp0k(k) dssk(k)=0.0 ! these layers are zero in sigma mode enddo ******c ******c --- array allocation ******c ****** kk = 0 ****** kkmax = max(kkin,kkout) ****** call plot_alloc ******c ******c ****** allocate( incoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) ******c ******c call zaiost ****** flnm_d='regional.depth' ****** call getdepth(flnm_d) ******c ****** call bigrid(depths) ******c ******c do j=1,jj ******c do i=1,ii ******c p(i,j,1)=0. ******c enddo ******c enddo ******c ******c --- read the archive file, from "*.[ab]". ******c --- Use hcom_mean MPI Version of Sept. 2009 ******c kk = kkin ****** call getdat(flnm_i,time3,iweight,mntype,lsteric,icegln,trcout, ****** & iexpt,iversn,yrflag,kkin) ! hycom input ****** time = time3(3) ****** artype=mntype+1 ****** ****** if (artype.eq.3) then ****** if(mnproc.eq.1)then ****** write(lp,*) ****** write(lp,*) 'error - cannot remap std.dev. archive' ****** write(lp,*) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call xcstop('ncoda_archv_inc') ****** stop ****** endif C C READ THE INPUT PROFILE, START THE OUTPUT HEADER C OPEN(UNIT=11, FILE=flnm_i, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(flnm_i) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF OPEN(UNIT=21, FILE=flnm_o, FORM='FORMATTED', STATUS='NEW', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(flnm_o) WRITE(6,*) 'ios = ',ios CALL EXIT(5) ENDIF DO K= 1,3 READ( 11,'(a)') CLINE WRITE(21,'(a)') TRIM(CLINE) ENDDO READ( 11,'(a)') CLINE READ(CLINE(2:),*) time, srfht(1,1), & skip,skip, !surflx,dpbl & dpmixl(i,j), & skip,skip,skip,skip,skip, !Xmix & ubaro(1,1),vbaro(1,1),steric(1,1) srfht(1,1) = srfht(1,1)*g/100.0 steric(1,1) = steric(1,1)*g/100.0 dpmixl(1,1) = dpmixl(1,1)*onem ubaro(1,1) = ubaro(1,1)*0.01 vbaro(1,1) = vbaro(1,1)*0.01 DO K= 1,99 READ( 11,'(a)') CLINE IF (CLINE(1:5).EQ.'# k ') then EXIT ENDIF ENDDO p(1,1,1) = 0.0 DO K= 1,KK READ(11,'(a)',IOSTAT=IOS) CLINE IF (IOS.NE.0) THEN WRITE(6,*) 'Error: inconsistent input profile' CALL EXIT(6) ENDIF READ(CLINE,*) KI, u(1,1,k), v(1,1,k), & temp(1,1,k),saln(1,1,k),th3d(1,1,k), & dp(1,1,k) u(1,1,k) = u(1,1,k)*0.01 v(1,1,k) = v(1,1,k)*0.01 th3d(1,1,k) = th3d(1,1,k)-thbase dp(1,1,k) = dp(1,1,k)*onem p(1,1,k+1) = p(1,1,k) + dp(1,1,k) ENDDO CLOSE(11) depths(1,1) = p(1,1,kk+1) ip(1,1) = 1 icegln = .false. !ignore ice sigver = 6 !hardwired 17-term sig2 c c --- maximum absolute T increment (frcrng * vertical range) c ****** allocate( tincmx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) c tincmx(:,:) = 999.0 if (frcrng.gt.0.0) then do j= 1,jj do i= 1,ii if (ip(i,j).eq.1) then tincmx(i,j) = frcrng*(maxval(temp(i,j,:)) - & minval(temp(i,j,:)) ) endif enddo !i enddo !j ****** call xctilr(tincmx,1,1,nbdy,nbdy, halo_ps) endif !frcrng c c --- full-region NCODA c ******c --- 'flnm_t' = name of ncoda temperature increment file ******c --- 'flnm_s' = name of ncoda salinity increment file ******c --- 'flnm_u' = name of ncoda u-velocity increment file, or "NONE" ******c --- 'flnm_v' = name of ncoda v-velocity increment file, or "NONE" ******c --- 'flnm_p' = name of ncoda density displacement file, or "NONE" c --- 'flnm_a' = name of ncoda u,v,t,s,p increment file c --- 'flnm_m' = name of ncoda region mask file, or "NONE" c --- 'flnm_c' = name of ncoda ice concentration file, or "NONE" c --- 'flnm_z' = name of ncoda cell interface depths (text file), or "NONE" c --- 'kncoda' = number of input ncoda levels c --- 'kncodz' = number of additional zero increment ncoda levels c if(mnproc.eq.1)then write(lp,*) write(lp,*) 'full-region NCODA' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block read (uoff+99,'(a)') flnm_t flnm_s = flnm_t flnm_p = flnm_t flnm_u = 'NONE' flnm_v = 'NONE' if(mnproc.eq.1)then write (lp,'(2a)') 'Tncoda file: ',trim(flnm_t) endif ! end if(mnproc.eq.1) block ****** read (uoff+99,'(a)') flnm_s if(mnproc.eq.1)then write (lp,'(2a)') 'Sncoda file: ',trim(flnm_s) endif ! end if(mnproc.eq.1) block ****** read (uoff+99,'(a)') flnm_u if(mnproc.eq.1)then write (lp,'(2a)') 'Uncoda file: ',trim(flnm_u) endif ! end if(mnproc.eq.1) block ****** read (uoff+99,'(a)') flnm_v if(mnproc.eq.1)then write (lp,'(2a)') 'Vncoda file: ',trim(flnm_v) endif ! end if(mnproc.eq.1) block ****** read (uoff+99,'(a)') flnm_p if(mnproc.eq.1)then write (lp,'(2a)') 'Pncoda file: ',trim(flnm_p) endif ! end if(mnproc.eq.1) block lprsout = .false. ****** l = len_trim(flnm_p) ****** lprsout = flnm_p(l-3:l).eq.'_OUT' !undocumented debug option ****** if (lprsout) then ****** flnm_p = flnm_p(1:l-4) ****** if(mnproc.eq.1)then ****** write (lp,'(2a)') 'Pncoda file: ',trim(flnm_p) ****** endif ! end if(mnproc.eq.1) block ****** endif if (flnm_p.eq."NONE") then if (nddflg.ne.0) then if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - must have nddflg==0 with this flnm_p' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif else if (nddflg.le.0) then if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - must have nddflg>0 with this flnm_p' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif endif read (uoff+99,'(a)') flnm_m if(mnproc.eq.1)then write (lp,'(2a)') 'Mncoda file: ',trim(flnm_m) endif ! end if(mnproc.eq.1) block read (uoff+99,'(a)') flnm_c if(mnproc.eq.1)then write (lp,'(2a)') 'Cncoda file: ',trim(flnm_c) endif ! end if(mnproc.eq.1) block if (flnm_c.eq."NONE") then icegln = .false. !no ice in output archive elseif (.not.icegln) then !flnm_c.ne."NONE" if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - input archive has no ice' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif read (uoff+99,'(a)') flnm_z if(mnproc.eq.1)then write (lp,'(2a)') 'Zncoda file: ',trim(flnm_z) write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call blkini(kncoda,'kncoda') call blkini(kncodz,'kncodz') if (flnm_z.eq."NONE") then if (kncoda.gt.kkin) then if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - kncoda > kdm for layer-space case' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif if (kncodz.ne.0) then if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error - kncodz must be 0 for layer-space case' write(lp,*) call flush(lp) endif ! end if(mnproc.eq.1) block call xcstop('ncoda_archv_inc') stop endif endif close(uoff+99) !finished with blkdat.input C C --- Arctic if j=jtdm is open. C --- Note that lartic can only be true on the last tile (with j0+jj=jtdm). C larctic = .false. ****** if(j0+jj.eq.jtdm)then ****** do i= 1,ii ****** larctic = larctic .or. ip(i,jj).eq.1 ****** enddo ****** endif if (larctic) then write(lp,*) write(lp,'(a)') 'region has an arctic boundary' write(lp,*) call flush(lp) endif c c --- read kncoda sets of ncoda fields and kncoda+kncodz levels. c --- extend from kncoda to kncoda+kncodz with zero increments. c if (flnm_z.ne."NONE") then if(mnproc.eq.1)then write(lp,*) 'open ',trim(flnm_z) call flush(lp) endif ! end if(mnproc.eq.1) block ****** call zhopnc(9, flnm_z, 'FORMATTED', 'OLD', 0) OPEN(UNIT=9, FILE=flnm_z, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(flnm_z) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF do k= 0,kncoda+kncodz read(9,*) zi(k) zi(k) = onem*zi(k) enddo !k close(unit=9) if(mnproc.eq.1)then write(lp,*) 'close ',trim(flnm_z) call flush(lp) endif ! end if(mnproc.eq.1) block zz(1) = 0.0 do k= 2,kncoda+kncodz zz(k) = 0.5*(zi(k-1)+zi(k)) enddo !k endif !flnm_z c ****** allocate(tncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kncoda+kncodz), ****** & sncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kncoda+kncodz), ****** & pncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kncoda+kncodz), ****** & cncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), ****** & mncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) c if(mnproc.eq.1)then write(lp,*) 'open ',trim(flnm_t) call flush(lp) endif ! end if(mnproc.eq.1) block OPEN(UNIT=9, FILE=flnm_t, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(flnm_t) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF do k= 1,kncoda READ(9,*) uncoda(1,1,k),vncoda(1,1,k), & tncoda(1,1,k),sncoda(1,1,k),pncoda(1,1,k) enddo close(unit=9) ****** call zaiopf(flnm_t, 'old', 9) ****** call zaiord3(tncoda,kncoda,ip,.false., hminn,hmaxn, 9) ****** call zaiocl(9) ****** call xctilr(tncoda(1-nbdy,1-nbdy,1),1,kncoda,nbdy,nbdy, halo_ps) if(mnproc.eq.1)then write(lp,*) 'close ',trim(flnm_t) call flush(lp) ****** write(lp,*) 'open ',trim(flnm_s) ****** call flush(lp) endif ! end if(mnproc.eq.1) block ****** call zaiopf(flnm_s, 'old', 9) ****** call zaiord3(sncoda,kncoda,ip,.false., hminn,hmaxn, 9) ****** call zaiocl(9) ****** call xctilr(sncoda(1-nbdy,1-nbdy,1),1,kncoda,nbdy,nbdy, halo_ps) ****** if(mnproc.eq.1)then ****** write(lp,*) 'close ',trim(flnm_s) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block if (kncodz.gt.0) then tncoda(:,:,kncoda+1:kncoda+kncodz) = 0.0 sncoda(:,:,kncoda+1:kncoda+kncodz) = 0.0 endif ******c ****** if (flnm_u.ne."NONE") then ****** allocate( uncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kncoda), ****** & vncoda(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kncoda) ) ****** allocate( uinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kkout), ****** & vinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kkout) ) ******c ****** if(mnproc.eq.1)then ****** write(lp,*) 'open ',trim(flnm_u) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call zaiopf(flnm_u, 'old', 9) ****** call zaiord3(uncoda,kncoda,ip,.false., hminn,hmaxn, 9) ****** call zaiocl(9) ****** call xctilr(uncoda(1-nbdy,1-nbdy,1),1,kncoda,nbdy,nbdy, halo_uv) ****** if(mnproc.eq.1)then ****** write(lp,*) 'close ',trim(flnm_u) ****** call flush(lp) ******c ****** write(lp,*) 'open ',trim(flnm_v) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call zaiopf(flnm_v, 'old', 9) ****** call zaiord3(vncoda,kncoda,ip,.false., hminn,hmaxn, 9) ****** call zaiocl(9) ****** call xctilr(vncoda(1-nbdy,1-nbdy,1),1,kncoda,nbdy,nbdy, halo_vv) ****** if(mnproc.eq.1)then ****** write(lp,*) 'close ',trim(flnm_v) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block if (kncodz.gt.0) then uncoda(:,:,kncoda+1:kncoda+kncodz) = 0.0 vncoda(:,:,kncoda+1:kncoda+kncodz) = 0.0 endif ****** endif !u&vncoda ******c ****** if (flnm_p.ne."NONE") then ****** if(mnproc.eq.1)then ****** write(lp,*) 'open ',trim(flnm_p) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call zaiopf(flnm_p, 'old', 9) ****** call zaiord3(pncoda,kncoda,ip,.false., hminn,hmaxn, 9) ****** call zaiocl(9) ****** call xctilr(pncoda(1-nbdy,1-nbdy,1),1,kncoda,nbdy,nbdy, halo_ps) ****** if(mnproc.eq.1)then ****** write(lp,*) 'close ',trim(flnm_p) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block if (kncodz.gt.0) then pncoda(:,:,kncoda+1:kncoda+kncodz) = 0.0 endif ****** endif !pncoda c c --- increase kncoda to allow for kncodz c kncoda = kncoda + kncodz c ****** if (flnm_c.ne."NONE") then ****** if(mnproc.eq.1)then ****** write(lp,*) 'open ',trim(flnm_c) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call zaiopf(flnm_c, 'old', 9) ****** call zaiord(cncoda,ip,.false., hmina,hmaxa, 9) ****** call zaiocl(9) ****** call xctilr(cncoda(1-nbdy,1-nbdy),1,1,nbdy,nbdy, halo_ps) ****** if(mnproc.eq.1)then ****** write(lp,*) 'close ',trim(flnm_c) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** endif !cncoda ******c ****** if (flnm_m.ne."NONE") then ****** if(mnproc.eq.1)then ****** write(lp,*) 'open ',trim(flnm_m) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** call zaiopf(flnm_m, 'old', 9) ****** call zaiord(mncoda,ip,.false., hmina,hmaxa, 9) ****** call zaiocl(9) ****** call xctilr(mncoda(1-nbdy,1-nbdy),1,1,nbdy,nbdy, halo_ps) ****** if(mnproc.eq.1)then ****** write(lp,*) 'close ',trim(flnm_m) ****** call flush(lp) ****** endif ! end if(mnproc.eq.1) block ****** do j= 1-nbdy,jj+nbdy ****** do i= 1-nbdy,ii+nbdy ****** if (mncoda(i,j).lt.2.0**99) then ****** incoda(i,j) = 0 !ncoda sea point, default is all zero increments ****** do k= 1,kncoda ****** if (min(tncoda(i,j,k), ****** & sncoda(i,j,k) ).le.-900.0) then ****** exit !ncoda sea point, all zero increments ****** elseif (tncoda(i,j,k).ne.0.0 .or. ****** & sncoda(i,j,k).ne.0.0 ) then ****** incoda(i,j) = 1 !ncoda sea point, non-zero increments ****** exit ****** endif ****** enddo !k, exit ******c ******c --- inforce maximum absolute T and S increment ******c ****** if (incoda(i,j).eq.1) then ****** do k= 1,kncoda ****** tncoda(i,j,k) = max( -tincmx(i,j), ****** & min( tincmx(i,j), tncoda(i,j,k) ) ) ****** enddo !k ****** endif !incoda ****** else ****** incoda(i,j) = 0 !ncoda land point ****** endif ****** enddo ****** enddo ****** else do j= 1-nbdy,jj+nbdy do i= 1-nbdy,ii+nbdy if (ip(i,j).eq.1) then incoda(i,j) = 0 !ncoda sea point, default is all zero increments do k= 1,kncoda if (min(tncoda(i,j,k), & sncoda(i,j,k) ).le.-900.0) then exit !ncoda sea point, all zero increments elseif (tncoda(i,j,k).ne.0.0 .or. & sncoda(i,j,k).ne.0.0 ) then incoda(i,j) = 1 !ncoda sea point, non-zero increments exit endif enddo !k, exit c c --- inforce maximum absolute T and S increment c if (incoda(i,j).eq.1) then do k= 1,kncoda tncoda(i,j,k) = max( -tincmx(i,j), & min( tincmx(i,j), tncoda(i,j,k) ) ) enddo !k endif !incoda else incoda(i,j) = 0 !ncoda land point endif enddo enddo ****** endif !mncoda:else ******c ******c --- Ice Concentration and Thickness. ******c ****** if (flnm_c.ne."NONE") then ****** do j= 1,jj ****** do i= 1,ii ****** if (incoda(i,j).eq.1) then ****** if (cncoda(i,j).gt.2.0**99) then ****** write(lp,*) ****** write(lp,*) 'error - cncoda data void over sea' ****** write(lp,*) 'mnproc,i,j = ',mnproc,i+i0,j+j0 ****** write(lp,*) ****** call flush(lp) ****** call xchalt('ncoda_archv_inc') ****** stop ****** else ******c --- percent to fraction covered, 0 to 1 ****** cncoda(i,j) = max(0.0, min(1.0, 0.01*cncoda(i,j))) ****** endif ****** if (abs(covice(i,j)-cncoda(i,j)).gt.0.03) then ****** if (cncoda(i, j) .le.0.00001) then ****** thkice(i,j) = 0.0 ****** covice(i,j) = 0.0 ****** elseif (covice(i,j).le.0.00001) then !hice=hicemn ****** thkice(i,j) = cncoda(i,j)*hicemn !average ice thickness ****** covice(i,j) = cncoda(i,j) ****** else !hice=thkice/covice ****** thkice(i,j) = cncoda(i,j)*thkice(i,j)/ ****** & covice(i,j) ****** covice(i,j) = cncoda(i,j) ******c --- covice will be reset to 1.0 by HYCOM if thkice/covice > hicemn ****** endif ****** endif !change ice coverage ****** endif !incoda ****** enddo !i ****** enddo !j ****** if (larctic) then !arctic patch ****** do i= 1,ii ****** ih = idm-mod(i-1,idm) ****** thkice(i,jj) = thkice(ih,jj-1) ****** covice(i,jj) = covice(ih,jj-1) ****** enddo !i ****** endif !larctic ****** endif !cncoda c c --- mixed-layer depth in layer space (1 < mldlay < kk+1). c --- must be at least as deep as the fixed vertical coordinate layers. c --- only calculate this once. c ****** if ((flnm_p.ne."NONE" .or. isoflg.ne.0) .and. ****** & .not. allocated(mldlay) ) then ****** allocate( pij(kkin+1), mldlay(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), ****** & work(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ) ****** allocate( theta3(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kkin) ) c c --- initialize theta3, from a file if necessary ****** if (vsigma) then ****** call zaiopf('iso.sigma.a', 'old', 922) ****** do k=1,kkin ****** call zaiord(work,ip,.false., hmina,hmaxa, 922) ****** call xctilr(work(1-nbdy,1-nbdy),1,1,nbdy,nbdy, halo_ps) ****** theta3(:,:,k) = work(:,:) - thbase ****** enddo ****** call zaiocl(922) ****** else do k=1,kkin theta3(:,:,k) = sigma(k) - thbase enddo ****** endif c do j= 1,jj do i= 1,ii if (ip(i,j).eq.1) then mldlay(i,j) = kkin !usually modified below mldij = max( dpmixl(i,j), isotop ) dp0cum(1) = 0.0 pij( 1) = 0.0 do k= 1,kkin c --- q is dp0k(k) when in surface fixed coordinates c --- q is dp00i when much deeper than surface fixed coordinates if (dp0k(k).le.dp00i .or. k.eq.1) then q = dp0k(k) else q = max( dp00i, & dp0k(k) * dp0k(k)/ & max( dp0k(k), & pij(k)-dp0cum(k) ) ) endif dp0ij=min(q,max(ds0k(k),dssk(k)*depths(i,j))) dp0cum(k+1) = dp0cum(k) + dp0ij pij( k+1) = pij( k) + dp(i,j,k) if (dp(i,j,k).lt.1.5*dp0ij) then mldij = max( mldij, pij(k+1) ) !fixed coordinates elseif (abs(th3d( i,j,k)- & theta3(i,j,k) ).gt.deniso) then mldij = max( mldij, pij(k+1) ) !non-isopycnal coordinates endif if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,th3d-target = ', & k,th3d(i,j,k), & abs(th3d(i,j,k)-theta3(i,j,k)) write(6,*) 'pij,mldij=',pij(k+1)*qonem,mldij*qonem call flush(lp) endif if (mldij.lt.pij(k+1)) then if (isotop.ge.pij(k)) then !always fixed zone mldlay(i,j) = k + 1 else mldlay(i,j) = k + (mldij - pij(k)) / dp(i,j,k) endif exit endif enddo !k if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(6,*) 'mld = ',dpmixl(i,j)*qonem !meters write(6,*) 'mldlay = ',mldlay(i,j) !layers endif endif !ip enddo !i enddo !j ******c ******c --- smooth three times ******c ****** call xctilr(mldlay,1,1,nbdy,nbdy, halo_ps) ****** call psmoo(mldlay, work) ****** call xctilr(mldlay,1,1,nbdy,nbdy, halo_ps) ****** call psmoo(mldlay, work) ****** call xctilr(mldlay,1,1,nbdy,nbdy, halo_ps) ****** call psmoo(mldlay, work) ****** call xctilr(mldlay,1,1,nbdy,nbdy, halo_ps) ****** if (min(itest,jtest).gt.0) then ****** if(i0.lt.itest.and.i0+ii.ge.itest.and. ****** & j0.lt.jtest.and.j0+jj.ge.jtest)then ****** write(lp,*)'mldlaysm = ',mldlay(itest-i0,jtest-j0) !layers ****** ****** endif ! itest,jtest lie on current block ****** endif ******c ****** deallocate( pij, work) !keep mldlay and theta3 ****** endif !mldlay, if needed ******c ******c --- are we saving prsout? ******c ****** if (lprsout) then ****** allocate( prsout(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kncoda) ) ****** prsout(:,:,:) = 0.0 !default is no displacement ****** endif c c --- form exisiting and target interface depths. c ****** if (.not. allocated(pout)) then ****** allocate( pout(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kkout+1) ) ****** endif c do j= 1,jj do i= 1,ii if (ip(i,j).eq.1) then * write(lp,'(a,2i4)') 'p - i,j = ',i,j * call flush(lp) p(i,j,1) = 0.0 do k= 1,kkin-1 p(i,j,k+1) = min(p(i,j,k) + dp(i,j,k), & depths(i,j)) enddo !k p(i,j,kkin+1) = depths(i,j) c c --- default is that interfaces remain the same c do k= 1,kkout+1 !kkout=kkin pout(i,j,k) = p(i,j,k) enddo !k endif !ip enddo !i enddo !j ****** call xctilr(p(1-nbdy,1-nbdy,1),1,kkout+1,nbdy,nbdy, halo_ps) c if (flnm_p.ne."NONE") then do j= 1,jj do i= 1,ii if (incoda(i,j).eq.1) then c c --- initialize pncoda c if (flnm_z.eq."NONE") then do k= 0,kncoda zi(k) = p(i,j,k+1) enddo !k zz(1) = 0.0 do k= 2,kncoda zz(k) = 0.5*(zi(k-1)+zi(k)) enddo !k endif !no flnm_z lpncz = .true. do k= 1,kncoda if (pncoda(i,j,k).ne.0.0) then lpncz = .false. exit endif enddo !k if (lpncz) then c --- input pncoda is zero, use zero pnc(1:kncoda) = 0.0 elseif (nddflg.eq.1) then c --- use the input pncoda pnc(1) = max(0.0, pncoda(i,j,1)*1.0198*onem ) do k= 2,kncoda pnimin = (zz(k-1) -zz(k)) ! negative pnmax = (zz(kncoda)-zz(k)) !non-negative c surface moves from zz(k) to zz(k)+pncoda(i,j,k), which c must be below zz(k-1)+pncoda(i,j,k-1) but above zz(kncoda) pnc(k) = pncoda(i,j,k)*1.0198*onem pnc(k) = max( pnc(k), & pnc(k-1)+pnimin ) pnc(k) = min( pnc(k), pnmax ) enddo !k else !nddflg.ge.2 c --- calculate a new "pncoda" from the change in density profile c --- nddflg==2; use the smaller of the input and new pncoda c --- nddflg==3; use the new pncoda if (flnm_z.eq."NONE") then !HYCOM layer structure do k= 1,kncoda !same as 1,kkout to(k) = temp(i,j,k) so(k) = saln(i,j,k) enddo !k else !remap the original layers to NCODA layers p1(0) = 0.0 do k= 1,kkout p1(k) = p(i,j,k+1) t1(k) = temp(i,j,k) s1(k) = saln(i,j,k) enddo do k= 1,kncoda to(k) = -999.0 !data void marker so(k) = -999.0 !data void marker enddo * if (i+i0.eq.itest .and. j+j0.eq.jtest) then * call remap_plm_1_debug(t1,p1,kkout, * & to,zi,kncoda, i+i0,j+j0) * call remap_plm_1_debug(s1,p1,kkout, * & so,zi,kncoda, i+i0,j+j0) * else call remap_plm_1(t1,p1,kkout, & to,zi,kncoda, i+i0,j+j0) call remap_plm_1(s1,p1,kkout, & so,zi,kncoda, i+i0,j+j0) * endif endif !flnm_z c --- calculate the orignal and new density profile do k= 1,kncoda ti(k) = tncoda(i,j,k) si(k) = sncoda(i,j,k) if (min(ti(k),si(k)).le.-900.0 .or. & zz(k).gt.depths(i,j) ) then ti(k) = ti(k-1) !inherit from above si(k) = si(k-1) endif if (min(to(k),so(k)).le.-900.0 .or. & zz(k).gt.depths(i,j) ) then to(k) = to(k-1) !inherit from above so(k) = so(k-1) endif tz(k) = to(k) + ti(k) sz(k) = so(k) + si(k) ro(k) = SIG_V(to(k),so(k),sigver) rz(k) = SIG_V(tz(k),sz(k),sigver) rz(k) = max( rz(k), rz(max(k-1,1)) ) !stable profile * if (i+i0.eq.itest .and. j+j0.eq.jtest) then * write(lp,*) 'k,ro,rz = ',k,ro(k),rz(k) * endif enddo !k * if (i+i0.eq.itest .and. j+j0.eq.jtest) then * write(lp,*) 'k,depth = ', * & k,min(zz(kncoda), * & depths(i,j))*qonem * call flush(lp) * endif * if (i+i0.eq.itest .and. j+j0.eq.jtest) then * call layprs_debug(ro,zz,rz,zz,pnc,kncoda, * & min(zz(kncoda),depths(i,j)) ) * else c --- calculate a new "pncoda" call layprs(ro,zz,rz,zz,pnc,kncoda, & min(zz(kncoda),depths(i,j)) ) * endif c --- smooth it twice pns(1:kncoda) = pnc(1:kncoda) do k= 2,kncoda-1 pnc(k) = 0.5*pns(k) + 0.25*(pns(k-1)+pns(k+1)) enddo !k pns(1:kncoda) = pnc(1:kncoda) do k= 2,kncoda-1 pnc(k) = 0.5*pns(k) + 0.25*(pns(k-1)+pns(k+1)) enddo !k if (nddflg.eq.2) then c --- select between layprs pnc and pncoda pnck = max(0.0, pncoda(i,j,1)*1.0198*onem ) if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,pnck = ',1,pnck*qonem write(lp,*) 'k,pnc = ',1,pnc(1)*qonem call flush(lp) endif if (abs(pnck).lt.abs(pnc(1))) then pnc(1) = pnck if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,pncNEW = ',1,pnc(1)*qonem call flush(lp) endif endif do k= 2,kncoda pnimin = (zz(k-1) -zz(k)) ! negative pnmax = (zz(kncoda)-zz(k)) !non-negative c surface moves from zz(k) to zz(k)+pncoda(i,j,k), which c must be below zz(k-1)+pncoda(i,j,k-1) but above zz(kncoda) pnck = pncoda(i,j,k)*1.0198*onem pnck = max( pnck, & pnc(k-1)+pnimin ) pnck = min( pnck, pnmax ) if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,pncoda = ', & k,pncoda(i,j,k)*1.0198*onem*qonem write(lp,*) 'k,pnmin = ', & k,(pnc(k-1)+pnimin)*qonem write(lp,*) 'k,pnmax = ',k,pnmax*qonem write(lp,*) 'k,pnck = ',k,pnck*qonem write(lp,*) 'k,pnc = ',k,pnc(k)*qonem call flush(lp) endif if (abs(pnck).lt.abs(pnc(k))) then pnc(k) = pnck if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,pncNEW = ',k,pnc(k)*qonem call flush(lp) endif endif enddo !k elseif (i+i0.eq.itest .and. j+j0.eq.jtest) then !nddflg==3 do k= 1,kncoda write(lp,*) 'k,pnck = ', & k,pncoda(i,j,k)*1.0198*onem write(lp,*) 'k,pncNEW = ', & k,pnc(k)*qonem call flush(lp) enddo !k endif !nddflg==2:3 endif !nddflg if (lprsout) then do k= 1,kncoda prsout(i,j,k) = pnc(k)/(1.0198*onem) enddo !k endif !lprsout c c --- move interfaces based on pncoda (pnc) c kmld = int(mldlay(i,j)) dp0cum(1)=0.0 pout(i,j,1) = 0.0 do k= 2,kkout !kkout=kkin c --- q is dp0k(k-1) when in surface fixed coordinates c --- q is dp00i when much deeper than surface fixed coordinates if (dp0k(k-1).le.dp00i) then q = dp0k(k-1) else q = max( dp00i, & dp0k(k-1) * dp0k(k-1)/ & max( dp0k( k-1), & p(i,j,k-1)-dp0cum(k-1) ) ) endif dp0ij=min(q,max(ds0k(k-1),dssk(k-1)*depths(i,j))) dp0cum(k)=dp0cum(k-1)+dp0ij pk = p(i,j,k) !no motion if (k-1.ge.kmld) then c --- below the mixed layer and isotop do l=2,kncoda if (zz(l).gt.pk) then qk = (zz(l)-pk)/(zz(l)-zz(l-1)) pnk = (1.0-qk)*pnc(l) + & qk *pnc(l-1) * if (k-1.eq.kmld) then * --- scale by fraction of layer deeper than the mixed-layer * pk = pk + pnk*(kmld+1.0-mldlay(i,j)) * else pk = pk + pnk * endif if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,l,zz = ', & k,l,zz(l-1)*qonem,zz(l)*qonem write(lp,*) 'k,l,qk = ', & k,l,qk write(lp,*) 'k,l,pn = ', & k,l,pnc(l-1)*qonem, & pnc(l) *qonem write(lp,*) 'k,l,pk = ', & k,l,p(i,j,k)*qonem,pk*qonem call flush(lp) endif exit endif enddo !l elseif (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,th3d,target = ', & k-1,th3d(i,j,k-1),thk(k-1) !in the mixed layer call flush(lp) endif !below mixed layer pout(i,j,k) = min( max( pk, & pout(i,j,k-1)+dp0ij ), & depths(i,j) ) c c --- rebalance layer when thick-thin-thick or thin-thick-thin c --- but only if ncoda is making the imbalance worse c if (k.gt.2 .and. & th3d( i,j,k-2).lt. & theta3(i,j,k-1)-epsil .and. & th3d( i,j,k) .gt. & theta3(i,j,k-1)+epsil ) then dptop = p(i,j,k-1) - p(i,j,k-2) dpmid = p(i,j,k) - p(i,j,k-1) dpbot = p(i,j,k+1) - p(i,j,k) if (dp0ij.lt.thnthk*min(dptop,dpbot) .and. & dpmid.lt.thnthk*min(dptop,dpbot) ) then dpnew = pout(i,j,k) - pout(i,j,k-1) if (dpnew.lt.dpmid) then c c --- thicken layer k-1 by taking fluid from above and below c --- in proportion so that the net is at k-1''s target density c if (i+i0.eq.itest .and. j+j0.eq.jtest) then write(lp,*) 'k,l,po = ', & k,-1, p(i,j,k)*qonem, & pout(i,j,k)*qonem call flush(lp) endif !test q = (theta3(i,j,k-1)-th3d(i,j,k-2))/ & (th3d( i,j,k) -th3d(i,j,k-2)) ! 0 depth. c ro.z at z=zi(k)+dz(k) should be ri(k), but this will only be c the case where ri is locally non-decreasing. c c 5) Alan J. Wallcraft, Naval Research Laboratory, May 2012. c* c********** c real thin parameter (thin=1.e-6) ! minimum layer thickness c integer k,ko,kom1 real q,z,zold c kom1 = 1 zold = zo(1) do k= 1,kk c --- find the location of ri(k) in ro(:). if (zi(k).ge.depth) then dz(k:kk) = 0.0 exit elseif (ro(kom1).ge.ri(k)) then dz(k) = zold - zi(k) !zold and kom1 unchanged else do ko= kom1+1,kk if (ro(ko).ge.ri(k)) then !also ro(ko-1).lt.ri(k) q = (ro(ko) - ri(k))/max(ro(ko) - ro(ko-1), thin) z = q*zo(ko-1) + (1.0-q)*zo(ko) z = min( z, depth ) zold = max( z, zold ) dz(k) = zold - zi(k) kom1 = ko-1 exit endif if (ko.eq.kk) then z = zo(ko) z = min( z, depth ) zold = max( z, zold ) dz(k) = zold - zi(k) kom1 = ko-1 exit endif enddo !ko endif !ro.kom1:else enddo !k return end subroutine layprs subroutine layprs_debug(ri,zi,ro,zo,dz,kk, depth) implicit none c integer kk real ri(kk),zi(kk), & ro(kk),zo(kk),dz(kk),depth c c********** c* c 1) calculate the isopycnal displacement between two density profiles. c c 2) input arguments: c ri - 1st density profile values c zi - 1st density profile depths c ro - 2nd density profile values c zo - 2nd density profile depths c kk - number of levels c depth - maximum depth, can be less than zi(kk) and zo(kk) c c 3) output arguments: c dz - displacement between the two density profiles c c 4) assumes that ro is a non-decreasing profile. c zi(:)+dz(:) will be non-decreasing, between zo(1) and depth. c dz will be zero where zi > depth. c ro.z at z=zi(k)+dz(k) should be ri(k), but this will only be c the case where ri is locally non-decreasing. c c 5) Alan J. Wallcraft, Naval Research Laboratory, May 2012. c* c********** c real thin parameter (thin=1.e-6) ! minimum layer thickness c integer k,ko,kom1 real q,z,zold c kom1 = 1 zold = zo(1) WRITE(6,*) 'kom1,zold,ro = ',kom1,zold,ro(kom1) do k= 1,kk c --- find the location of ri(k) in ro(:). if (zi(k).ge.depth) then dz(k:kk) = 0.0 WRITE(6,*) 'k,zi,depth = ',k,zi(k),depth exit elseif (ro(kom1).ge.ri(k)) then dz(k) = zold - zi(k) !zold and kom1 unchanged WRITE(6,*) 'K,DZ,RI = ',k,dz(k),ri(k) else WRITE(6,*) 'k,ri = ',k,ri(k) do ko= kom1+1,kk if (ro(ko).ge.ri(k)) then !also ro(ko-1).lt.ri(k) q = (ro(ko) - ri(k))/max(ro(ko) - ro(ko-1), thin) z = q*zo(ko-1) + (1.0-q)*zo(ko) z = min( z, depth ) zold = max( z, zold ) dz(k) = zold - zi(k) kom1 = ko-1 WRITE(6,*) 'ko,ro = ',ko,ro(ko) WRITE(6,*) 'k,dz,ri = ',k,dz(k),ri(k) WRITE(6,*) 'kom1,zold,ro = ',kom1,zold,ro(kom1) exit endif if (ko.eq.kk) then z = zo(ko) z = min( z, depth ) zold = max( z, zold ) dz(k) = zold - zi(k) kom1 = ko-1 WRITE(6,*) 'KO,RO = ',ko,ro(ko) WRITE(6,*) 'K,DZ,RI = ',k,dz(k),ri(k) WRITE(6,*) 'KOM1,ZOLD,RO = ',kom1,zold,ro(kom1) exit endif enddo !ko endif !ro.kom1:else enddo !k return end subroutine layprs_debug subroutine remap_plm_1(t, p, kk, & tz,pz,kz, i,j) ****** use mod_xc ! Hycom mpi-2 interface implicit none c integer kk,kz,i,j real t( kk),p( kk+1), & tz(kz),pz(kz+1) c c********** c* c 1) remap from one set of vertical cells to another. c method: piecewise linear across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c t - scalar field in p-layer space c p - layer interface depths (non-negative m) c p( 1) is the surface c p(kk+1) is the bathymetry c kk - dimension of t (number of input layers) c pz - target interface depths (non-negative m) c pz(k+1) >= pz(k) c kz - dimension of tz (number of output layers) c i,j - horizontal grid location, for debugging only c c 3) output arguments: c tz - scalar field in pz-layer space c c 4) must have: c 0 = p(1) <= pz(l) <= pz(l+1) c 0 <= pz(k) <= pz(k+1) c output layer spaning p(kk+1) uses input tz unchanged c output layers below p(kk+1) return input tz unchanged c c 5) Tim Campbell, Mississippi State University, October 2002. c* c********** c real,parameter :: thin=9806.0e-6 !minimum layer thickness c integer k,l,lf real q,qc,qzbt,zb,zc,zt,tzk real ts(kk),pt(kk+1) c c --- compute PLM slopes for input layers do k=1,kk pt(k)=max(p(k+1)-p(k),thin) enddo call plm1(pt, t, ts, kk) c --- compute output layer averages lf=1 zb=pz(1) do k= 1,kz zt = zb zb = pz(k+1) * WRITE(6,*) 'k,zt,zb = ',k,zt/9806.0,zb/9806.0 if (zb-zt.lt.thin) then c c --- thin layer, values taken from layer above c tz(k) = tz(k-1) elseif (zb.gt.p(kk+1)) then c c --- layer spanning or below the "bottom", input is unchanged c * tz(k) = tz(k) else c c form layer averages. c if (p(lf).gt.zt) then WRITE(6,*) 'bad lf = ',lf,p(lf),zt,p(lf)-zt call remap_plm_1_debug(t, p, kk, & tz,pz,kz, i,j) call xchalt('ncoda_archv_inc') stop endif qzbt = 1.0/(zb-zt) !zb-zt > thin tzk = 0.0 do l= lf,kk if (p(l).ge.zb) then c WRITE(6,*) 'l,lf= ',l,lf,l-1 lf = l-1 exit elseif (p(l).ge.zt .and. p(l+1).le.zb) then c c the input layer is completely inside the output layer c q = (p(l+1)-p(l))*qzbt tzk = tzk + q*t(l) * WRITE(6,'(a,i5,3f12.7)') 'L,q,t = ', * & l,q,t(l),tzk else c c the input layer is partially inside the output layer c average of linear profile is its value at the layer center c q = (min(p(l+1),zb)-max(p(l),zt))*qzbt if (q.gt.0.0) then zc = 0.5*(min(p(l+1),zb)+max(p(l),zt)) qc = (zc-p(l))/pt(l) - 0.5 tzk = tzk + q*(t(l) + qc*ts(l)) * WRITE(6,'(a,i5,3f12.7)') 'L,q,t* = ', * & l,q,(t(l)+qc*ts(l)),tzk endif endif enddo !l tz(k) = tzk endif enddo !k return end subroutine remap_plm_1 subroutine remap_plm_1_debug(t, p, kk, & tz,pz,kz, i,j) ****** use mod_xc ! Hycom mpi-2 interface implicit none c integer kk,kz,i,j real t( kk),p( kk+1), & tz(kz),pz(kz+1) c c********** c* c 1) remap from one set of vertical cells to another. c method: piecewise linear across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c t - scalar field in p-layer space c p - layer interface depths (non-negative m) c p( 1) is the surface c p(kk+1) is the bathymetry c kk - dimension of t (number of input layers) c pz - target interface depths (non-negative m) c pz(k+1) >= pz(k) c kz - dimension of tz (number of output layers) c i,j - horizontal grid location, for debugging only c c 3) output arguments: c tz - scalar field in pz-layer space c c 4) must have: c 0 = p(1) <= pz(l) <= pz(l+1) c 0 <= pz(k) <= pz(k+1) c output layer spaning p(kk+1) uses input tz unchanged c output layers below p(kk+1) return input tz unchanged c c 5) Tim Campbell, Mississippi State University, October 2002. c* c********** c real,parameter :: thin=9806.0e-6 !minimum layer thickness c logical lfatal integer k,l,lf real q,qc,qzbt,zb,zc,zt,tzk real ts(kk),pt(kk+1) c c --- check pz lfatal = .false. do k= 2,kz if (pz(k).lt.pz(k-1)) then WRITE(6,*) 'bad pz = ',k-1,pz(k-1)/9806.,k,pz(k)/9806. lfatal = .true. elseif (lfatal) then WRITE(6,*) ' pz = ',k-1,pz(k-1)/9806.,k,pz(k)/9806. endif enddo !k if (lfatal) then WRITE(6,*) 'fatal error at global',i,j call xchalt('ncoda_archv_inc') stop endif c c --- compute PLM slopes for input layers do k=1,kk pt(k)=max(p(k+1)-p(k),thin) enddo call plm1(pt, t, ts, kk) c --- compute output layer averages lf=1 zb=pz(1) do k= 1,kz zt = zb zb = pz(k+1) WRITE(6,*) 'k,zt,zb = ',k,zt/9806.0,zb/9806.0 if (zb-zt.lt.thin) then c c --- thin layer, values taken from layer above c tz(k) = tz(k-1) WRITE(6,*) 'k,THIN = ',k,tz(k) elseif (zb.gt.p(kk+1)) then c c --- layer spanning or below the "bottom", input is unchanged * tz(k) = tz(k) WRITE(6,*) 'k,BOTT = ',k,tz(k) else c c form layer averages. c if (p(lf).gt.zt) then WRITE(6,*) 'bad lf = ',lf,' i,j =',i,j call xchalt('ncoda_archv_inc') stop endif qzbt = 1.0/(zb-zt) !zb-zt > thin tzk = 0.0 do l= lf,kk if (p(l).ge.zb) then WRITE(6,*) 'l,lf= ',l,lf,l-1 lf = l-1 exit elseif (p(l).ge.zt .and. p(l+1).le.zb) then c c the input layer is completely inside the output layer c q = (p(l+1)-p(l))*qzbt tzk = tzk + q*t(l) WRITE(6,'(a,i5,3f12.7)') 'L,q,t = ', & l,q,t(l),tzk else c c the input layer is partially inside the output layer c average of linear profile is its value at the layer center c q = (min(p(l+1),zb)-max(p(l),zt))*qzbt if (q.gt.0.0) then zc = 0.5*(min(p(l+1),zb)+max(p(l),zt)) qc = (zc-p(l))/pt(l) - 0.5 tzk = tzk + q*(t(l) + qc*ts(l)) WRITE(6,'(a,i5,3f12.7)') 'L,q,t* = ', & l,q,(t(l)+qc*ts(l)),tzk endif endif enddo !l tz(k) = tzk endif enddo !k return end subroutine remap_plm_1_debug subroutine plm1(pt, t, ts,kk) implicit none c integer kk real t(kk),pt(kk),ts(kk) c c********** c* c 1) generate a monotonic PLM interpolation of a layered field c c 2) input arguments: c pt - layer interface thicknesses (non-zero) c t - scalar field in layer space c kk - dimension of a (number of layers) c c 3) output arguments: c ts - scalar field slopes for PLM interpolation c c 4) except at data voids, must have: c p( 1) == zero (surface) c p( l+1) >= p(:,:,l) c p(kk+1) == bathymetry c c 5) Tim Campbell, Mississippi State University, September 2002. c* c********** c integer l real ql(kk),qc(kk),qr(kk) c !compute grid spacing ratios for slope computations ql(1)=0.0 qc(1)=0.0 qr(1)=0.0 do l=2,kk-1 ql(l)=2.0*pt(l)/(pt(l-1)+pt(l)) qc(l)=2.0*pt(l)/(pt(l-1)+2.0*pt(l)+pt(l+1)) qr(l)=2.0*pt(l)/(pt(l)+pt(l+1)) enddo ql(kk)=0.0 qc(kk)=0.0 qr(kk)=0.0 !compute normalized layer slopes call slope(ql,qc,qr,t,ts,kk) return end subroutine plm1 subroutine slope(rl,rc,rr,a,s,n) implicit none c integer,intent(in) :: n real, intent(in) :: rl(n),rc(n),rr(n),a(n) real, intent(out) :: s(n) c c********** c* c 1) generate slopes for monotonic piecewise linear distribution c c 2) input arguments: c rl - left grid spacing ratio c rc - center grid spacing ratio c rr - right grid spacing ratio c a - scalar field zone averages c n - number of zones c c 3) output arguments: c s - zone slopes c c 4) Tim Campbell, Mississippi State University, September 2002. c* c********** c integer,parameter :: ic=2, im=1, imax=100 real,parameter :: fracmin=1e-6, dfac=0.5 c integer i,j real sl,sc,sr real dnp,dnn,dl,dr,ds,frac c c Compute zone slopes c Campbell Eq(15) -- nonuniform grid c s(1)=0.0 do j=2,n-1 sl=rl(j)*(a(j)-a(j-1)) sr=rr(j)*(a(j+1)-a(j)) if (sl*sr.gt.0.) then s(j)=sign(min(abs(sl),abs(sr)),sl) else s(j)=0.0 endif enddo s(n)=0.0 c c Minimize discontinuities between zones c Apply single pass discontinuity minimization: Campbell Eq(19) c do j=2,n-1 if(s(j).ne.0.0) then dl=-0.5*(s(j)+s(j-1))+a(j)-a(j-1) dr=-0.5*(s(j+1)+s(j))+a(j+1)-a(j) ds=sign(min(abs(dl),abs(dr)),dl) s(j)=s(j)+2.0*ds endif enddo return end subroutine slope subroutine remap_isopyc(t, r, kk, & tz,rz,kz) implicit none c integer kk,kz real t( kk),r( kk), & tz(kz),rz(kz) c c********** c* c 1) remap from one set of vertical cells to another. c method: sample at isopycnal depths c c 2) input arguments: c t - scalar field in p-layer space c r - density field in p-layer space c kk - dimension of t (number of input layers) c rz - target densities c rz(k+1) >= rz(k) c kz - dimension of tz (number of output layers) c c 3) output arguments: c tz - scalar field in density space c c 4) except at data voids, must have: c r(1) <= r(l) <= r( l+1) c a target density lighter than r(1) returns tz unchanged c a target density heavier than r(kk) returns tz unchanged c c 5) Alan Wallcraft, NRL, April 2005. c* c********** c real,parameter :: rtiny=0.001 !minimum density change c integer k,l,lf real q c lf = 1 do k= 1,kz if (r( 1).le.rz(k) .and. & r(kk).ge.rz(k) ) then if (r(lf).eq.rz(k)) then tz(k) = t(l) else do l= lf,kk if (r(l).ge.rz(k)) then q = (r(l)-rz(k))/max(r(l)-r(l-1),rtiny) tz(k) = t(l) - q*(t(l)-t(l-1)) lf = l exit !l endif !found target density enddo !l endif !r(lf)==rz(k) endif !rz(k) in range enddo !k return end subroutine remap_isopyc subroutine remap_isopyc_debug(t, r, kk, & tz,rz,kz) ****** use mod_xc ! Hycom mpi-2 interface implicit none c integer kk,kz real t( kk),r( kk), & tz(kz),rz(kz) c c********** c* c 1) remap from one set of vertical cells to another. c method: sample at isopycnal depths c c 2) input arguments: c t - scalar field in p-layer space c r - density field in p-layer space c kk - dimension of t (number of input layers) c rz - target densities c rz(k+1) >= rz(k) c kz - dimension of tz (number of output layers) c c 3) output arguments: c tz - scalar field in density space c c 4) except at data voids, must have: c r(1) <= r(l) <= r( l+1) c a target density lighter than r(1) returns tz unchanged c a target density heavier than r(kk) returns tz unchanged c c 5) Alan Wallcraft, NRL, April 2005. c* c********** c real,parameter :: rtiny=0.001 !minimum density change c integer k,l,lf real q c lf = 1 do k= 1,kz if (r( 1).le.rz(k) .and. & r(kk).ge.rz(k) ) then if (r(lf).eq.rz(k)) then tz(k) = t(lf) WRITE(6,*) 'k,lf,rz,t,tz = ',k,lf,rz(k),r(l),tz(k),t(lf) else WRITE(6,*) 'k,lf= ',k,lf do l= lf,kk if (r(l).ge.rz(k)) then q = (r(l)-rz(k))/max(r(l)-r(l-1),rtiny) tz(k) = t(l) - q*(t(l)-t(l-1)) WRITE(6,*) 'k,lm,rz,t,tz = ',k,l-1,rz(k),r(l-1), & tz(k),t(l-1),q WRITE(6,*) 'k,l, rz,t,tz = ',k,l, rz(k),r(l), & tz(k),t(l) lf = l exit !l endif !found target density enddo !l endif !r(lf)==rz(k) endif !rz(k) in range enddo !k return end subroutine remap_isopyc_debug REAL*4 FUNCTION SIG_V(TT,SS,SIGVER) IMPLICIT NONE INTEGER SIGVER REAL*4 TT,SS C C SIGVER WRAPPER FOR SIG C REAL*8 SS8,TT8 REAL*8 SIG_1,SIG_2,SIG_3,SIG_4,SIG_5,SIG_6,SIG_7,SIG_8 C TT8 = TT SS8 = SS IF (MOD(SIGVER,2).EQ.1) THEN IF (SIGVER.EQ.1) THEN SIG_V = SIG_1(TT8,SS8) ELSEIF (SIGVER.EQ.3) THEN SIG_V = SIG_3(TT8,SS8) ELSEIF (SIGVER.EQ.5) THEN SIG_V = SIG_5(TT8,SS8) ELSEIF (SIGVER.EQ.7) THEN SIG_V = SIG_7(TT8,SS8) ENDIF ELSE IF (SIGVER.EQ.2) THEN SIG_V = SIG_2(TT8,SS8) ELSEIF (SIGVER.EQ.4) THEN SIG_V = SIG_4(TT8,SS8) ELSEIF (SIGVER.EQ.6) THEN SIG_V = SIG_6(TT8,SS8) ELSEIF (SIGVER.EQ.8) THEN SIG_V = SIG_8(TT8,SS8) ENDIF ENDIF RETURN END REAL*4 FUNCTION SIGLOC_V(TT,SS,PRS,SIGVER) IMPLICIT NONE INTEGER SIGVER REAL*4 TT,SS,PRS C C SIGVER WRAPPER FOR SIGLOC C REAL*8 SS8,TT8,PRS8 REAL*8 SIGLOC_1,SIGLOC_2,SIGLOC_3,SIGLOC_4, & SIGLOC_5,SIGLOC_6,SIGLOC_7,SIGLOC_8 C TT8 = TT SS8 = SS PRS8 = PRS IF (MOD(SIGVER,2).EQ.1) THEN IF (SIGVER.EQ.1) THEN SIGLOC_V = SIGLOC_1(TT8,SS8,PRS8) ELSEIF (SIGVER.EQ.3) THEN SIGLOC_V = SIGLOC_3(TT8,SS8,PRS8) ELSEIF (SIGVER.EQ.5) THEN SIGLOC_V = SIGLOC_5(TT8,SS8,PRS8) ELSEIF (SIGVER.EQ.7) THEN SIGLOC_V = SIGLOC_7(TT8,SS8,PRS8) ENDIF ELSE IF (SIGVER.EQ.2) THEN SIGLOC_V = SIGLOC_2(TT8,SS8,PRS8) ELSEIF (SIGVER.EQ.4) THEN SIGLOC_V = SIGLOC_4(TT8,SS8,PRS8) ELSEIF (SIGVER.EQ.6) THEN SIGLOC_V = SIGLOC_6(TT8,SS8,PRS8) ELSEIF (SIGVER.EQ.8) THEN SIGLOC_V = SIGLOC_8(TT8,SS8,PRS8) ENDIF ENDIF RETURN END REAL*4 FUNCTION SOFSIG_V(RR,TT,SIGVER) IMPLICIT NONE INTEGER SIGVER REAL*4 RR,TT C C SIGVER WRAPPER FOR SOFSIG C REAL*8 RR8,TT8 REAL*8 SOFSIG_1,SOFSIG_2,SOFSIG_3,SOFSIG_4, & SOFSIG_5,SOFSIG_6,SOFSIG_7,SOFSIG_8 C RR8 = RR TT8 = TT IF (MOD(SIGVER,2).EQ.1) THEN IF (SIGVER.EQ.1) THEN SOFSIG_V = SOFSIG_1(RR8,TT8) ELSEIF (SIGVER.EQ.3) THEN SOFSIG_V = SOFSIG_3(RR8,TT8) ELSEIF (SIGVER.EQ.5) THEN SOFSIG_V = SOFSIG_5(RR8,TT8) ELSEIF (SIGVER.EQ.7) THEN SOFSIG_V = SOFSIG_7(RR8,TT8) ENDIF ELSE IF (SIGVER.EQ.2) THEN SOFSIG_V = SOFSIG_2(RR8,TT8) ELSEIF (SIGVER.EQ.4) THEN SOFSIG_V = SOFSIG_4(RR8,TT8) ELSEIF (SIGVER.EQ.6) THEN SOFSIG_V = SOFSIG_6(RR8,TT8) ELSEIF (SIGVER.EQ.8) THEN SOFSIG_V = SOFSIG_8(RR8,TT8) ENDIF ENDIF RETURN END REAL*4 FUNCTION TOFSIG_V(RR,SS,SIGVER) IMPLICIT NONE INTEGER SIGVER REAL*4 RR,SS C C SIGVER WRAPPER FOR TOFSIG C REAL*8 RR8,SS8 REAL*8 TOFSIG_1,TOFSIG_2,TOFSIG_3,TOFSIG_4, & TOFSIG_5,TOFSIG_6,TOFSIG_7,TOFSIG_8 C RR8 = RR SS8 = SS IF (MOD(SIGVER,2).EQ.1) THEN IF (SIGVER.EQ.1) THEN TOFSIG_V = TOFSIG_1(RR8,SS8) ELSEIF (SIGVER.EQ.3) THEN TOFSIG_V = TOFSIG_3(RR8,SS8) ELSEIF (SIGVER.EQ.5) THEN TOFSIG_V = TOFSIG_5(RR8,SS8) ELSEIF (SIGVER.EQ.7) THEN TOFSIG_V = TOFSIG_7(RR8,SS8) ENDIF ELSE IF (SIGVER.EQ.2) THEN TOFSIG_V = TOFSIG_2(RR8,SS8) ELSEIF (SIGVER.EQ.4) THEN TOFSIG_V = TOFSIG_4(RR8,SS8) ELSEIF (SIGVER.EQ.6) THEN TOFSIG_V = TOFSIG_6(RR8,SS8) ELSEIF (SIGVER.EQ.8) THEN TOFSIG_V = TOFSIG_8(RR8,SS8) ENDIF ENDIF RETURN END REAL*8 FUNCTION SIG_1(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_7term.h' SIG_1 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_1(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA0_7term.h' SIGLOC_1 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_1(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 INCLUDE '../include/stmt_fns_SIGMA0_7term.h' SOFSIG_1 = SOFSIG(RR8,TT8) END REAL*8 FUNCTION TOFSIG_1(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_7term.h' TOFSIG_1 = TOFSIG(RR8,SS8) END REAL*8 FUNCTION SIG_3(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_9term.h' SIG_3 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_3(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA0_9term.h' SIGLOC_3 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_3(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 INCLUDE '../include/stmt_fns_SIGMA0_9term.h' SOFSIG_3 = SOFSIG(RR8,TT8) END REAL*8 FUNCTION TOFSIG_3(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_9term.h' TOFSIG_3 = TOFSIG(RR8,SS8) END REAL*8 FUNCTION SIG_5(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_17term.h' SIG_5 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_5(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA0_17term.h' SIGLOC_5 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_5(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 REAL*8, PARAMETER :: TOL=1.D-6 INTEGER NN REAL*8 SN,SO REAL*8 SOFSIG_7 INCLUDE '../include/stmt_fns_SIGMA0_17term.h' C sofsig via Newton iteration from a 12-term 1st guess SN = SOFSIG_7(RR8,TT8) !non-negative DO NN= 1,10 SO = SN SN = SO - (SIG(TT8,SO)-RR8)/DSIGDS(TT8,SO) IF (NN.EQ.10 .OR. ABS(SN-SO).LT.TOL) THEN EXIT ENDIF ENDDO !nn SOFSIG_5 = SN END REAL*8 FUNCTION TOFSIG_5(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 REAL*8, PARAMETER :: TOL=1.D-6 INTEGER NN REAL*8 TN,TO REAL*8 TOFSIG_7 INCLUDE '../include/stmt_fns_SIGMA0_17term.h' C sofsig via Newton iteration from a 12-term 1st guess TN = TOFSIG_7(RR8,SS8) !non-negative DO NN= 1,10 TO = TN TN = TO - (SIG(TO,SS8)-RR8)/DSIGDT(TO,SS8) IF (NN.EQ.10 .OR. ABS(TN-TO).LT.TOL) THEN EXIT ENDIF ENDDO !nn TOFSIG_5 = TN END REAL*8 FUNCTION SIG_7(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_12term.h' SIG_7 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_7(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA0_12term.h' SIGLOC_7 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_7(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 INCLUDE '../include/stmt_fns_SIGMA0_12term.h' SOFSIG_7 = SOFSIG(RR8,TT8) END REAL*8 FUNCTION TOFSIG_7(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 INCLUDE '../include/stmt_fns_SIGMA0_12term.h' TOFSIG_7 = TOFSIG(RR8,SS8) END REAL*8 FUNCTION SIG_2(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_7term.h' SIG_2 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_2(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA2_7term.h' SIGLOC_2 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_2(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 INCLUDE '../include/stmt_fns_SIGMA2_7term.h' SOFSIG_2 = SOFSIG(RR8,TT8) END REAL*8 FUNCTION TOFSIG_2(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_7term.h' TOFSIG_2 = TOFSIG(RR8,SS8) END REAL*8 FUNCTION SIG_4(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_9term.h' SIG_4 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_4(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA2_9term.h' SIGLOC_4 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_4(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 INCLUDE '../include/stmt_fns_SIGMA2_9term.h' SOFSIG_4 = SOFSIG(RR8,TT8) END REAL*8 FUNCTION TOFSIG_4(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_9term.h' TOFSIG_4 = TOFSIG(RR8,SS8) END REAL*8 FUNCTION SIG_6(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_17term.h' SIG_6 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_6(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA2_17term.h' SIGLOC_6 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_6(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 REAL*8, PARAMETER :: TOL=1.D-6 INTEGER NN REAL*8 SN,SO REAL*8 SOFSIG_8 INCLUDE '../include/stmt_fns_SIGMA2_17term.h' C sofsig via Newton iteration from a 12-term 1st guess SN = SOFSIG_8(RR8,TT8) !non-negative DO NN= 1,10 SO = SN SN = SO - (SIG(TT8,SO)-RR8)/DSIGDS(TT8,SO) IF (NN.EQ.10 .OR. ABS(SN-SO).LT.TOL) THEN EXIT ENDIF ENDDO !nn SOFSIG_6 = SN END REAL*8 FUNCTION TOFSIG_6(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 REAL*8, PARAMETER :: TOL=1.D-6 INTEGER NN REAL*8 TN,TO REAL*8 TOFSIG_8 INCLUDE '../include/stmt_fns_SIGMA2_17term.h' C sofsig via Newton iteration from a 12-term 1st guess TN = TOFSIG_8(RR8,SS8) !non-negative DO NN= 1,10 TO = TN TN = TO - (SIG(TO,SS8)-RR8)/DSIGDT(TO,SS8) IF (NN.EQ.10 .OR. ABS(TN-TO).LT.TOL) THEN EXIT ENDIF ENDDO !nn TOFSIG_6 = TN END REAL*8 FUNCTION SIG_8(TT8,SS8) IMPLICIT NONE REAL*8 TT8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_12term.h' SIG_8 = SIG(TT8,SS8) END REAL*8 FUNCTION SIGLOC_8(TT8,SS8,PRS8) IMPLICIT NONE REAL*8 TT8,SS8,PRS8 INCLUDE '../include/stmt_fns_SIGMA2_12term.h' SIGLOC_8 = SIGLOC(TT8,SS8,PRS8) END REAL*8 FUNCTION SOFSIG_8(RR8,TT8) IMPLICIT NONE REAL*8 RR8,TT8 INCLUDE '../include/stmt_fns_SIGMA2_12term.h' SOFSIG_8 = SOFSIG(RR8,TT8) END REAL*8 FUNCTION TOFSIG_8(RR8,SS8) IMPLICIT NONE REAL*8 RR8,SS8 INCLUDE '../include/stmt_fns_SIGMA2_12term.h' TOFSIG_8 = TOFSIG(RR8,SS8) END C### C--------------------------------------------------------------------------- C --- routines not in ALLmpi/src/ncoda_archv_inc.f C--------------------------------------------------------------------------- C subroutine blkinr(rvar,cvar,cfmt) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c real rvar character cvar*6,cfmt*(*) c c read in one real value c character*6 cvarin c read(uoff+99,*) rvar,cvarin if (mnproc.eq.1) then write(lp,cfmt) cvarin,rvar call flush(lp) endif !1st tile c if (cvar.ne.cvarin) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkinr - input ',cvarin, + ' but should be ',cvar write(lp,*) call flush(lp) endif !1st tile call xcstop('(blkinr)') stop endif return end subroutine blkinr2(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c real rvar integer nvar character cvar1*6,cvar2*6,cfmt1*(*),cfmt2*(*) c c integer lp c common/linepr/lp c c read in one real value from stdin, c identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) c character*6 cvarin c read(uoff+99,*) rvar,cvarin c if (cvar1.eq.cvarin) then nvar = 1 if(mnproc.eq.1)then write(lp,cfmt1) cvarin,rvar call flush(lp) endif elseif (cvar2.eq.cvarin) then nvar = 2 if(mnproc.eq.1)then write(lp,cfmt2) cvarin,rvar call flush(lp) endif else if(mnproc.eq.1)then write(lp,*) write(lp,*) 'error in blkinr2 - input ',cvarin, + ' but should be ',cvar1,' or ',cvar2 write(lp,*) call flush(lp) endif call xcstop('(blkinr2)') stop endif return end subroutine blkinr9(rvar,nvar, & cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3, & cvar4,cfmt4,cvar5,cfmt5,cvar6,cfmt6, & cvar7,cfmt7,cvar8,cfmt8,cvar9,cfmt9) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c real rvar integer nvar character cvar1*6,cvar2*6,cvar3*6,cfmt1*(*),cfmt2*(*),cfmt3*(*) character cvar4*6,cvar5*6,cvar6*6,cfmt4*(*),cfmt5*(*),cfmt6*(*) character cvar7*6,cvar8*6,cvar9*6,cfmt7*(*),cfmt8*(*),cfmt9*(*) c c integer lp c common/linepr/lp c c read in one of nine possible real values from stdin, c identified as either cvar1 (return nvar=1) or c cvar2 (return nvar=2) or c ... c cvar9 (return nvar=9) c character*6 cvarin c read(uoff+99,*) rvar,cvarin c if (cvar1.eq.cvarin) then nvar = 1 if (mnproc.eq.1) then write(lp,cfmt1) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar2.eq.cvarin) then nvar = 2 if (mnproc.eq.1) then write(lp,cfmt2) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar3.eq.cvarin) then nvar = 3 if (mnproc.eq.1) then write(lp,cfmt3) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar4.eq.cvarin) then nvar = 4 if (mnproc.eq.1) then write(lp,cfmt4) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar5.eq.cvarin) then nvar = 5 if (mnproc.eq.1) then write(lp,cfmt5) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar6.eq.cvarin) then nvar = 6 if (mnproc.eq.1) then write(lp,cfmt6) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar7.eq.cvarin) then nvar = 7 if (mnproc.eq.1) then write(lp,cfmt7) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar8.eq.cvarin) then nvar = 8 if (mnproc.eq.1) then write(lp,cfmt8) cvarin,rvar call flush(lp) endif !1st tile elseif (cvar9.eq.cvarin) then nvar = 9 if (mnproc.eq.1) then write(lp,cfmt9) cvarin,rvar call flush(lp) endif !1st tile else if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkinr9 - input ',cvarin, + ' but should be:' write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or' write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or' write(lp,*) cvar7,' or ',cvar8,' or ',cvar9 write(lp,*) call flush(lp) endif !1st tile call xcstop('(blkinr9)') endif return end subroutine blkini(ivar,cvar) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c integer ivar character*6 cvar c c read in one integer value c character*6 cvarin c read(uoff+99,*) ivar,cvarin if (mnproc.eq.1) then write(lp,6000) cvarin,ivar call flush(lp) endif !1st tile c if (cvar.ne.cvarin) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkini - input ',cvarin, + ' but should be ',cvar write(lp,*) call flush(lp) endif !1st tile call xcstop('(blkini)') stop endif return 6000 format('blkini: ',a6,' =',i10) end subroutine blkini2(ivar,nvar,cvar1,cvar2) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c integer ivar,nvar character*6 cvar1,cvar2 c c integer lp c common/linepr/lp c c read in one integer value from stdin c identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) c character*6 cvarin c read(uoff+99,*) ivar,cvarin if (mnproc.eq.1) then write(lp,6000) cvarin,ivar call flush(lp) endif c if (cvarin.eq.cvar1) then nvar = 1 elseif (cvarin.eq.cvar2) then nvar = 2 else if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkini2 - input ',cvarin, + ' but should be ',cvar1,' or ',cvar2 write(lp,*) call flush(lp) endif call xcstop('(blkini2)') stop endif return 6000 format('blkini: ',a6,' =',i6) end subroutine blkini3(ivar,nvar,cvar1,cvar2,cvar3) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c integer ivar,nvar character*6 cvar1,cvar2,cvar3 c c integer lp c common/linepr/lp c c read in one of three possible integer values from stdin c identified as either cvar1 (return nvar=1) or c cvar2 (return nvar=2) or c cvar3 (return nvar=3) c call blkini9(ivar,nvar,cvar1,cvar2,cvar3, & 'XXXXXX','XXXXXX','XXXXXX', & 'XXXXXX','XXXXXX','XXXXXX') return end subroutine blkini9(ivar,nvar,cvar1,cvar2,cvar3, & cvar4,cvar5,cvar6, & cvar7,cvar8,cvar9) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c integer ivar,nvar character*6 cvar1,cvar2,cvar3 character*6 cvar4,cvar5,cvar6 character*6 cvar7,cvar8,cvar9 c c integer lp c common/linepr/lp c c read in one of nine possible integer values from stdin c identified as either cvar1 (return nvar=1) or c cvar2 (return nvar=2) or c ... c cvar9 (return nvar=9) c cvar3 (return nvar=3) or c character*6 cvarin c read(uoff+99,*) ivar,cvarin if (mnproc.eq.1) then write(lp,6000) cvarin,ivar call flush(lp) endif !1st tile c if (cvarin.eq.cvar1) then nvar = 1 elseif (cvarin.eq.cvar2) then nvar = 2 elseif (cvarin.eq.cvar3) then nvar = 3 elseif (cvar4.eq.cvarin) then nvar = 4 elseif (cvar5.eq.cvarin) then nvar = 5 elseif (cvar6.eq.cvarin) then nvar = 6 elseif (cvar7.eq.cvarin) then nvar = 7 elseif (cvar8.eq.cvarin) then nvar = 8 elseif (cvar9.eq.cvarin) then nvar = 9 else if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkini9 - input ',cvarin, + ' but should be:' write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or' write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or' write(lp,*) cvar7,' or ',cvar8,' or ',cvar9 write(lp,*) call flush(lp) endif !1st tile call xcstop('(blkini9)') endif return 6000 format('blkini: ',a6,' =',i6) end subroutine blkinl(lvar,cvar) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c logical lvar character*6 cvar c c read in one logical value c due to a SGI bug for logical I/O: read in an integer 0=F,1=T c character*6 cvarin integer ivar c read(uoff+99,*) ivar,cvarin lvar = ivar .ne. 0 if (mnproc.eq.1) then write(lp,6000) cvarin,lvar call flush(lp) endif !1st tile c if (cvar.ne.cvarin) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkinl - input ',cvarin, + ' but should be ',cvar write(lp,*) call flush(lp) endif !1st tile call xcstop('(blkinl)') stop endif return 6000 format('blkinl: ',a6,' =',l10) end subroutine blkinl_98(lvar,cvar) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c logical lvar character*6 cvar c c read in one logical value from unit 98 c due to a SGI bug for logical I/O: read in an integer 0=F,1=T c character*6 cvarin integer ivar c read(uoff+98,*) ivar,cvarin lvar = ivar .ne. 0 if (mnproc.eq.1) then write(lp,6000) cvarin,lvar call flush(lp) endif !1st tile c if (cvar.ne.cvarin) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in blkinl - input ',cvarin, + ' but should be ',cvar write(lp,*) call flush(lp) endif !1st tile call xcstop('(blkinl)') stop endif return 6000 format('blkinl: ',a6,' =',l10) end subroutine blkinc(line) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c character line*80 c character cvar*6,cfmt*(*) c c read in one line up to 80 characters long c c character*6 cvarin c read(uoff+99,'(a)') line c if (mnproc.eq.1) then c write(lp,cfmt) cvarin,rvar c call flush(lp) c endif !1st tile c c if (cvar.ne.cvarin) then c if (mnproc.eq.1) then c write(lp,*) c write(lp,*) 'error in blkinr - input ',cvarin, c + ' but should be ',cvar c write(lp,*) c call flush(lp) c endif !1st tile c call xcstop('(blkinc)') c stop c endif return end subroutine xcstop(cerror) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c character*(*), intent(in) :: cerror c c********** c* c 1) stop all processes. c c 2) all processes must call this routine. c use 'xchalt' for emergency stops. c c 3) parameters: c name type usage description c ---------- ---------- ------- ---------------------------- c cerror char*(*) input error message c* c********** c c shared memory version, just stop. c if (cerror.ne.' ') then write(lp,*) '**************************************************' write(lp,*) cerror write(lp,*) '**************************************************' call flush(lp) endif stop '(xcstop)' end subroutine xcstop subroutine xchalt(cerror) implicit none c integer i0,j0,nbdy,mnproc,lp,uoff common/xcm/ i0,j0,nbdy,mnproc,lp,uoff c character*(*), intent(in) :: cerror c c********** c* c 1) stop all processes. c c 2) only one processes need call this routine, i.e. it is for c emergency stops. use 'xcstop' for ordinary stops called c by all processes. c c 3) parameters: c name type usage description c ---------- ---------- ------- ---------------------------- c cerror char*(*) input error message c* c********** c c shared memory version, just stop. c if (cerror.ne.' ') then write(lp,*) '**************************************************' write(lp,*) cerror write(lp,*) '**************************************************' call flush(lp) endif stop '(xchalt)' end subroutine xchalt