************************************************************************ c Lamont Ocean AML Model (LOAM), Version 2.0 implicit real(a-h,o-z),integer(i-n) c************************************************************************ include 'comm_para.h' include 'comm_main.h' include 'comm_new.h' include 'comm_data.h' include 'diffiso.h' include 'comm_tracer.h' include 'comm_nao.h' include 'comm_bio.h' include 'comm_tios.h' logical non_stable dimension en(MPTEN*(MAXNZ+1)) dimension f_lor(0:3) c after a complete lorenz 4-cycle for a LINEAR scheme c u(n+1) = u(n) + delt * (f1 -f2 + f3) c where fi denotes the forcing for istep=i data f_lor /0.,1.,-1.,1./ data ioin,iout,ioerr /1,2,2/ call init_all flor = f_lor(0) iday_new = int(nstart*delt) iday_curr = iday_new write(iout, *)'Total allocated memory = ',mem_get(),' bytes.' c...............compute w from the initial fields. call ddiv (nxp,nyp,nz,npt,nzi,uc,vc,w,fhd,isk,tp,h,hupi, * dx2i,dy2i,ycos,inx,iny) call wtop (npt, w, fhd, fh, h, ncyc, istep, dnt) call dwcal (npt,nz,nsig,nzi,w,fhd,sigma) print*,'LOAM: beginning time loop' c......................................................................... c.....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 tmonth = mod(tenso,12.) ! month of the year, [0,12) if (tmonth.lt.0) tmonth = tmonth + 12. call set_tanom(tenso) call set_ramp(delt,nstep,w_day) binv = dnt/(ncyc-istep) istep = mod(istep+1,ncyc) abinv = -(istep/dnt)*binv if (use_trac) call force_tracer(tenso,npt,nz,ntrac,nstep,nxp,nyp, * iv_bot,rjuljar,juljar,dnt,nzi,tr,ftr,t,h,wnd,sal,pdens,ym,iox,tp) if (diffiso_alpha.gt.0) then if (i_slope_step.eq.0) then call cons_shap_dens(nstep,npt,nz,1,nz,h,pdens,tp,lxxk,lyyk,snxk, * snyk,isyk,isk,lpbcwk,lpbcek,saxmk,saxpk,saymk,saypk,sxp,syp) call shap_scl(nstep,npt,nz,1,nz,pdens,lxxk,lyxk,isyk,isk, * ifxk,ifpxk,ifyk,tp) call slope(npt,pdens,dens,nz,nzi,h,lxxk,lyyk,lxyk,lyxk,snxk, * snyk,isyk,isk,lok,tp,lpbcwk,lpbcek,saxmk,saxpk,saymk,saypk,xmld) endif i_slope_step = mod(i_slope_step+ 1, nslope) endif if (mtau .eq. 1) then call tau_lin(nstep, inao,tscl_anom, npt,iox,xm,ym, * taux,tauy,dtx,dty,taux_anom,tauy_anom,tmonth) elseif (mtau.eq.4) then call tau_special(nstep,nstart,delt,npt,taux,tauy) endif call clim_updt(npt,nz,nstep,sigma,dzin,hclim,tclim,sclim,dclim) if (use_ice) then call amlice_flux(nstep, delt, tmonth, tscl_anom, npt, nxp, nyp, * sst, cld, solr, wnd, sal, sss, prcp, qb, flor) else if (initq.eq.4) then call qforc_special(nstep,nstart,delt,npt,q,qr,qsr) else call qforc(nstep, npt, nxp,nyp, sst,cld,solr, wnd, qb) endif call epforc(nstep, npt, sal, sss, prcp, qb) endif if (ihfprt .gt. 0) then call hflx_pert(npt,nz,nxp,nyp,nstep,ym) endif if (use_dyice) then call dyice (nxp,nyp,npt,ntrac_sur,trs,uice,vice,pice, * etaC,etaB,wnd(npt+1),wnd(2*npt+1),u,v,tauwx,tauwy) endif if (use_bio) call force_bio(npt,nz,nstep,flor,tmonth) c if (use_dyice) then c call force_u_dyice (npt,taux,tauy,trs(i_cice),tauwx,tauwy,fu,fv) c else call force_u (npt,taux,tauy,fu,fv) c endif call pgfu(npt,nzi,h,dens,fu,fv,saxmk,saxpk,saymk,saypk,lxxk,lyyk, * lxyk,lyxk,snxk,snyk,isyk,isk,lok,tp,dpres,ddens,lpbcwk,lpbcek, * pgfx,pgfy,gradE) call vertu (npt,nz,nsig,nzi,bint,u,v,w,h,fu,fv,fh) call horizu (nxp,nyp,nz,npt,nsig,u,v,uc,vc,f,fu,fv,fhd,tp,isk, * lxxk,lyyk,nzibx,nziby,corx,cory,fh,hupi,dx2i,dy2i,ycos,inx,iny) call vertts (npt,nz,nzi,cint,q,qr,ep,w,h,t,ft,sal,fsal, * ft_va,ft_vd,fsal_va,fsal_vd,flor) if (use_river) call river_flux (npt,nz,nzi,tmonth,emx,emy, * fsal,fsal_ri,flor) call horizt (nxp,nyp,nz,npt,h,uc,vc,u,v,t,ft,isk,tp, * 1,ft_ha,flor,hupi_ts,dx2i,dy2i,ycos,inx,iny) call horizt (nxp,nyp,nz,npt,h,uc,vc,u,v,sal,fsal,isk,tp, * 1,fsal_ha,flor,hupi_ts,dx2i,dy2i,ycos,inx,iny) if (use_trac.or.use_bio) then call verttr (npt,nz,ntrac_tot,nzi,cint,w,h, * tr,ftr,ftr_va,ftr_vd,flor) do i = 1, ntrac_tot it = npt*nz*(i-1)+1 it1 = npt*(nz+1)*(i-1)+1 call horizt (nxp,nyp,nz,npt,h,uc,vc,u,v,tr(it),ftr(it),isk,tp, * if_tr_adv,ftr_ha(it1),flor,hupi_tr,dx2i,dy2i,ycos,inx,iny) enddo endif if (use_dyice_old) call ice_adv(nxp,nyp,npt,u,v, * trs(i_cice),trs(i_hice),trs(i_thice),isk,tp, * fcice,fhice,fthice,dx2i,dy2i,ycos,inx,iny) c add momentum diffusion, if desired if (use_modiff) then call diff_hor(diff_coef_mo,diff_dz_mo,npt,nz,nzi,h,u,fu, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk) call diff_hor(diff_coef_mo,diff_dz_mo,npt,nz,nzi,h,v,fv, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk) endif if (use_div_damp) then call diff_div(diff_coef_dd,npt,nz,nzi,fhd,fu,fv, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk,nxp,nyp,iox,mask,isteer) endif call shapu_diff (dtshap,npt,nz,u,v,h,fu,fv,lxxk,lyxk, * isyk,isk,ifxk,ifpxk,ifyk,tp) c add tracer diffusion, if desired if (use_trdiff) then c update variable diffusion coefficient, when desired if (initkmap.eq.5) then if (i_kmap_step.eq.0) then call Kmap_compute (npt, nz, nxp, nyp) endif i_kmap_step= mod(i_kmap_step+ 1, nkmap) endif call diff_isogm(diffiso_alpha,diff_coef_tr,npt,nz,nzi,h,t,ft, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk,ft_hd,ft_vd,flor,xmld) call diff_isogm(diffiso_alpha,diff_coef_tr,npt,nz,nzi,h,sal,fsal, * lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk,fsal_hd,fsal_vd,flor,xmld) do i = 1, ntrac_tot it = npt*nz*(i-1)+1 it1 = npt*(nz+1)*(i-1)+1 call diff_isogm(diffiso_alpha,diff_coef_tr,npt,nz,nzi,h, * tr(it),ftr(it),lxxk,lyyk,lxyk,lyxk,snxk,snyk,isyk,isk,tp,lpbcwk, * lpbcek,saxmk,saxpk,saymk,saypk,ftr_hd(it1),ftr_vd(it1),flor,xmld) enddo endif if (use_dyice_old) * call ice_diff(npt,nz,trs(i_cice),trs(i_hice),trs(i_thice), * fcice,fhice,fthice, * lxxk,lyyk,snxk,snyk,isyk,isk,tp,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk) call bc_relax (dnt,bc_coef,mbc,lxxk,lyyk,npt,uc,vc,fu,fv,nzi) call baro_split(npt,nz,nzi,h,dept,fu,fv,u,v,tfu,tfv) c call bc_steer(npt,fu,fv,nzi,lxxk,lyyk,iox,mask,isteer,saxpk,saypk) call vel_updat (t_fac,npt,nz,nzi,binv,abinv,uc,vc,fu,fv) if (ibaro .ne. 0) then call baro_comp(npt,nz,nzi,uforc,vforc,tfu,tfv, * corx,cory,pgfx,pgfy,cor,pgf,gradE) call baro_rhs (nxp,nyp,npt,isk,uforc,vforc,rhs,tp,gradE,gradH, * binv,dt_sub,cor,pgf,dx2i,dy2i,ycos,inx,iny,lxxk,lyyk) call baro_solv (nxp,nyp,npt,iox,rhs,uforc,vforc,psin) call curlH (nxp,nyp,npt,psin,dept,dx2i,dy2i,ycos,inx,iny,u,v) call baro_updat1(npt,nz,nzi,h,binv,abinv,uc,vc,u,v,ubar,vbar) call make_small(nxp,nyp,npt,psi,psin,inx,iny,isk) ! for output and save file only endif c****enforce boundary conditions call bcset (mbc,lxxk,lyyk,npt,uc,vc) call decap (npt, nz, nzi, u,v,uc,vc,h) c****update the tracer fields if (idens.ne.0) then call tupdat(ttr_fac,npt,nz,nzi,binv,abinv,h,t,ft) call tupdat(ttr_fac,npt,nz,nzi,binv,abinv,h,sal,fsal) endif do i = 1, ntrac_tot it = npt*nz*(i-1)+1 call tupdat(ttr_fac,npt,nz,nzi,binv,abinv,h,tr(it),ftr(it)) enddo if (use_dyice) then call ice_limit (100., npt, t, trs(i_cice), trs(i_hice), trs(i_thice)) endif if (use_dyice_old) then call ice_limit (8., npt, t, trs(i_cice), trs(i_hice), trs(i_thice)) call ice_updat(npt,binv,abinv,h, * trs(i_cice),trs(i_hice),trs(i_thice),fcice,fhice,fthice) endif if (.not.use_ice) call t_limit (npt, nzi, -1.7, 100, t) c****filter the tracer fields if (idens.ne.0) then call shap_scl(nstep,npt,nz,1,nz,t,lxxk,lyxk,isyk,isk, * ifxk,ifpxk,ifyk,tp) call cons_shap_scl(nstep,npt,nz,1,nz,h,t,tp,lxxk,lyyk,snxk, * snyk,isyk,isk,lpbcwk,lpbcek,saxmk,saxpk,saymk,saypk,sxp,syp) call shap_scl(nstep,npt,nz,1,nz,sal,lxxk,lyxk,isyk,isk, * ifxk,ifpxk,ifyk,tp) call cons_shap_scl(nstep,npt,nz,1,nz,h,sal,tp,lxxk,lyyk, * snxk,snyk,isyk,isk,lpbcwk,lpbcek,saxmk,saxpk,saymk,saypk, * sxp,syp) endif if (use_trac.or.use_bio) then do i = 1, ntrac_tot it = npt*nz*(i-1)+1 call shap_scl(nstep,npt,nz,1,nz,tr(it),lxxk,lyxk, * isyk,isk,ifxk,ifpxk,ifyk,tp) call cons_shap_scl(nstep,npt,nz,1,nz,h,tr(it),tp,lxxk, * lyyk,snxk,snyk,isyk,isk,lpbcwk,lpbcek,saxmk,saxpk, * saymk,saypk,sxp,syp) enddo endif call situ_dens (npt, nz, nzi, t, sal, dens, h) call potn_dens (npt, nz, nzi, t, sal, pdens) call mixing(i_mix_step,dtmix,mix_count,bint,cint) call clim_relax (npt,nz,ntrac_tot,h,t,sal,u,v,ubar,vbar, * hclim,tclim,sclim,tr,trclim,flor,dt_sub,ft_sp,fsal_sp,ftr_sp) c call shap_vec(nstep,npt,nz,u,v,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) call baro_reset(nstep,npt,nz,nzi,nzibx,nziby,dept,h,ubar,vbar,u,v, * lxxk,lyyk,uc,vc) call ddiv (nxp,nyp,nz,npt,nzi,uc,vc,w,fhd,isk,tp,h,hupi, * dx2i,dy2i,ycos,inx,iny) call wtop (npt, w, fhd, fh, h, ncyc, istep, dnt) call dwcal (npt,nz,nsig,nzi,w,fhd,sigma) flor = f_lor(istep) 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 (use_trac) call tr_integrate(npt,nz,nstep,nstart,nzi,area, * tr,h,tramt,trfirst,vollev,trlev,sqrlev,itrac_cons) if (non_stable(iout, npt, nxp, nz, nzi, iox, t, u, v)) then call cpulog (fcpu, nstep, iday_curr) write (iout, *) 'Stable:ERROR, step', nstep, ' t,u or v is crazy' print*, 'Stable:ERROR, step', nstep, ' t,u or v is crazy' goto 333 else if ( istep .eq. 0 ) then ! end of Lorenz cycle iday_new = int(nstep*delt) call add_mean(nstep, npt, nxp, nyp) call keep_rstrt(1,nstep, nskip) call keep_rstrt(0,nstep, nskip) call data_out (tenso, nxp, nyp, npt, nz, en) if (NEWRUN .or. (iday_new .ne. iday_curr)) then iday_curr = iday_new call cpulog (fcpu, nstep, iday_curr) if (NEWRUN) NEWRUN = .false. call misc_out (tenso) endif call bud_reset endif endif call baro_updat2(npt,nz,binv,abinv,cor,pgf,gradE,tfu,tfv) 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, nz, 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