************************************************************************ program !MCPG implicit real(a-h,o-z),integer(i-n) c************************************************************************ include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' dimension * en(MPTEN*(MAXNZ+1)), * lxxk(MXBDY*MAXNZ),lyyk(MXBDY*MAXNZ), * lxyk(MAXNB*MAXNZ),lyxk(MAXNB*MAXNZ),lsponge(MAXSP), * snxk(MAXNB*MAXNZ),snyk(MAXNB*MAXNZ),lok(4*MAXSID*MAXNZ), * lpbcwk(MAXSID*MAXNZ),lpbcek(MAXSID*MAXNZ), * ifxk(9*MAXSID*MAXNZ), ifpxk(5*MAXSID*MAXNZ), ifyk(9*MAXSID*MAXNZ) logical NEWRUN, non_stable common /param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20), * itherm,mlc,limp common /run/ nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch common /winds/ mtau,matau,tausc,atau,froude,icloud common /shapi/ nord,nshapu,nshaph,nshapt,mshapu,mshaph,mshapt common /vert/ zin(MAXNZ+1), hin(MAXNZ), t_in(MAXNZ+1), s_in(MAXNZ+1), * bint(MAXNZ), cint(MAXNZ), dzin(MAXNZ+1), sigma(MAXNZ) common /strech/ xs(MAXXS), alpha(MAXXS), beta(MAXXS) common /errors/ ioerr, nstep common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) common /baro_input/ n_def_cor, mod_scheme, mod_solver, BAR_DELTA, * BAR_DSINK, ibar_key, nbaro, rayl common /main/npt data ioin,iout,ioerr /1,2,2/, k15flag /0/ npten = MPTEN c...............getting input & output filenames from the command line call inout (ioin) c...............input model parameters. call model_input(npt) call model_memory (nxp, nyp, nz, npt) call make_iox (nxp, nyp, mask, iox, nlok, lok, nsponge, lsponge) if (nsponge.gt.MAXSP) then print*,'increase MAXSP' stop endif call scaset (iox,xm,ym,xp,yp,f,emx,emy,emxy,emx2,emy2,tp) call depth_init (npt, zin) call new_topo (nxp, nyp, nz, npt, zin, dzin, hin, nsig, sigma, dept, * h, nptk, nzi) call bndrys(npt,iox,tp,isxk,isyk,mask,h,nzi, * isk,iyk,lxxk,lyyk,lxyk,lyxk,snxk,snyk,lok, * lpbcwk,lpbcek,ifxk,ifpxk,ifyk,dept) call baro_dept(npt,nz,nzi,nzi_b,h,lxxk,lyyk,mbc,dept,tp,tp(npt2)) call data_init (npt,nptk,nz,isk,u,v,uc,vc,fu,fv,ft,fsal,bdiv,ubar,vbar,use_salt) c..............compute some geometry stuff call aarea (lxxk,lyxk,emx,emy,emxy,area,basin) if (irest .eq. 0) then NEWRUN = .true. call init_rstrt (nxp, nyp, nz, npt, zin, tp) else NEWRUN = irest .eq. 3 if (irest.eq.1 .or. irest.eq.2) CALL TIOS_CNTRL (3, 1) if (ibaro .ne. 0) ibaro = 1 call read_rstrt (nxp, nyp, nz, npt) endif c..............initialize Temperature/Salinity Climatology: call clim_init(npt,nstart,hin,sigma,dzin,h, * hclim,tclim,sclim,dclim,tp,nsponge,lsponge) if (irest .eq. 0) then call h_init (npt, nz, nstart, h, hclim) call temp_init (npt,nz,nzi,nstart,t_in,t,tclim) if (use_salt) call salt_init (npt,nz,nzi,nstart,s_in,sal,sclim) endif c..............initialize Heat/EP forcing: call hflx_init (nstart, npt, nxp, nyp, t, sst, cld, solr) if (use_salt) then call evpr_init (nstart, npt, sal, sss, ep) call dens_init (npt, nz, nzi, t, sal, dens, h) endif c..............initialize Wind forcing: call tau_init (nstart, npt, dtx, dty) dnt = delt * real(ncyc) * D2SEC dtmix = delt * real(limp) * D2SEC DLT_MIX = 2.0 * dtmix c..............initialize TIOS io-system call init_data_out (ftios, fbt, nxp, nyp, npt, xm, ym, en) c istep = 0 iday_new = int(nstart*delt) iday_curr = iday_new if (ibaro .ne. 0) then eps1 = 1./dnt eps2 = eps1 if (ibaro .eq. 2) eps1 = 0. call baro_sum (npt, nz, nzi_b, uc, vc, ubar, vbar) call baro_scale (npt, ubar, vbar, dept) mem0 = mem_get() call baro_init (iglob,eps1,nxp,nyp,nxyc,iox,nbxk,lxxk, * nbyk,lyxk,alon,blon,alat,blat,xm,ym,dept,dep_max,mgrid) write(iout, *)'Barotropic solver memory = ',mem_get()-mem0,' bytes.' if (ibaro .eq. 2) then call baro_tau (npt, uforc, vforc, taux, tauy) call baro_rhs (npt,emx,emy,emxy,lxxk,lyyk,lxyk,lyxk,snxk,snyk, * isyk,isk,mbc,lpbcwk,lpbcek,uforc,vforc,tp,tp(npt2),dept) call baro_solv (nxp,nyp,npt,iox,tp,uforc,vforc,psi) call curl_of_psi (npt,emx,emy,emxy,lxxk,lyyk,lxyk,lyxk,snxk,snyk, * isyk,isk,mbc,lpbcwk,lpbcek,psi,tp(npt1),tp(npt2),tp(npt3),dept) call baro_updat(npt,nz,nzi_b,h,uc,vc, * tp(npt1),tp(npt2),uforc,vforc,ubar,vbar) call decap (npt, nz, nzi, u,v,uc,vc,h) endif call baro_div (npt,emx,emy,emxy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,mbc,lpbcwk,lpbcek,ubar,vbar,bdiv,tp,dept) call baro_rinit (eps2) endif write(iout, *)'Total allocated memory = ',mem_get(),' bytes.' call bcset (mbc,lxxk,lyyk,npt,u,v,nzi,nzi_b) call baro_shap(nstep,npt,nz,nzi,nzi_b,dept,h,uc,vc,ubar,vbar,u,v,lxxk,lyyk) c...............compute w from the initial fields. call ddiv(npt,nzi_b,uc,vc,emx,emy,emxy,w,fhd,lxxk,lyyk,lxyk,lyxk,snxk,snyk, * isyk,isk,tp,mbc,lpbcwk,lpbcek,h,bdiv) call wtop (npt, w, fhd, fh, h, ncyc, istep, dnt) call dwcal (npt,nz,nsig,nzi,w,fhd,sigma) c......................................................................... c.....MAIN LOOP .......use an n-cycle Lorentz scheme for the timestep loop. call cpulog (fcpu, 0, iday_curr) DO NSTEP = NSTART, NLAST tenso = enso_start + enso_scale * nstep c if (nstep.eq.4) then c call data_out (tenso, nxp, nyp, npt, en) c stop c endif if (mtau .eq. 1) then call tau_lin(nstep,npt, ixd,im2d,blcf, taux,tauy,dtx,dty,tp) endif call clim_updt(npt,nz,nstep,hin,sigma,dzin,hclim,tclim,sclim,dclim) call qforc(nstep, npt, nxp,nyp,sst,cld,solr, tp, qb) if (use_salt) call epforc(nstep, npt, sal, sss, ep, qb) call vertu (npt,nz,nsig,nzi,nzi_b,bint,taux,tauy,u,v,w,h,fu,fv,fh, * vertx,verty,zfu,zfv) c print* c print*, t(6*npt+23), ft(6*npt+23), w(5*npt+23) if (use_salt) then call vertts (npt,nz,nzi,cint,q,qr,ep,w,h,t,ft,sal,fsal) else call vertt (npt,nz,nzi,cint,q,qr,w,h,t,ft) endif c print*, t(6*npt+23), ft(6*npt+23), w(5*npt+23) if (use_trac) call verttr (npt,nz,nzi,cint,w,h,tr,ftr) call dhoriz (npt,u,v,uc,vc,f,fu,fv,fhd,emx,emy,emxy,tp,mbc,zfu,zfv, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,lpbcwk,lpbcek,nzi_b, * corx,cory,xnl,ynl,fh) call btpgf (npt, nzi_b, h, t,dens,fu,fv,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,u,v,lpbcwk,lpbcek,zfu,zfv,pgfx,pgfy) call thoriz (npt,uc,vc,t,ft,fhd,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,mbc,lpbcwk,lpbcek) call capt (npt,nz,nzi,t,h) if (use_salt) then call thoriz (npt,uc,vc,sal,fsal,fhd,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,mbc,lpbcwk,lpbcek) call capt(npt,nz,nzi,sal,h) endif do i = 1, ntrac it = npt*nz*(i-1)+1 call thoriz (npt,uc,vc,tr(it),ftr(it),fhd,emx,emy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,lok,tp,mbc,lpbcwk,lpbcek) call capt (npt,nz,nzi,tr(it),h) enddo binv = dnt/(ncyc-istep) istep = mod(istep+1,ncyc) abinv = -(istep/dnt)*binv call baro_sum(npt,nz, nzi_b, fu,fv,u,v) call baro_scale (npt, u, v, dept) call fixed_dep(npt,nzi_b,h,fu,fv,u,v,rhsx,rhsy,crhsx,crhsy) call vel_updat (npt,nz,nzi_b,binv,abinv,uc,vc,fu,fv) if (ibaro .ne. 0) then call baro_comp(npt,dnt,abinv,binv,nbaro,uforc,vforc,u,v,zfu,zfv,dept) if (mod(nstep, nbaro*ncyc) .eq. 0) then call baro_rhs(npt,emx,emy,emxy,lxxk,lyyk,lxyk,lyxk,snxk,snyk, * isyk,isk,mbc,lpbcwk,lpbcek,uforc,vforc,u,tp,dept) call baro_solv (nxp,nyp,npt,iox,u,uforc,vforc,psi) call curl_of_psi (npt,emx,emy,emxy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,mbc,lpbcwk,lpbcek,psi,u,v,tp,dept) call baro_updat(npt,nz,nzi_b,h,uc,vc,u,v,uforc,vforc,ubar,vbar) call baro_div (npt,emx,emy,emxy,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,mbc,lpbcwk,lpbcek,ubar,vbar,bdiv,tp,dept) endif endif call bcset (mbc,lxxk,lyyk,npt,uc,vc,nzi,nzi_b) call shap_vec (nstep,npt,nz,uc,vc,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) call decap (npt, nz, nzi, u,v,uc,vc,h) if (nsig.gt.0) then call h_updat (npt,nsig, binv,abinv,h,fh) call hbcset (npt, nz, nsig, lok, h, hclim) call shap_scl(nstep,npt,nsig,h,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) endif call tupdat(npt,nz,nzi,binv,abinv,t,ft) call tbcset(npt,nz,lok,t_in,h,t,tclim) !reset temp at open boundaries call shap_scl(nstep,npt,nz,t,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) call tdecap (npt, nzi, t, h) if (use_salt) then call tupdat(npt,nz,nzi,binv,abinv,sal,fsal) call tbcset(npt,nz,lok,s_in,h,sal,sclim) call shap_scl(nstep,npt,nz,sal,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) call tdecap (npt, nzi, sal, h) endif do i = 1, ntrac it = npt*nz*(i-1)+1 call tupdat(npt,nz,nzi,binv,abinv,tr(it),ftr(it)) call tdecap (npt,nzi, tr(it), h) enddo if (imix.ne.0 .and. mod(nstep, limp).eq.0) then if (use_salt) call potn_dens (npt, nzi, t, sal, dens) if (imix .eq. 1) then !! Convective Adjustment call dconv (npt,nz,nzi,u,v,uc,vc,h,t,sal,dens,tr,convn) elseif (imix .eq. 2) then !! Ri Number Mixing call decap (npt, nz, nzi, u,v,uc,vc,h) call drich_mix (npt, nz, nzi, h, u,v,uc,vc,t,sal,dens) elseif (imix .eq. 3) then !! Combination of (1) & (2) call decap (npt, nz, nzi, u,v,uc,vc,h) call drich_mix (npt, nz, nzi, h, u,v,uc,vc,t,sal,dens) call dconv (npt,nz,nzi,u,v,uc,vc,h,t,sal,dens,tr,convn) elseif (imix .eq. 4) then !! Dake Chen's Scheme call decap (npt, nz, nzi, u,v,uc,vc,h) call comp_bncy(npt,nzi,dens,tp) call cvmix(npt,nzi,h,t,sal,tp,u,v) call jpmix(npt,nz,nzi,h,t,sal,tp,u,v) call ktmix(npt,nsig,dtmix,h,t,sal,tp,u,v,q,qr,ep,taux,tauy,sigma,uc) call capfrm(npt,nz,nzi,u,v,uc,vc,h) elseif (imix .eq. 5) then !! Dake Chen's Conv. Adjustment only call decap (npt, nz, nzi, u,v,uc,vc,h) call comp_bncy(npt,nzi,dens,tp) call cvmix(npt,nzi,h,t,sal,tp,u,v) call capfrm(npt,nz,nzi,u,v,uc,vc,h) endif endif call clim_relax (npt,nz,h,t,sal,hclim,tclim,sclim) if (use_salt) call situ_dens (npt, nz, nzi, t, sal, dens, h) call baro_shap(nstep,npt,nz,nzi,nzi_b,dept,h,uc,vc,ubar,vbar,u,v,lxxk,lyyk) call ddiv(npt,nzi_b,uc,vc,emx,emy,emxy,w,fhd,lxxk,lyyk,lxyk,lyxk, * snxk,snyk,isyk,isk,tp,mbc,lpbcwk,lpbcek,h,bdiv) call wtop (npt, w, fhd, fh, h, ncyc, istep, dnt) call dwcal (npt,nz,nsig,nzi,w,fhd,sigma) if (NEWRUN .or. (nergy.ne.0 .and. mod(nstep, nergy).eq.0)) then call knergy(npt,nz,nptk,isk,area,basin,h,u,v,en) call pnergy (NEWRUN,npt,nptk,nz,isk,h,area,t,dens,w,basin,en,tp) endif if (non_stable(iout, npt, nxp, nz, iox, t, u, v)) then call cpulog (fcpu, nstep, iday_curr) write (iout, *) 'Stable:ERROR, step', nstep, ' temp or velocity is bizarre' print*, 'Stable:ERROR, step', nstep, ' temp or velocity is bizarre' goto 333 else if ( mod(nstep,ncyc) .eq. 0 ) then iday_new = int(nstep*delt) call add_mean call keep_rstrt(nstep, nskip) if (NEWRUN .or. (iday_new .ne. iday_curr)) then iday_curr = iday_new call data_out (tenso, nxp, nyp, npt, en) call cpulog (fcpu, nstep, iday_curr) if (NEWRUN) NEWRUN = .false. endif endif endif call flush(iout) if ( first_step ) first_step = .false. c.....END of the MAIN LOOP ENDDO goto 444 c.............ABnormally finished run: 333 CALL TIOS_CNTRL (1, 1) call data_out (tenso, nxp, nyp, npt, en) c.............normally finished run: 444 call close_rstrt write (iout, *) 'Finished at step =', nstep call enso2date (tenso, id, im, iy) write (iout, *) ' <', tenso, '>' write (iout, *) ' <',id,':',im,':',iy,'>' stop end