program ocean c c======================================================================= c === c OCEAN is the primary calling routine. It performs all === c operations which need be done only once at the === c beginning of each run of the model, calls STEP === c once per timestep, and attends to operations === c which must be done only at the end of each run. === c === c======================================================================= c === c This code uses as a base the GFDL (Cox,1984) PE model. === c === c * Open boundaries on all sides === c * Shapiro filter in horizontal replaces Laplacian diffusion === c * The model domain may be rotated with respect to Cartesian === c system of coordinates === c Michael A. Spall 13/01/87 === c * Hybrid coordinates in the vertical === c Michael A. Spall 18/11/87 === c * Intermittent OI Assimilation === c Carlos J. Lozano 03/05/92 === c Hernan G. Arango 01/07/92 === c * Lagrangian drifters trajectories === c Carlos J. Lozano === c Patrick J. Haley === c Hernan G. Arango 30/07/93 === c * Modularized with C-preprocessing options === c Hernan G. Arango 01/11/93 === c * Generalized terrain-following vertical coordinates === c Patrick J. Haley === c Carlos J. Lozano === c Hernan G. Arango 01/11/93 === c * Open boundary conditions === c Carlos J. Lozano === c N. Quinn Sloan === c Hernan G. Arango 01/11/93 === c * New hydrostatic pressure gradient formulation === c Patrick J. Haley === c Carlos J. Lozano 06/01/94 === c * Coasts and Islands setup === c Carlos J. Lozano 27/01/94 === c * I/O Data Management using NetCDF === c Hernan G. Arango 25/02/94 === c Hernan G. Arango 17/08/94 === c * McGillicuddy et al. (1995) biological model === c Carlos J. Lozano === c Hernan G. Arango 06/07/94 === c * Vertical mixing revisited: Mixed-layer parameterization === c Laurence A. Anderson 02/03/95 === c Carlos J. Lozano 06/12/99 === c Patrick J. Haley 06/12/99 === c * Reinstated Spherical Coordinates === c Patrick J. Haley 15/09/95 === c * Biological Model revisited: === c (1) McGillicuddy et al. (1990) biological model === c (2) Anderson expansion of McGillicuddy biological model === c (3) Fasham et al. (1990) biological model === c Laurence A. Anderson 15/09/95 === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef oias # include #endif #include #include #include c c----------------------------------------------------------------------- c Define local and equivalence local data c----------------------------------------------------------------------- c integer i,j,k,l,m,n,nbuf #ifdef oias integer idisk,idiska,idiskb,jdisk,jdiska #endif #ifdef sunfpe integer ieeer,my_handler,ieee_handler #endif FLOAT * fkmp(imt,jmt),fkmq(imt,jmt),fins(ndices),tnext equivalence (cfn,fins),(cfs,fkmp),(cfe,fkmq) #ifdef sunfpe external my_handler #endif c c======================================================================= c Begin executable code. c======================================================================= c #ifdef sunfpe c----------------------------------------------------------------------- c Set trap for floating point execeptions. c----------------------------------------------------------------------- c ieeer=ieee_handler('set','common',my_handler) if(ieeer.ne.0) * write(stdout,*) ' IEEE_HANDLER cannot set my_handler' c #endif #if ( defined nest2larger | defined nest2smaller ) & !defined nestnultest c----------------------------------------------------------------------- c Initialize two-way nesting communications. c----------------------------------------------------------------------- c call nest_init c #endif c----------------------------------------------------------------------- c Write header to echo file. c----------------------------------------------------------------------- c call headln (stdout) c c======================================================================= c Begin introductory section which is needed for each run ============ c (including re-starts) ============ c======================================================================= c c----------------------------------------------------------------------- c Read in input parameters c----------------------------------------------------------------------- c call inparm c #if defined nest2smaller & !defined nestnultest c----------------------------------------------------------------------- c Spawn two-way nested sub-domain. c----------------------------------------------------------------------- c call nest_spawn c #endif c----------------------------------------------------------------------- c Call INFLD to get model domain information #ifdef oias c and read in OI assimilation parameters #endif c----------------------------------------------------------------------- c itt=-1 call infld(0) c #if ( defined nest2larger | defined nest2smaller ) & !defined nestnultest c----------------------------------------------------------------------- c Transfer domain & flag information c----------------------------------------------------------------------- c call nest_domain call nest_flags c #endif c----------------------------------------------------------------------- c Check C-Preprocessing options and input parameters. c----------------------------------------------------------------------- c call chkparm c c----------------------------------------------------------------------- c Set horizontal and vertical grid c----------------------------------------------------------------------- c call hvgrid c c----------------------------------------------------------------------- c Initialize output NetCDF files: define dimensions, variables, and c attributes. Write out geometry fields. (Only if number of output c levels NLEV is greater than zero). c----------------------------------------------------------------------- c wnetcdf=.false. if(nlev.gt.0) then wnetcdf=.true. call defout endif nlast=nlast+1 c c----------------------------------------------------------------------- c Open the disk datasets c----------------------------------------------------------------------- c call ostart(kontrl,20,20,1) call ostart(kflds,nkflds*nwds,nwds,1) nbuf=2 call ostart(labs(1),jmt*nslab,nslab,nbuf) call ostart(labs(2),jmt*nslab,nslab,nbuf) call ostart(labs(3),jmt*nslab,nslab,nbuf) c c If re-start mode, advance INFLD file to the beginning of the c boundary conditions record and then read re-start file c if(nfirst.eq.0)then itt=0 ttsec=c0 do 100 j=1,jmt #ifdef analytical call setvert(j) #endif call infld(j) 100 continue call ord(21) endif c c======================================================================= c End introductory section =========================================== c======================================================================= c c======================================================================= c Begin section which is executed only when starting ================= c a run from scratch ================= c======================================================================= c if(nfirst.eq.0) goto 260 c c----------------------------------------------------------------------- c Set boundary indicators and auxiliary arrays. c----------------------------------------------------------------------- c call bdryindx(fkmp,fkmq) c c Set time counters prior to call INFLD c itt=0 ttsec=c0 #ifdef oias goto 120 c c Upgrade permuting disk I/O units on an assimilation timestep c 110 continue idiskb=mod(itt+1+2,3)+1 idisk =mod(itt+1+0,3)+1 idiska=mod(itt+1+1,3)+1 c c If assimilating, fetch P field c if(ioip.eq.1) then jdiska=mod(itt+1,3)+1 jdisk =mod(itt+0,3)+1 call oget(kflds,nwds,(jdiska-1)*nwds+1,p) endif 120 continue #endif c c----------------------------------------------------------------------- c Initialize slab data on disk c----------------------------------------------------------------------- c do 180 j=1,jmt #ifdef oias c c Read in J slab on an assimilation step c if((noi.gt.0).and.(ittoi.eq.2)) then call oget(labs(idiskb),nslab,(j-1)*nslab+1,tb ) call oget(labs(idisk ),nslab,(j-1)*nslab+1,t ) else do 130 i=1,imt c c Set maximum level indicators to values computed above c fkmt(i)=fkmp(i,j) fkmu(i)=fkmq(i,j) 130 continue endif #else do 130 i=1,imt c c Set maximum level indicators to values computed above c fkmt(i)=fkmp(i,j) fkmu(i)=fkmq(i,j) 130 continue #endif c c----------------------------------------------------------------------- #ifdef oias c Read in initial conditions or begin row-by-row assimilation #else c Read in initial conditions #endif c----------------------------------------------------------------------- c #if defined analytical | defined oias call setvert(j) #endif call infld(j) #ifdef oias c c Move TAU data to TAU-1 level on an assimilation timestep c if((noi.gt.0).and.(ittoi.eq.2)) then if(ioiuv.eq.1) then do 140 k=1,km do 140 i=1,imt ub(i,k)=u(i,k) vb(i,k)=v(i,k) 140 continue endif do 160 m=1,nt if(ioits(m).eq.1) then do 150 k=1,km do 150 i=1,imt c tb(i,k,m)=t(i,k,m) 150 continue endif 160 continue else #endif do 170 k=1,km do 170 i=1,imt #ifdef coast if (fkmu(i).lt.c1) then u(i,k)=c0 v(i,k)=c0 endif #endif ub(i,k)=u(i,k) vb(i,k)=v(i,k) do 170 m=1,nt #ifdef coast if (fkmt(i).lt.c1) then t(i,k,m)=c0 endif #endif tb(i,k,m)=t(i,k,m) 170 continue #ifdef oias endif c c Send the initial or assimilation slabs to disk c if((noi.gt.0).and.(ittoi.eq.2)) then call oput(labs(idiskb),nslab,(j-1)*nslab+1,tb) call oput(labs(idisk ),nslab,(j-1)*nslab+1,t ) else call oput(labs(1),nslab,(j-1)*nslab+1,tb) call oput(labs(2),nslab,(j-1)*nslab+1,t ) call oput(labs(3),nslab,(j-1)*nslab+1,t ) endif #else c c Send the initial slabs to disk c call oput(labs(1),nslab,(j-1)*nslab+1,tb) call oput(labs(2),nslab,(j-1)*nslab+1,t ) call oput(labs(3),nslab,(j-1)*nslab+1,t ) #endif 180 continue #ifdef oias if((noi.gt.0).and.(ittoi.eq.2)) then if(ioip.eq.1) then call oput(kflds,nwds,(jdiska-1)*nwds+1,p) endif endif c c Skip remainder of disk initialization on an assimilation timestep c if((noi.gt.0).and.(ittoi.eq.2)) goto 380 #endif c c----------------------------------------------------------------------- c Initialize remainder of disk c----------------------------------------------------------------------- c c Set streamfunction at TAU-1 to initial streamfunction c Set ZTD and ZTDB to 0.0 (initialize time change of P to 0.0 on disk) c do 190 j=1,jmt do 190 i=1,imt pb(i,j)=p(i,j) ztd(i,j)=c0 ztdb(i,j)=c0 190 continue c c Send initialized data to disc c call oput(kflds,nwds,1,pb) call oput(kflds,nwds,nwds+1,p) call oput(kflds,nwds,3*nwds+1,hr) call oput(kflds,nwds,4*nwds+1,ztd) call oput(kflds,nwds,5*nwds+1,ztd) c c Convert start and end indices to real (this is done to accomodate c the option of running the model in halfword mode. c n=0 do 200 l=1,lseg do 200 j=1,jmt n=n+1 fins(n)=FLoaT(isz(j,l)) 200 continue do 210 l=1,lseg do 210 j=1,jmt n=n+1 fins(n)=FLoaT(iez(j,l)) 210 continue #ifdef islands do 220 m=1,misle n=n+1 fins(n)=FLoaT(isis(m)) 220 continue do 230 m=1,misle n=n+1 fins(n)=FLoaT(ieis(m)) 230 continue do 240 m=1,misle n=n+1 fins(n)=FLoaT(jsis(m)) 240 continue do 250 m=1,misle n=n+1 fins(n)=FLoaT(jeis(m)) 250 continue #endif c c Send remainder of initialized data to disc c call oput(kontrl,20,1,itt) call oput(kflds,ndices,6*nwds+1,fins) c #ifdef peprf c----------------------------------------------------------------------- c Read in input parameters for sampling profiles. c----------------------------------------------------------------------- c call readprf c #endif c======================================================================= c End section to start from scratch ================================== c======================================================================= c c======================================================================= c Begin section to timestep the model ================================ c======================================================================= c c----------------------------------------------------------------------- c Read disk data into memory for startup c----------------------------------------------------------------------- c c Read in timestep counter, total elapsed time, area and volume c 260 call oget(kontrl,20,1,itt) c c If re-start mode, reset time to force reading of boundary conditions c in the first call to INFLD in STEP c if(nfirst.eq.0)then write(stdout,*) ' Re-start mode - TNEXT, TTSEC = ',tnext,ttsec tnext=ttsec-p5 write(stdout,*) ' TNEXT reset to = ',tnext endif c c Read in start and end indices and convert to integers c call oget(kflds,ndices,6*nwds+1,fins) n=0 do 270 l=1,lseg do 270 j=1,jmt n=n+1 isz(j,l)=int(fins(n)) 270 continue do 280 l=1,lseg do 280 j=1,jmt n=n+1 iez(j,l)=int(fins(n)) 280 continue #ifdef islands do 290 m=1,misle n=n+1 isis(m)=int(fins(n)) 290 continue do 300 m=1,misle n=n+1 ieis(m)=int(fins(n)) 300 continue do 310 m=1,misle n=n+1 jsis(m)=int(fins(n)) 310 continue do 320 m=1,misle n=n+1 jeis(m)=int(fins(n)) 320 continue #endif c c Compute permuting disc indicators and read in 2 levels of c stream function as well as reciprocal depth. c ndisk =mod(itt+0,3)+1 ndiska=mod(itt+1,3)+1 call oget(kflds,nwds,(ndisk-1)*nwds+1,pb) call oget(kflds,nwds,(ndiska-1)*nwds+1,p) call oget(kflds,nwds,3*nwds+1,hr) c c If re-start, initialize the remainder of the /HYBRID/ common block c and build array HV from its inverse depth. c if(nfirst.eq.0) then do 330 j=1,jmt do 330 i=1,imt hd(i,j)=c1/hr(i,j) 330 continue do 340 j=1,jmt do 340 i=1,imt hdv(i,j)=hv(i,j) hv(i,j)=c1/hv(i,j) 340 continue c c If re-start, set boundary indicators auxiliary arrays. c call bdryindx(fkmp,fkmq) #ifdef rmdenbar c c If re-start, calculate mean background density profile c call meanrho #endif endif c c----------------------------------------------------------------------- c Initialize several variables to zero to avoid an uninitialized c variable type of error later where, for purposes of vectorization, c the computation proceeds across land points c----------------------------------------------------------------------- c do 370 i=1,imt do 350 j=1,jmt ztd(i,j)=c0 350 continue do 360 k=1,kmp1 tempa(i,k)=c0 tempb(i,k)=c0 360 continue 370 continue #ifdef oias c c End of disk initialization on an assimilation timestep. c 380 continue #endif c c----------------------------------------------------------------------- c Proceed with timestepping until the specified number c of steps have been taken c----------------------------------------------------------------------- c 390 call step #ifdef ldrifters c c----------------------------------------------------------------------- c Compute Lagrangian trajectories c----------------------------------------------------------------------- c call drifters #endif #ifdef extraction c c----------------------------------------------------------------------- c Extract information for sub-domain(s) and write to NetCDF file(s) c----------------------------------------------------------------------- c call xtrsubdom #endif c c----------------------------------------------------------------------- c Save restart data every NWRITE timesteps, and at end of run if NA=1 c----------------------------------------------------------------------- c if(mod(itt,nwrite).eq.0) call owrt(22) if((itt.eq.nlast).and.(na.eq.1)) call owrt(21) #ifdef oias c c Activate assimilation switch ITTOI on an assimilation timestep c if((noi.gt.0).and.(ittoi.eq.1)) then ittoi=2 write(stdout,*) ' Begin assimilation cycle = ',icoi, * ' at day = ',ttsec*sec2day endif if((noi.gt.0).and.(ittoi.eq.2)) goto 110 #endif #ifdef sunflush c c----------------------------------------------------------------------- c Flush output buffers. c----------------------------------------------------------------------- c call flush(stdout) #endif c if(itt.lt.nlast) goto 390 c c======================================================================= c End timestepping of the model ====================================== c======================================================================= c call exitus('OCEAN DONE') end