subroutine anaflds c c======================================================================= c === c This subroutine provides Analytical fields to the PE model. The === c Analytical fields are defined by user. === c === c Entries: === c === c ANADINFO domain information and initialization parameters. === c ANATOPO bottom topography. === c ANAUVPT initial conditions for internal mode velocity === c components, transport streamfunction, and tracers. === c ANAMRHO mean density. === c ANABNDRY lateral boundary conditions. === c DPTHBNDRY depths of boundary points. === c ANAVBC vertical boundary conditions. === c ANADIAG diagnostics. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #if defined bioAnder | defined bioFasham | defined bioMcGillic | defined bioDuse # include #endif #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c logical first integer i,ibdy,j,jj,k,m,kz FLOAT * a0,a1,b0,b1,delt,depth,sbfscl,sqfscl,ssur,swfscl,tsur,uvmag #ifdef forcing FLOAT * fx #endif #ifdef bioMcGillic FLOAT * cnh41,cnh42,cno31,cno32,cphy1,cphy2,czoo1,czoo2 parameter (cnh41=c1em3,cnh42=c0,cno31=p5,cno32=c0,cphy1=p5, * cphy2=c0,czoo1=0.05,czoo2=c0) #endif #ifdef bioAnder FLOAT * cdet1,cdet2,cnh41,cnh42,cno31,cno32,cphy1,cphy2,czoo1,czoo2 parameter (cdet1=c1em3,cdet2=c0,cnh41=c1em3,cnh42=c0,cno31=p5, * cno32=c0,cphy1=p5,cphy2=c0,czoo1=0.05,czoo2=c0) #endif #ifdef bioDuse FLOAT * cdet1,cdet2,cnh41,cnh42,cno31,cno32,czoo1,czoo2,cqn31, * cqn32,cqn41,cqn42,cchl1,cchl2 parameter (cdet1=c0,cdet2=c0,cnh41=c0,cnh42=c0,cno31=0.75, * cno32=c0,czoo1=c0,czoo2=c0,cqn31=0.018,cqn32=c0, * cqn41=c0,cqn42=c0,cchl1=0.04,cchl2=c0) FLOAT * wsnkphy2 #endif #ifdef bioFasham FLOAT * cbac1,cbac2,cdon1,cdon2,cnh41,cnh42,cno31,cno32,cphy1,cphy2, * cpon1,cpon2,czoo1,czoo2 parameter (cbac1=c1em3,cbac2=c0,cdon1=c1em3,cdon2=c0,cnh41=c1em3, * cnh42=c0,cno31=p5,cno32=c0,cphy1=p5,cphy2=c0, * cpon1=c1em3,cpon2=c0,czoo1=0.05,czoo2=c0) #endif FLOAT * rho(imt,km),sal(imt,km),tem(imt,km) parameter (delt=0.6799,a0=c0,a1=1.0e5,b0=delt,b1=1.0e5, * ssur=35.0,tsur=18.85) save first,sbfscl,sqfscl,swfscl data first /.true./ c c----------------------------------------------------------------------- c Begin excutable code. c----------------------------------------------------------------------- c c======================================================================= c Set domain information. c======================================================================= c entry anadinfo #ifdef grids c c Domain configuration and bottom topography is read from GRIDS NetCDF c file. c call readgrids #else c c Domain information is initialized in the named common blocks via c a BLOCK DATA sub-program (See blkdat.F). c c Compute flat level thicknesses HZ from them the flat level depths. c hz(1)=c2*refz(1) # ifndef barotropic do 20 k=2,km hz(k)=c2*(refz(k)-refz(k-1))-hz(k-1) if(hz(k).le.c0) then write(stdout,10) hz(k) 10 format(/,' ANAFLDS: too small thickness, HZ(k) = ',f12.2) call exitus('ANAFLDS') endif 20 continue # endif #endif #ifdef dblsigma c write(stdout,30) zc1,zc2,zslope,zref 30 format(/'Double-Sigma Parameters:'// * f10.2,' ZC1 minimum depth of the coordinate interface ', * '(cm)',/, * f10.2,' ZC2 maximum depth of the coordinate interface ', * '(cm)',/, * f10.4,' ZSLOPE slope parameter of the coordinate interface',/, * f10.2,' ZREF reference depth for the coordinate interface ', * '(cm)'/) #endif c return #ifdef coast c c======================================================================= c Setup coastal segments and land/sea mask at tracer points. c======================================================================= c entry anacoast # ifdef grids c c Read in Land/Sea mask data from GRIDS NetCDF file. c call landsea(ncinpid) # endif return #endif c c======================================================================= c Analytical Bottom topography (cm) at the tracer points. c======================================================================= c entry anatopo #ifndef grids c c Set analytical bottom topography at the tracer points (centimeters, c positive). c do 100 jj=1,jmt do 100 i=1,imt hr(i,jj)=10000.0 ??? 100 continue # ifdef cyclic c c User needs to insure Cyclic and continuous conditions on analytical c topography. c do 110 jj=1,jmt hr(1 ,jj)=hr(imtm1,jj) hr(imt,jj)=hr(2 ,jj) 110 continue # endif #endif return c c======================================================================= c Analytical initial conditions. c======================================================================= c entry anauvpt(j) c c Internal mode velocity comnponents (cm/s) c call depthslab(j,tgrid,tdepth(1,1,0)) call depthslab(j,vgrid,vdepth(1,1,0)) do 210 k=1,km do 200 i=1,imt depth=vdepth(i,k,0) u(i,k)=c0*depth ??? v(i,k)=c0*depth ??? 200 continue 210 continue c c Transport streamfunction (cm^3/s) c do 220 i=1,imt p(i,j)=c0 ??? 220 continue c c Tracers: temperature (deg C) and salinity (PSS) - SMEAN (PSS) c (Set an exponential temperature profile and salinity constant). c Add additional tracers as necessary. c do 240 k=1,km do 230 i=1,imt depth=tdepth(i,k,0) t(i,k,1)=tsur+a0*exp(-depth/a1) ??? t(i,k,2)=ssur-smean ??? #ifdef bioMcGillic t(i,k,ino3)=cno31+depth*cno32 ??? t(i,k,iphy)=cphy1+depth*cphy2 ??? t(i,k,izoo)=czoo1+depth*czoo2 ??? t(i,k,inh4)=cnh41+depth*cnh42 ??? #endif #ifdef bioAnder t(i,k,ino3)=cno31+depth*cno32 ??? t(i,k,iphy)=cphy1+depth*cphy2 ??? t(i,k,izoo)=czoo1+depth*czoo2 ??? t(i,k,inh4)=cnh41+depth*cnh42 ??? t(i,k,idet)=cdet1+depth*cdet2 ??? #endif #ifdef bioDuse c if(depth.lt.5000) then c t(i,k,ino3)=c0 c else t(i,k,ino3)=cno31+depth*cno32 ??? c endif t(i,k,izoo)=czoo1+depth*czoo2 ??? t(i,k,inh4)=cnh41+depth*cnh42 ??? t(i,k,idet)=cdet1+depth*cdet2 ??? t(i,k,iqn3)=cqn31+depth*cqn32 ??? t(i,k,iqn4)=cqn41+depth*cqn42 ??? t(i,k,ichl)=cchl1+depth*cchl2 ??? #endif #ifdef bioFasham t(i,k,ino3)=cno31+depth*cno32 ??? t(i,k,iphy)=cphy1+depth*cphy2 ??? t(i,k,izoo)=czoo1+depth*czoo2 ??? t(i,k,inh4)=cnh41+depth*cnh42 ??? t(i,k,ipon)=cpon1+depth*cpon2 ??? t(i,k,idon)=cdon1+depth*cdon2 ??? t(i,k,ibac)=cbac1+depth*cbac2 ??? #endif 230 continue 240 continue #ifdef cyclic c c Set Cyclic boundary conditions. c do 250 k=1,km u(1 ,k)=u(imtm1,k) u(imt,k)=u(2 ,k) v(1 ,k)=v(imtm1,k) v(imt,k)=v(2 ,k) do 250 m=1,nt t(1 ,k,m)=t(imtm1,k,m) t(imt,k,m)=t(2 ,k,m) 250 continue p(1 ,j)=p(imtm1,j) p(imt,j)=p(2 ,j) #endif return c c======================================================================= c Analytical background density. c======================================================================= c entry anamrho c do 350 jj=1,jmt call depthslab(jj,tgrid,tdepth(1,1,0)) do 320 k=1,km do 310 i=1,imt depth=tdepth(i,k,0) tem(i,k)=tsur+b0*exp(-depth/b1) ??? sal(i,k)=ssur-smean ??? 310 continue 320 continue call state(tem,sal,tdepth(1,1,0),rho) do 340 k=1,km do 330 i=1,imt rhobar(i,jj,k)=rho(i,k) 330 continue 340 continue 350 continue return c c======================================================================= c Analytical lateral boundary conditions. c======================================================================= c entry anabndry c do 440 ibdy=1,4 do 410 i=1,imx(ibdy) c c Transport streamfunction (cm^3/s) and time rate of change of c barotropic vorticity (1/s^2). c po(i,1,ibdy)=c0 ??? po(i,2,ibdy)=c0 ??? c c Internal mode velocity components (cm/s) c do 400 k=1,km depth=bvdepth(i,k,ibdy) uo(i,k,ibdy)=c0 ??? vo(i,k,ibdy)=c0 ??? 400 continue 410 continue c c Tracers. Add additional tracers as necessary. c do 430 i=1,imx(ibdy) do 420 k=1,km depth=btdepth(i,k,ibdy) to(i,k,ibdy,1)=tsur+a0*exp(-depth/a1) ??? to(i,k,ibdy,2)=ssur-smean ??? #ifdef bioMcGillic to(i,k,ibdy,ino3)=cno31+depth*cno32 ??? to(i,k,ibdy,iphy)=cphy1+depth*cphy2 ??? to(i,k,ibdy,izoo)=czoo1+depth*czoo2 ??? to(i,k,ibdy,inh4)=cnh41+depth*cnh42 ??? #endif #ifdef bioAnder to(i,k,ibdy,ino3)=cno31+depth*cno32 ??? to(i,k,ibdy,iphy)=cphy1+depth*cphy2 ??? to(i,k,ibdy,izoo)=czoo1+depth*czoo2 ??? to(i,k,ibdy,inh4)=cnh41+depth*cnh42 ??? to(i,k,ibdy,idet)=cdet1+depth*cdet2 ??? #endif #ifdef bioDuse to(i,k,ibdy,ino3)=cno31+depth*cno32 ??? to(i,k,ibdy,izoo)=czoo1+depth*czoo2 ??? to(i,k,ibdy,inh4)=cnh41+depth*cnh42 ??? to(i,k,ibdy,idet)=cdet1+depth*cdet2 ??? to(i,k,ibdy,iqn3)=cqn31+depth*cqn32 ??? to(i,k,ibdy,iqn4)=cqn41+depth*cqn42 ??? to(i,k,ibdy,ichl)=cchl1+depth*cchl2 ??? #endif #ifdef bioFasham to(i,k,ibdy,ino3)=cno31+depth*cno32 ??? to(i,k,ibdy,iphy)=cphy1+depth*cphy2 ??? to(i,k,ibdy,izoo)=czoo1+depth*czoo2 ??? to(i,k,ibdy,inh4)=cnh41+depth*cnh42 ??? to(i,k,ibdy,ipon)=cpon1+depth*cpon2 ??? to(i,k,ibdy,idon)=cdon1+depth*cdon2 ??? to(i,k,ibdy,ibac)=cbac1+depth*cbac2 ??? #endif 420 continue 430 continue 440 continue return c c======================================================================= c Set the depths at boundary points. c======================================================================= c entry dpthbndry c do 550 jj=1,jmt call depthslab(jj,tgrid,tdepth(1,1,0)) call depthslab(jj,vgrid,vdepth(1,1,0)) c c Depths at southern boundary. c if(jj.eq.1) then do 510 k=1,km do 500 i=1,imt bvdepth(i,k,south)=vdepth(i,k,0) btdepth(i,k,south)=tdepth(i,k,0) 500 continue 510 continue c c Depths at northern boundary. c elseif(jj.eq.jmtm1) then do 530 k=1,km do 520 i=1,imt bvdepth(i,k,north)=vdepth(i,k,0) btdepth(i,k,north)=tdepth(i,k,0) 520 continue 530 continue c c Depths at western and eastern boundaries. c else do 540 k=1,km bvdepth(jj,k,west)=vdepth(1,k,0) btdepth(jj,k,west)=tdepth(1,k,0) bvdepth(jj,k,east)=vdepth(imt,k,0) btdepth(jj,k,east)=tdepth(imt,k,0) 540 continue endif 550 continue return c c======================================================================= c Set vertical boundary conditions. c======================================================================= c entry anavbc(j) c c On first call, compute scale factors for fluxes. A scale factor from c W/m^2 to degC cm/s is provided as an example. User needs to select c the appropriate scale for the analytical forcing. c if(first) then sqfscl=c1e6/(rho0*cp) swfscl=sec2day sbfscl=sec2day first=.false. endif c c Set surface momentum flux: wind stress (dynes/cm^2). c #ifdef forcing do 600 i=1,imt smf(i,1)=FLoaT(j)*c0 ??? smf(i,2)=FLoaT(j+1)*c0 ??? 600 continue #else do 600 i=1,imt smf(i,1)=c0 smf(i,2)=c0 600 continue #endif c c Set bottom momentum flux: bottom stress. c do 610 i=1,imt kz=kmu(i) if(kz.ne.0) then uvmag=sqrt(ub(i,kz)**2+vb(i,kz)**2) bmf(i,1)=cdbot*ub(i,kz)*uvmag bmf(i,2)=cdbot*vb(i,kz)*uvmag else bmf(i,1)=c0 bmf(i,2)=c0 endif 610 continue c c Set surface (net) heat flux at rows J and J+1 (convert from W/m^2 c to degC cm/seg). c #ifdef forcing do 620 i=1,imt stf(i,1)=sqfscl*FLoaT(j)*c0 ??? stfp(i,1)=sqfscl*FLoaT(j+1)*c0 ??? 620 continue #else do 620 i=1,imt stf(i,1)=c0 stfp(i,1)=c0 620 continue #endif c c Set surface water flux at rows J and J+1: E-P (convert from cm/day c to ppt cm/s. c #ifdef forcing do 630 i=1,imt fx=swfscl*FLoaT(j)*c0 ??? stf(i,2)=(tb(i,1,2)+smean)*fx fx=swfscl*FLoaT(j+1)*c0 ??? stfp(i,2)=(tbp(i,1,2)+smean)*fx 630 continue #else do 630 i=1,imt stf(i,2)=c0 stfp(i,2)=c0 630 continue #endif #if defined bioAnder | defined bioFasham | defined bioMcGillic | defined bioDuse c c Set surface flux for biological tracers. Convert to appropriate c units. c # ifdef forcing c c Set dilution flux due to surface water flux, if any. c do 640 m=3,nt do 640 i=1,imt fx=stf(i,2)/(tb(i,1,2)+smean) stf(i,m)=tb(i,1,m)*fx 640 continue c c Surface nitrate flux. c # ifndef bioDuse do 650 i=1,imt fx=sbfscl*FLoaT(j)*c0 ??? stf(i,ino3)=tb(i,1,ino3)*fx ??? 650 continue # endif # else do 650 m=3,nt do 650 i=1,imt stf(i,m)=c0 650 continue # endif #endif c c Set bottom tracer flux. c do 660 m=1,nt do 660 i=1,imt btf(i,m)=c0 #ifdef bioAnder if(m.eq.iphy) btf(i,m)=wsnkphy*tb(i,km,iphy) if(m.eq.idet) btf(i,m)=wsnkdet*tb(i,km,idet) if(m.eq.ino3) btf(i,m)=-fracrmn*(wsnkphy*tb(i,km,iphy)+ * wsnkdet*tb(i,km,idet)) #endif #ifdef bioDuse if(wsnkphy.lt.c0) then if(tb(i,km,ichl).gt.c0) then wsnkphy2=c2*(chl2rad*tb(i,km,ichl)**chl2rade)* * grav*cm2m* * (rhophy-rhos(i,km) # ifdef rmdenbar * -rhobar(i,j,km) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if(biopos.eq.2) then if(m.eq.iqn3) btf(i,m)=wsnkphy2*max(c0,tb(i,km,iqn3)) if(m.eq.iqn4) btf(i,m)=wsnkphy2*max(c0,tb(i,km,iqn4)) if(m.eq.ichl) btf(i,m)=wsnkphy2*max(c0,tb(i,km,ichl)) if(m.eq.idet) btf(i,m)=wsnkdet*max(c0,tb(i,km,idet)) if(m.eq.ino3) btf(i,m)=-fracrmn*(wsnkphy2* * (max(c0,tb(i,km,iqn3)) * +max(c0,tb(i,km,iqn4))) * +wsnkdet*max(c0,tb(i,km,idet))) else if(m.eq.iqn3) btf(i,m)=wsnkphy2*tb(i,km,iqn3) if(m.eq.iqn4) btf(i,m)=wsnkphy2*tb(i,km,iqn4) if(m.eq.ichl) btf(i,m)=wsnkphy2*tb(i,km,ichl) if(m.eq.idet) btf(i,m)=wsnkdet*tb(i,km,idet) if(m.eq.ino3) btf(i,m)=-fracrmn*(wsnkphy2*( * tb(i,km,iqn3) * +tb(i,km,iqn4)) * +wsnkdet*tb(i,km,idet)) endif #endif #ifdef bioFasham if(m.eq.ipon) btf(i,m)=wsnkpon*tb(i,km,ipon) if(m.eq.inh4) btf(i,m)=-wsnkpon*tb(i,km,ipon) #endif 660 continue c c Set shortwave radiation flux at rows J and J+1 (convert from W/m^2 c to degC cm/sec). c #ifdef forcing do 670 i=1,imt srf(i)=sqfscl*217.0 ??? srfp(i)=sqfscl*217.0 ??? 670 continue #else do 670 i=1,imt srf(i)=c0 srfp(i)=c0 670 continue #endif c c======================================================================= c Compute, monitor, and report other User depended diagnostics. c======================================================================= c entry anadiag(j) return end