module physics
use const
use comm
use var
use io
implicit real*8 (a-h,o-z)
contains
! ----------------------------------------------------------------------
subroutine windglo(u,v,days,dt,odz)
integer nld
dimension u(i0,j0),v(i0,j0),odz(*)
common/klimatglo/fnew,fold,new,nld
common/windsglo/taux(i2,j2,12),tauy(i2,j2,12)
! the velocity arrays u,v have a perimeter "ghost zone".
! thus, we set wind forcing only in the interior zones.
! taux,tauy are surface wind stress components.
! taux,tauy units are force per unit area (i.e., energy per unit volume).
! tmp=odz(1)/rho
! all units are cgs, so we use rho=1.
if(myid.eq.0) write(*,*) 'days = ',days
tmp=dt*odz(1)
do j=2,j1
do i=2,i1
u(i,j)=u(i,j)+tmp*(fold*taux(i-1,j-1,nld)+fnew*taux(i-1,j-1,new))
v(i,j)=v(i,j)+tmp*(fold*tauy(i-1,j-1,nld)+fnew*tauy(i-1,j-1,new))
enddo
enddo
end subroutine
! ----------------------------------------------------------------------
subroutine cfcsurf(c,t,s,days,dt,odz)
dimension c(i0,j0),t(i0,j0),s(i0,j0),odz(*)
!local
logical,save :: first_call=.true.
integer, parameter :: cfc_nt = 2
integer, parameter :: nyrbeg = 31
integer, parameter :: nyrend = 99
integer,parameter :: imt=i2
integer,parameter :: jmt=j2
real*8,save :: year(nyrend),p11(nyrend,cfc_nt),p12(nyrend,cfc_nt)
real*8,save :: fice(i2,j2,12),xkw(i2,j2,12),patm(i2,j2,12)
real*8,save :: ratio,day0,nday
real*8,save :: p11_new(cfc_nt),p12_new(cfc_nt)
real*8,save :: pnorth(cfc_nt),psouth(cfc_nt)
real*8 :: rlon(imt,jmt),rlat(imt,jmt),CFCatm(imt,jmt,cfc_nt)
real*8 :: Kw,CFCsat
real*8 :: cfc_sc,cfc_alpha
real*8, save :: CMAX,CMIN,STMAX,STMIN,SSMAX,SSMIN
real*8, save :: KwMAX,scMAX,CFCMAX,alphaMAX,xkwMAX
real*8, save :: KwMIN,scMIN,CFCMIN,alphaMIN,xkwMIN
real*8 :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9
real*8 :: tmp10,tmp11,tmp12,tmp13,tmp14,tmp15,tmp16
INTEGER :: I,J,K,kn
integer,save :: indx1,indx2,month
!local
if(first_call) then
! if(myid.eq.0)write(*,*)'days,dt=',days,dt
if(myid.eq.0)write(*,*)'first call Read CFC gas and atm!'
if(myid.eq.0)write(*,*)'i0,i1,i2=',i0,i1,i2
if(myid.eq.0)write(*,*)'j0,j1,j2=',j0,j1,j2
call read_gas(fice,xkw,patm)
call read_cfcatm(cfc_nt,nyrbeg,nyrend,year,p11,p12)
!Initial year in file
indx1 = 90 !90 !31 !45 !1 !90 !we can use 29 =>spin up 2 yrs and run 37yrs
indx2 = indx1+1
month=1
day0 = 0. !10800.
if(myid.eq.0)write(*,*)'initial year,month:',1900+indx1,month
STMAX = -9999.0
SSMAX = -9999.0
CMAX = -9999.0
KwMAX = -9999.0
scMAX = -9999.0
CFCMAX = -9999.0
alphaMAX = -9999.0
xkwMAX = -9999.0
STMIN = 9999.0
SSMIN = 9999.0
CMIN = 9999.0
KwMIN = 9999.0
scMIN = 9999.0
CFCMIN = 9999.0
alphaMIN = 9999.0
xkwMIN = 9999.0
do j=1,j2
do i=1,i2
rlon(i,j) = xdeg(i+1)
rlat(i,j) = ydeg(j+1)
! if(myid.eq.0) write(99,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.1) write(98,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.2) write(97,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.3) write(96,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.4) write(95,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.5) write(94,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.6) write(93,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.7) write(92,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.8) write(91,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.9) write(90,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.10) write(89,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.11) write(88,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.12) write(87,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.13) write(86,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.14) write(85,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
! if(myid.eq.15) write(84,*)rlon(i,j),rlat(i,j),xkw(i,j,12)
enddo
enddo
first_call=.false.
endif !first call
nday = days-day0;
if (mod(days,30.).eq.0.) then
month = month +1
month = min(month,12)
if(myid.eq.0)write(*,*)'month=',month
end if
if (mod(days,360.).eq.0.) then
month = 1
day0 = days
indx1 = indx2
indx2 = indx1+1
if(myid.eq.0)write(*,*)'New year: month,year1,year2=',month,indx1,indx2
end if
ratio = nday/360. !360.!360.0
p11_new(1) = p11(indx1,1)*(1-ratio)+p11(indx2,1)*ratio
p11_new(2) = p11(indx1,2)*(1-ratio)+p11(indx2,2)*ratio
p12_new(1) = p11(indx1,1)*(1-ratio)+p11(indx2,1)*ratio
p12_new(2) = p11(indx1,2)*(1-ratio)+p11(indx2,2)*ratio
pnorth(1)=p11_new(1)
pnorth(2)=p12_new(1)
psouth(1)=p11_new(2)
psouth(2)=p12_new(2)
call cfc_interp(cfc_nt,pnorth, psouth, rlat, imt, jmt, CFCatm)
! if(mod(days,1.).eq.0) then
! do j=1,J2
! do i=1,I2
! if(myid.eq.0)write(99,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1)
! if(myid.eq.1)write(98,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1)
! if(myid.eq.2)write(97,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1)
! if(myid.eq.3)write(96,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1)
! end do
! end do
! endif
kn = 11 !11: CFC-11 12:CFC-12
! if(myid.eq.0) write(*,*)'T,S=',T(70,60),S(70,60)
! if(myid.eq.0) call sc_cfc(T(70,60),kn,cfc_sc)
! if(myid.eq.0) call sol_cfc(T(70,60),S(70,60),kn,cfc_alpha)
! if(myid.eq.0)write(*,*)'alpha=',cfc_sc,cfc_alpha
do j=1,j2
do i=1,i2
call sc_cfc(t(I+1,J+1),kn,cfc_sc)
call sol_cfc(t(I+1,J+1),s(I+1,J+1),kn,cfc_alpha)
CFCsat = cfc_alpha*CFCatm(I,J,1) !1=CFC-11
! convert to pmol/kg for comparison with GLODAP data set
CFCsat = CFCsat*1e12/1030. !pmol/kg
Kw = (1-fice(I,J,month))*xkw(I,J,month)*(cfc_sc/660.0)**(-0.5) !m/s
Kw = Kw*ODZ(1)*100. !the latter *100 because ODZ is in 1/cm rather than 1/m
Kw = Kw*dt
! if(myid.eq.0) then
! if(i.eq.50.and.j.eq.50) then
! write(*,*)'check kw',myid,month,xkw(I+1,J+1,month)
! endif
! endif
c(I+1,J+1) = c(I+1,J+1) + Kw*(CFCsat-c(I+1,J+1))*IN(I+1,J+1,1)
STMAX=MAX(STMAX,T(I+1,J+1))
SSMAX=MAX(SSMAX,S(I+1,J+1))
CMAX=MAX(CMAX,C(I+1,J+1))
KwMAX=MAX(KwMAX,Kw)
scMAX=MAX(scMAX,cfc_sc)
CFCMAX=MAX(CFCMAX,CFCsat)
alphaMAX=MAX(alphaMAX,cfc_alpha)
xkwMAX=MAX(xkwMAX,xkw(I,J,month))
STMIN=MIN(STMIN,T(I+1,J+1))
SSMIN=MIN(SSMIN,S(I+1,J+1))
CMIN=MIN(CMIN,C(I+1,J+1))
KwMIN=MIN(KwMIN,Kw)
scMIN=MIN(scMIN,cfc_sc)
CFCMIN=MIN(CFCMIN,CFCsat)
alphaMIN=MIN(alphaMIN,cfc_alpha)
xkwMIN=MIN(xkwMIN,xkw(I,J,month))
end do
end do
tmp1=cmax
tmp2=cmin
tmp3=cfcmax
tmp4=cfcmin
tmp5=kwmax
tmp6=kwmin
tmp7=scmax
tmp8=scmin
tmp9=alphamax
tmp10=alphamin
tmp11=stmax
tmp12=stmin
tmp13=ssmax
tmp14=ssmin
tmp15=xkwmax
tmp16=xkwmin
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp1,cmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp2,cmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp3,cfcmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp4,cfcmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp5,kwmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp6,kwmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp7,scmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp8,scmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp9,alphamax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp10,alphamin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp11,stmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp12,stmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp13,ssmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp14,ssmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp15,xkwmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR)
CALL MPI_BARRIER(M_CART,IERR)
CALL MPI_ALLREDUCE(tmp16,xkwmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR)
if(myid.eq.0) then
if (mod(days,1.0).eq.0.) then
write(*,*)'CFC check! DT=',DT
write(*,*)days,month,nday,day0,ratio,indx1,indx2
write(*,*)stmax,ssmax,cmax,cfcmax,kwmax,scmax,alphamax
write(*,*)stmin,ssmin,cmin,cfcmin,kwmin,scmin,alphamin
write(*,*)xkwmax,xkwmin,odz(1),odz(1)*100.0*dt
end if
endif
end subroutine
subroutine sc_cfc(t,kn,cfc_sc)
!---------------------------------------------------
! CFC 11 and 12 schmidt number
! as a fonction of temperature.
!
! ref: Zheng et al (1998), JGR, vol 103,No C1
!
! t: temperature (degree Celcius)
! kn: = 11 for CFC-11, 12 for CFC-12
!
! J-C Dutay - LSCE
!---------------------------------------------------
! use cfc_glbl, only: rkind
implicit none
integer, intent(in) :: kn
real*8 :: a1( 11: 12), a2( 11: 12), a3( 11: 12), a4( 11: 12)
real*8,intent(in) :: t
real*8,intent(out) :: cfc_sc
!
! coefficients with t in degre Celcius
! ------------------------------------
a1(11) = 3501.8
a2(11) = -210.31
a3(11) = 6.1851
a4(11) = -0.07513
!
a1(12) = 3845.4
a2(12) = -228.95
a3(12) = 6.1908
a4(12) = -0.067430
!
cfc_sc = a1(kn)+a2(kn)*t+a3(kn)*t*t+a4(kn)*t*t*t
end subroutine
subroutine sol_cfc(pt,ps,kn,cfc_alpha)
implicit none
real*8,intent(in) :: pt, ps
real*8,intent(out) :: cfc_alpha
real*8 :: ta,d
real*8 :: a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12)
real*8 :: b1 ( 11: 12), b2 ( 11: 12), b3 ( 11: 12)
integer, intent(in) :: kn
!!
!! coefficient for solubility in mol/l/atm
!! ----------------------------------------
!
! for CFC 11
! ----------
a1 ( 11) = -229.9261
a2 ( 11) = 319.6552
a3 ( 11) = 119.4471
a4 ( 11) = -1.39165
b1 ( 11) = -0.142382
b2 ( 11) = 0.091459
b3 ( 11) = -0.0157274
!
! for CFC/12
! ----------
a1 ( 12) = -218.0971
a2 ( 12) = 298.9702
a3 ( 12) = 113.8049
a4 ( 12) = -1.39165
b1 ( 12) = -0.143566
b2 ( 12) = 0.091015
b3 ( 12) = -0.0153924
!
ta = ( pt + 273.16)* 0.01
d = ( b3 ( kn)* ta + b2 ( kn))* ta + b1 ( kn)
!
!
cfc_alpha = exp ( a1 ( kn) + a2 ( kn)/ ta + a3 ( kn)* log ( ta ) + &
a4 (kn)* ta * ta + ps* d )
!
! conversion from mol/(l * atm) to mol/(m^3 * atm)
! ------------------------------------------------
cfc_alpha = 1000. * cfc_alpha
!
! conversion from mol/(m^3 * atm) to mol/(m3 * pptv)
! --------------------------------------------------
cfc_alpha = 1.0e-12 * cfc_alpha
end subroutine
! ----------------------------------------------------------------------
subroutine cfc_interp(cfc_nt,pnorth,psouth,rlat,imt,jmt,CFCatm)
integer, parameter :: ijmax=200*200
integer, intent(in) :: cfc_nt, imt, jmt
integer,save :: ientry
integer :: i, j, ij, n
real*8,save :: ys, yn
real*8,save :: xintp(ijmax)
real*8, intent(in) :: Pnorth(cfc_nt), Psouth(cfc_nt)
real*8, intent(in) :: rlat(imt,jmt)
real*8, intent(out) :: CFCatm(imt,jmt,cfc_nt)
DATA ys, yn / -10., 10./
DATA ientry /0/
ientry = ientry + 1
! Test to see if ijmax is large enough
if (ijmax .lt. imt*jmt) then
print *," cfc_interp: ERROR -> ijmax must be at least",imt*jmt
stop
endif
! IF Block to be executed ONLY on first call to this routine
if (ientry .eq. 1) then
do j = 1,jmt
do i = 1,imt
ij = (j-1)*imt + i
if (rlat(i,j) .ge. yn) then
xintp(ij) = 1.
else if (rlat(i,j) .le. ys) then
xintp(ij) = 0.
else
xintp(ij) = (rlat(i,j) - ys)/(yn-ys)
endif
end do
end do
endif
! Block to be executed every pass
do n=1,cfc_nt
do j=1,jmt
do i=1,imt
ij = (j-1)*imt + i
CFCatm(i,j,n) = xintp(ij) * Pnorth(n) + &
(1.0 - xintp(ij)) * Psouth(n)
end do
end do
end do
end subroutine
! ----------------------------------------------------------------------
subroutine read_cfcatm(cfc_nt,nyrbeg,nyrend,year,p11,p12)
! needs modification to dbl
!local
integer :: i, n
integer, parameter :: iu = 10
integer,intent(in) :: nyrbeg,nyrend,cfc_nt
real*8,intent(out) :: year(nyrend), p11(nyrend,cfc_nt), p12(nyrend,cfc_nt)
character(len=80) :: filen
year = 0.
p11 = 0.
p12 = 0.
!
! OPEN FILE
! ---------
filen='data/cfc1112_new.atm'
open(UNIT=iu,FILE=filen,STATUS='old')
!
! Skip over 1st six descriptor lines
! ----------------------------------
do i=1,6
read(iu,99)
end do
!
! READ FILE
! ---------
do n = nyrbeg, nyrend
READ(iu,*)year(n), p11(n,1), p12(n,1), p11(n,2), p12(n,2)
year(n) = year(n) - 0.5
if(myid.eq.0) then
write(*,*) 'read cfcatm file', year(n),p11(n,1), p12(n,1), p11(n,2), p12(n,2)
endif
end do
close(iu)
99 FORMAT(1x)
100 FORMAT(f7.2, 4f8.2)
end subroutine
! ----------------------------------------------------------------------
subroutine pp82glo
! ----------------------------------------------------------------------
! vertical mixing coefficients
! approach:
! set background vertical diffusivities to dmz0 (o(0.1) cm-cm/sec) and
! add richardson number based vertical diffusivity plus numerical
! contribution according to vertical cell reynolds number.
! vertical cell re is o(100), except during winter cooling conditions
! and in wind-blown surface mixed layer, because internal waves, which
! dominate w below the sml in the model results during summer, do not
! mix t or s in nature.
! note: one cannot depend on diffusive closure by itself if one wants to
! model possible contra-diffusive effe!
parameter(i01=i0+1)
common/windmxglo/vbk(k2),hbk(k2)
common/tauglo/add(i2,j2,k2)
do k=1,k2
l=k+1
dzwl=orzmx/odzw(l)
! stability limit
emax=.4/(dt*odzw(l)**2)
do j=2,j1
do i=2,i1
! pacanowski and philander (1981): evisc=a*(1/(1+b*ri))**n
! evisc=turbulent eddy viscosity; diffuse=turbulent diffusivity
! ri=max(0.,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2*
ri=max(-0.02d0,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2*(0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2)))
rrinv=1./(1.+5.*ri)
tmp=rrinv**2
! ri- and cell re- dependent vertical mixing
! numerically kosher because, for large ri, laminar flow dominates
! and time mean overshoot effe! are limited by time mean w being small
temp=tmp*abs(w(i,j,l))*dzwl
evisc=5.*tmp
diffuse=evisc*rrinv
! add parameterizes extra momentum dissipation along steep slopes,
! high and low latitudes, and equatorial regions due to breaking of
! waves that were generated by big storm events that are missing when
! using climatological winds. add includes vbk(k) term.
ev(i-1,j-1,k)=evisc+add(i-1,j-1,k)+temp
hv(i-1,j-1,k)=diffuse+hbk(k)+temp
! apply stability limit and odzw factor
tmp=odzw(l)*iw(i,j,l)
ev(i-1,j-1,k)=tmp*min(emax,ev(i-1,j-1,k))
hv(i-1,j-1,k)=tmp*min(emax,hv(i-1,j-1,k))
enddo
enddo
enddo
end subroutine
! ----------------------------------------------------------------------
subroutine pp82glonew
! ----------------------------------------------------------------------
! vertical mixing coefficients
! approach:
! set background vertical diffusivities to dmz0 (o(0.1) cm-cm/sec) and
! add richardson number based vertical diffusivity plus numerical
! contribution according to vertical cell reynolds number.
! vertical cell re is o(100), except during winter cooling conditions
! and in wind-blown surface mixed layer, because internal waves, which
! dominate w below the sml in the model results during summer, do not
! mix t or s in nature.
! note: one cannot depend on diffusive closure by itself if one wants to
! model possible contra-diffusive effe!
parameter(i01=i0+1)
common/windmxglo/vbk(k2),hbk(k2)
rate=-.05
jt=mylat*j2
do k=1,k2
l=k+1
dzwl=orzmx/odzw(l)
! deepen and augment mixing in high lats to emulate big storms
tmpk=10.*exp(-.0001*z(2*k+1))
! stability limit
emax=.4/(dt*odzw(l)**2)
elim=min(200.,emax)
do j=2,j1
! decrease background mixing except in high lats to emulate big
! storms within 20 grid intervals (~2 deg lat) of lat boundaries
tmp=tmpk*(exp(rate*abs(2-j-jt)))+exp(rate*abs(j1t-j-jt))
vb=vbk(k)+tmp
hb=hbk(k)+tmp
do i=2,i1
! pacanowski and philander (1981): evisc=a*(1/(1+b*ri))**n
! evisc=turbulent eddy viscosity; diffuse=turbulent diffusivity
tmp=odzw(l)
! starting at day 0 of yr 16, decrease ri-based vertical mixing
! gap in izu ridge is also deepened to emulate unresolved deep kuroshio paths
! thus allowing it to pass over the ridgeline with less anomalous vortex
! squashing and thus taking a non-meander path when appropriate.
oldfact=980.*(rho(i,j,l)-rho(i,j,k))*tmp/(tmp**2*(0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2))
! big ri decreases vertical mixing
ri=max(-0.19,5.*oldfact)
rrinv=1./(1.+5.*ri)
tmp=rrinv**2
! ri- and cell re- dependent vertical mixing
! numerically kosher because, for large ri, laminar flow dominates
! and time mean overshoot effe! are limited by time mean w being small
temp=.01*tmp*abs(w(i,j,l))*dzwl
evisc=tmp
diffuse=evisc*rrinv
! evisc, diffuse, hb and vb units are cm-cm/sec
ev(i-1,j-1,k)=evisc+vb+temp
hv(i-1,j-1,k)=diffuse+hb+0.1*temp
! apply stability limit and odzw factor
tmp=odzw(l)*iw(i,j,l)
ev(i-1,j-1,k)=tmp*min(elim,ev(i-1,j-1,k))
hv(i-1,j-1,k)=tmp*min(elim,hv(i-1,j-1,k))
enddo
enddo
enddo
! ======================================================================
! add vertical mixing to emulate effect of equatorial vortices and
! penetrative convection in acc.
! ======================================================================
! stacked vortices in 10-point equator band (hardwired to 1/4 deg res)
! strong velocity oscillations having ~20-day time scale are observed
! within a few deg of equator
! 33-point equatorial and acc bands
!
! do 770 i=1,i2
! do 770 j=1,j2
! do 770 k=1,k2
! 770 scr(i+1,j+1,k)=ev(i,j,k)
sc1=.01
sc2=.01
do k=1,k2
l=k+1
! 100m scale height of added vertical mixing
tmp=exp(-.0001*z(2*k+1))
do j=1,j2
tmp1=10./(1.+sc1*((0-ydeg(j))/(ydeg(j+1)-ydeg(j)))**2)
addd=tmp1*odzw(l)*tmp
do i=1,i2
! insulated bottom, nonlinear drag invoked in fs
ad=addd*iw(i+1,j+1,l)
hv(i,j,k)=hv(i,j,k)+ad
ev(i,j,k)=ev(i,j,k)+ad
enddo
enddo
enddo
end subroutine
!--------------------------------------------------------------------------
subroutine kppglo
end subroutine kppglo
!--------------------------------------------------------------------------
! ----------------------------------------------------------------------
subroutine smagglo
! ----------------------------------------------------------------------
! horizontal mixing coefficients
! approach:smagorinsky, j. general circulation experiments with the
! primitive equations, i. the basi! experiment. monthly weather rev. 1963,
! 91, 99-164.
parameter(i01=i0+1)
do k=1,k1
do j=2,j1
temp=ocs(j)*ody(j)
do i=2,i1
scr(i,j,k)=0.12*dx(j)*dy(j)*in(i,j,k)*sqrt(((u(i,j,k)-u(i-1,j,k))*odx(j))**2+((csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)**2)
enddo
enddo
enddo
call mpi_sendrecv(scr(2,1,1),1,m_vlon,m_w,1,scr(i0,1,1),1,m_vlon,m_e,1,m_cart,istat,ierr)
call mpi_sendrecv(scr(i1,1,1),1,m_vlon,m_e,1,scr(1,1,1),1,m_vlon,m_w,1,m_cart,istat,ierr)
call mpi_sendrecv(scr(1,2,1),1,m_vlat,m_s,1,scr(1,j0,1),1,m_vlat,m_n,1,m_cart,istat,ierr)
call mpi_sendrecv(scr(1,j1,1),1,m_vlat,m_n,1,scr(1,1,1),1,m_vlat,m_s,1,m_cart,istat,ierr)
do k=1,k1
do j=1,j0
do i=1,i1
dmx(i,j,k)=dm0+.5*(scr(i,j,k)+scr(i+1,j,k))
enddo
enddo
do i=1,i0
dmy(i,1,k)=dm0+scr(i,2,k)
dmy(i,j1,k)=dm0+scr(i,j2,k)
do j=1,jf
! note: except in hi lats, dmy slightly violates ang. mom. conservation
dmy(i,j,k)=dm0+.5*(scr(i,j,k)+scr(i,j+1,k))
enddo
enddo
enddo
end subroutine
! ----------------------------------------------------------------------
subroutine qsurfglo(t,dt,odz)
dimension t(i0,j0),odz(*)
! the temperature array t has a permiter "ghost zone".
! thus, we set atmospheri! heat exchange only in the interior zones.
! typical latitudinal range: -50 to 50 watts per square meter.
! 1 watt per square meter = 1 erg per second per equare cm.
! example: surface heating qdot = 50 watts per square meter
qdot=50.
! tmp is conversion factor from watts per square meter
! to deg ! per time step in top model layer.
! tmp=dt*odz(1)/rho/cp
! all units are cgs, so we use rho=1.
! cp for water is 2.5e4 ergs/gram/deg c.
!check this cp value!!!
tmp=dt*odz(1)/2.5e4
dtemp=tmp*qdot
do j=1,j2
do i=1,i2
t(i+1,j+1)=t(i+1,j+1)+dtemp
enddo
enddo
end subroutine
! ----------------------------------------------------------------------
subroutine sunglo(t,dt,odz,qprof,kb)
integer*2 kb
dimension t(i0,j0,k1),odz(k1),qprof(k1),kb(i0,j0)
! this is like subroutine qsurf except
! we use predetermined input array qprof to vertically distribute qdot.
! qprof(k) is the cumulative fraction of total heat absorbed to depth
! of bottom of layer k for a pure water column.
! example: surface radiation flux of 10 watts per square meter
qdot=10.
fract=qprof(1)
do k=1,k1
dtemp=fract*qdot*dt/2.5e4
! odz(k) is not the correct factor for layers deeper than the local bottom.
! thus, odz factor must be added inside loop 99.
! this quickly programmed logi! must be checked.
do j=1,j2
do i=1,i2
kk=min(kb(i+1,j+1),k)
t(i+1,j+1,kk)=t(i+1,j+1,kk)+dtemp*odz(kk)
enddo
enddo
if (k.lt.k1) fract=qprof(k+1)-fract
enddo
end subroutine
end module