[2015年 12月 4日 金曜日 15:40:04 JST]
[~/2015.Ibnu.Indian.Ocean/ECCO.dr080g.Fortran/make_mld2]
[am@aofd165]
$ csh make_mld2.csh
#!/bin/csh -f
mkdir -vp dat
cat > a.f << EOF
program make_data
implicit none
!=======================================================================
!
!=======================================================================
integer, parameter :: nlon = 101
integer, parameter :: nlat = 126
integer, parameter :: nlev = 21
integer, parameter :: nvar = 11
real :: rlon(nlon)
real :: rlat(nlat)
real :: rlev(nlev)
real :: t(nlon,nlat,nlev)
real :: s(nlon,nlat,nlev)
real :: r(nlon,nlat,nlev)
real :: u(nlon,nlat,nlev)
real :: v(nlon,nlat,nlev)
real :: w(nlon,nlat,nlev)
real :: wet(nlon,nlat)
real :: wes(nlon,nlat)
real :: tm (nlon,nlat)
real :: tm1(nlon,nlat)
real :: tm2(nlon,nlat)
real :: sm (nlon,nlat)
real :: sm1(nlon,nlat)
real :: sm2(nlon,nlat)
real :: um (nlon,nlat)
real :: vm (nlon,nlat)
real :: buf(nlon,nlat,nlev)
real :: t1(nlon,nlat,nlev)
real :: t2(nlon,nlat,nlev)
real :: s1(nlon,nlat,nlev)
real :: s2(nlon,nlat,nlev)
real :: tent(nlon,nlat)
real :: sent(nlon,nlat)
real :: term1(nlon,nlat,nvar+3)
real :: term2(nlon,nlat,nvar+3)
real :: advx1(nlon,nlat)
real :: advy1(nlon,nlat)
real :: advz1(nlon,nlat)
real :: subx1(nlon,nlat)
real :: suby1(nlon,nlat)
real :: advx2(nlon,nlat)
real :: advy2(nlon,nlat)
real :: advz2(nlon,nlat)
real :: subx2(nlon,nlat)
real :: suby2(nlon,nlat)
integer :: imld(nlon,nlat)
integer :: imld1(nlon,nlat)
integer :: imld2(nlon,nlat)
real :: hm1(nlon,nlat)
real :: hm0(nlon,nlat)
real :: hm2(nlon,nlat)
real :: mld(nlon,nlat)
real :: ild(nlon,nlat)
real :: ivd(nlon,nlat)
real :: blt(nlon,nlat)
real :: d20(nlon,nlat)
real :: zw(0:nlev)
real :: zt(nlev)
real :: dz(nlev)
real :: fcor(nlat)
real :: xt(nlon,nlat)
real :: xu(nlon,nlat)
real :: xv(nlon,nlat)
real :: yt(nlon,nlat)
real :: yu(nlon,nlat)
real :: yv(nlon,nlat)
real :: dxt(nlon,nlat)
real :: dxu(nlon,nlat)
real :: dxv(nlon,nlat)
real :: dyt(nlon,nlat)
real :: dyu(nlon,nlat)
real :: dyv(nlon,nlat)
real :: umask(nlon,nlat,nlev)
real :: vmask(nlon,nlat,nlev)
real :: tmask(nlon,nlat,nlev)
real(8) :: t20 = 20.0
real(4) :: hm
real(4) :: dt, dt1, dt2, dt3
real(4) :: pi, deg2rad
real(4) :: radius
real :: mue, muw, mvn, mvs, mwt, mwb, div
real :: fue, fuw, fvn, fvs, fwt, fwb
real :: gue, guw, gvn, gvs, gwt, gwb
real :: hue, huw, hvn, hvs
real :: fume, fumw, fvmn, fvms
real :: fube, fubw, fvbn, fvbs
real :: hume, humw, hvmn, hvms
real :: hube, hubw, hvbn, hvbs
real :: gume, gumw, gvmn, gvms
real :: gube, gubw, gvbn, gvbs
real :: hme, hmw, hms, hmn
real :: hbe, hbw, hbs, hbn
real :: tmld, tsub
real :: smld, ssub
character(len=120) :: var(nvar) =
& (/'uadv ','vadv ','wadv ',
& 'xdif ','ydif ','zdif ','gmdif ','kppmix',
& 'kf_cor','sflux ','sdamp '/) !,'tice'/)
character(len=120) :: dir1,dir2
character(len=120) :: ifile, ifile1, ifile2, gfile
character(len=120) :: tfile
character(len=120) :: sfile
character(len=120) :: ofile
logical :: yes
real :: undef = -999.9
character(len=72) :: time_unit = 'days'
real(8) :: time
integer :: id(4), kt
integer :: iyr, mon, jyr
integer :: i, j, k, kb, kv
integer :: n
integer :: kle, klw, kln, kls, kp1, km1
integer :: irec, jrec
! dir = '/mnt/hgfs/Downloads/Ibnu/ECCO/'
dir1='../make_data/'
dir2='../make_budget/'
pi = 4.*atan(1.)
deg2rad = pi/180.0
radius = 6370.e3 ! [m]
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
do iyr = 1993, 2014
! do iyr = 2014 !1993, 1993
print *,'...',iyr
! ... open file
write(ifile,'(a,a,i4.4,a)')
& trim(dir1),'ecco-jpl.',iyr,'.nc'
inquire(file=trim(ifile),exist=yes)
if (.not. yes) then
print *,' error : no file : ',trim(ifile)
stop
end if
write(tfile,'(a,a,i4.4,a)')
& trim(dir2),'ecco-jpl-hbudget.',iyr,'.nc'
inquire(file=trim(tfile),exist=yes)
if (.not. yes) then
print *,' error : no file : ',trim(tfile)
stop
end if
write(sfile,'(a,a,i4.4,a)')
& trim(dir2),'ecco-jpl-sbudget.',iyr,'.nc'
inquire(file=trim(sfile),exist=yes)
if (.not. yes) then
print *,' error : no file : ',trim(sfile)
stop
end if
! ... read x,y,z
call read_nc_dim(ifile,'lon',nlon,rlon)
call read_nc_dim(ifile,'lat',nlat,rlat)
call read_nc_dim(ifile,'level',nlev,rlev)
zt = rlev
zw = 0
do k = 1, nlev-1
zw(k) = 0.5*(rlev(k)+rlev(k+1))
end do
do k = nlev, nlev
zw(k) = rlev(k)+(rlev(k)-zw(k-1))*2.0
end do
do k = 1, nlev
dz(k) = zw(k)-zw(k-1)
end do
! gfile = '/mnt/hgfs/Downloads/Ibnu/ECCO/ecco-jpl.grid.nc'
gfile = '../make_data/ecco-jpl.grid.nc'
call read_nc_data(gfile,'xt',nlon,nlat,xt(:,:),1,1,undef)
call read_nc_data(gfile,'yt',nlon,nlat,yt(:,:),1,1,undef)
call read_nc_data(gfile,'xv',nlon,nlat,xv(:,:),1,1,undef)
call read_nc_data(gfile,'yt',nlon,nlat,yt(:,:),1,1,undef)
call read_nc_data(gfile,'yu',nlon,nlat,yu(:,:),1,1,undef)
call read_nc_data(gfile,'yv',nlon,nlat,yv(:,:),1,1,undef)
do j = 1, nlat
do i = 1, nlon
dxt(i,j) = radius*1.0*deg2rad*cos(yt(i,j)*deg2rad)
dxu(i,j) = radius*1.0*deg2rad*cos(yu(i,j)*deg2rad)
dxv(i,j) = radius*1.0*deg2rad*cos(yv(i,j)*deg2rad)
end do
end do
do i = 1, nlon
do j = 1, nlat-1
dyt(i,j) = (yv(i,j+1)-yv(i,j))
dyu(i,j) = (yv(i,j+1)-yv(i,j))
end do
do j = nlat, nlat
dyt(i,j) = dyt(i,j-1)
dyu(i,j) = dyu(i,j-1)
end do
do j = 2, nlat
dyv(i,j) = (yt(i,j)-yt(i,j-1))
end do
do j = 1, 1
dyv(i,j) = dyv(i,j+1)
end do
do j = 1, nlat
dyt(i,j) = radius*dyt(i,j)*deg2rad
dyu(i,j) = radius*dyu(i,j)*deg2rad
dyv(i,j) = radius*dyv(i,j)*deg2rad
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
do mon = 1, 12
print *,'... ',mon
if (mon == 1) then
jyr = iyr-1
jrec = mon
if (iyr == 1993) then
jyr = iyr
jrec = 1
end if
else
jyr = iyr
jrec = mon-1
end if
write(ifile1,'(a,a,i4.4,a)')
& trim(dir1),'ecco-jpl.',jyr,'.nc'
inquire(file=trim(ifile1),exist=yes)
if (.not. yes) then
print *,' error #1: no file : ',trim(ifile1)
stop
end if
do k = 1, nlev
call read_nc_data(ifile1,'t',nlon,nlat,t1(:,:,k),k,jrec,undef)
call read_nc_data(ifile1,'s',nlon,nlat,s1(:,:,k),k,jrec,undef)
end do
call mixed_layer(nlon,nlat,nlev,zw,dz,t1,s1,r,
& imld1,hm1,ild,ivd,blt,tm1,sm1,undef)
if (mon == 12) then
jyr = iyr+1
jrec = 1
else
jyr = iyr
jrec = mon+1
end if
write(ifile2,'(a,a,i4.4,a)')
& trim(dir1),'ecco-jpl.',jyr,'.nc'
inquire(file=trim(ifile2),exist=yes)
if (.not. yes) then
print *,' error #2: no file : ',trim(ifile2)
stop
end if
do k = 1, nlev
call read_nc_data(ifile2,'t',nlon,nlat,t2(:,:,k),k,jrec,undef)
call read_nc_data(ifile2,'s',nlon,nlat,s2(:,:,k),k,jrec,undef)
end do
call mixed_layer(nlon,nlat,nlev,zw,dz,t2,s2,r,
& imld2,hm2,ild,ivd,blt,tm2,sm2,undef)
! ... read data
jrec = mon
umask = 1
vmask = 1
tmask = 1
do k = 1, nlev
call read_nc_data(ifile,'t',nlon,nlat,t(:,:,k),k,jrec,undef)
call read_nc_data(ifile,'s',nlon,nlat,s(:,:,k),k,jrec,undef)
call read_nc_data(ifile,'u',nlon,nlat,u(:,:,k),k,jrec,undef)
call read_nc_data(ifile,'v',nlon,nlat,v(:,:,k),k,jrec,undef)
call read_nc_data(ifile,'w',nlon,nlat,w(:,:,k),k,jrec,undef)
do j = 1, nlat
do i = 1, nlon
if (u(i,j,k) <= undef+1) umask(i,j,k) = 0.0
if (v(i,j,k) <= undef+1) vmask(i,j,k) = 0.0
if (t(i,j,k) <= undef+1) tmask(i,j,k) = 0.0
if (u(i,j,k) <= undef+1) u(i,j,k) = 0.0
if (v(i,j,k) <= undef+1) v(i,j,k) = 0.0
if (w(i,j,k) <= undef+1) v(i,j,k) = 0.0
end do
end do
end do
call mixed_layer(nlon,nlat,nlev,zw,dz,t,s,r,
& imld,mld,ild,ivd,blt,tm,sm,undef)
hm0 = mld
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
d20 = undef
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
kb = 0
do k = 1, nlev
if (t(i,j,k) <= undef+1) exit
if (t(i,j,k) >= t20) cycle
kb = k
exit
end do
if (kb <= 1) then
d20(i,j) = undef
else
d20(i,j) = zt(kb)
& + (t20-t(i,j,kb))/(t(i,j,kb-1)-t(i,j,kb))
& * (zt(kb)-zt(kb-1))
end if
end do
end do
um = undef
vm = undef
do j = 1, nlat
do i = 1, nlon
if (imld(i,j) == 0) cycle
um(i,j) = 0
vm(i,j) = 0
hm = 0
do k = 1, imld(i,j)
um(i,j) = um(i,j) + u(i,j,k)*dz(k)
vm(i,j) = vm(i,j) + v(i,j,k)*dz(k)
hm = hm + dz(k)
end do
um(i,j) = um(i,j)/hm
vm(i,j) = vm(i,j)/hm
end do
end do
!---------------------------------------------------------------------
! heat and salt budget
!---------------------------------------------------------------------
if (mon == 12) then
if (mod(iyr,4) == 0) then
dt2 = (30.+36.)*86400
dt = (36.)*86400
else
dt2 = (30.+35.)*86400
dt = (35.)*86400
end if
else
dt2 = (30.+30.)*86400
dt = (30.)*86400
end if
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
advx1 = undef
advy1 = undef
advz1 = undef
subx1 = undef
suby1 = undef
advx2 = undef
advy2 = undef
advz2 = undef
subx2 = undef
suby2 = undef
do j = 2, nlat-1
do i = 2, nlon-1
if (t(i,j,1) <= undef) cycle
kle = imld(i+1,j)
klw = imld(i-1,j)
kln = imld(i,j+1)
kls = imld(i,j-1)
hme = 0.0 ; hmw = 0.0; hmn = 0.0; hms = 0.0
hbe = 0.0 ; hbw = 0.0; hbn = 0.0; hbs = 0.0
fume = 0.0 ; fumw = 0.0; fvmn = 0.0; fvms = 0.0
fube = 0.0 ; fubw = 0.0; fvbn = 0.0; fvbs = 0.0
hm = 0
advx1(i,j) = 0
advy1(i,j) = 0
advz1(i,j) = 0
subx1(i,j) = 0
suby1(i,j) = 0
advx2(i,j) = 0
advy2(i,j) = 0
advz2(i,j) = 0
subx2(i,j) = 0
suby2(i,j) = 0
do k = 1, imld(i,j)
if (k == nlev) cycle
kp1 = min(k+1,nlev)
km1 = max(k-1,1)
mue = u(i+1,j,k)*dyu(i+1,j)*dz(k)*umask(i+1,j,k)
muw = u(i ,j,k)*dyu(i ,j)*dz(k)*umask(i,j,k)
mvn = v(i,j+1,k)*dxv(i,j+1)*dz(k)*vmask(i,j+1,k)
mvs = v(i,j ,k)*dxv(i,j )*dz(k)*vmask(i,j,k)
mwt = w(i,j,k )*dxt(i,j)*dyt(i,j)
mwb = w(i,j,kp1)*dxt(i,j)*dyt(i,j)
! ... u(T)x = (uT)x - T(u)x
! dut_dx = fue-fuw
! tdu_dx = t(i,j,k)*(mue-muw)
! udt_dx = dut_dx - tdu_dx
! advx1(i,j) = advx1(i,j) - udt_dx
hue = 0.5*(mue+abs(mue))*hm0(i,j)
& + 0.5*(mue-abs(mue))*hm0(i+1,j)
huw = 0.5*(muw+abs(muw))*hm0(i-1,j)
& + 0.5*(muw-abs(muw))*hm0(i,j)
hue = hue - mue*hm0(i,j)
huw = huw - muw*hm0(i,j)
fue = 0.5*(mue+abs(mue))*t(i,j,k)
& + 0.5*(mue-abs(mue))*t(i+1,j,k)
fuw = 0.5*(muw+abs(muw))*t(i-1,j,k)
& + 0.5*(muw-abs(muw))*t(i,j,k)
fue = fue - mue*t(i,j,k)
fuw = fuw - muw*t(i,j,k)
gue = 0.5*(mue+abs(mue))*s(i,j,k)
& + 0.5*(mue-abs(mue))*s(i+1,j,k)
guw = 0.5*(muw+abs(muw))*s(i-1,j,k)
& + 0.5*(muw-abs(muw))*s(i,j,k)
gue = gue - mue*s(i,j,k)
guw = guw - muw*s(i,j,k)
if (k <= kle) then
hme = hme + dz(k)
hume = hume + hue
fume = fume + fue
gume = gume + gue
else
hbe = hbe + dz(k)
hube = hube + hue
fube = fube + fue
gube = gube + gue
end if
if (k <= klw) then
hmw = hmw + dz(k)
humw = humw + huw
fumw = fumw + fuw
gumw = gumw + guw
else
hbw = hbw + dz(k)
hubw = hubw + huw
fubw = fubw + fuw
gubw = gubw + guw
end if
! ... v(T)y = (vT)y - T(v)y
! dvt_dy = fvn-fvs
! tdv_dy = t(i,j,k)*(mvn-mvs)
! vdt_dy = dvt_dy - tdv_dy
! advy1(i,j) = advy1(i,j) - vdt_dy
hvn = 0.5*(mvn+abs(mvn))*hm0(i,j)
& + 0.5*(mvn-abs(mvn))*hm0(i,j+1)
hvs = 0.5*(mvs+abs(mvs))*hm0(i,j-1)
& + 0.5*(mvs-abs(mvs))*hm0(i,j)
hvn = hvn - mvn*hm0(i,j)
hvs = hvs - mvs*hm0(i,j)
fvn = 0.5*(mvn+abs(mvn))*t(i,j,k)
& + 0.5*(mvn-abs(mvn))*t(i,j+1,k)
fvs = 0.5*(mvs+abs(mvs))*t(i,j-1,k)
& + 0.5*(mvs-abs(mvs))*t(i,j,k)
fvn = fvn - mvn*t(i,j,k)
fvs = fvs - mvs*t(i,j,k)
gvn = 0.5*(mvn+abs(mvn))*s(i,j,k)
& + 0.5*(mvn-abs(mvn))*s(i,j+1,k)
gvs = 0.5*(mvs+abs(mvs))*s(i,j-1,k)
& + 0.5*(mvs-abs(mvs))*s(i,j,k)
gvn = gvn - mvn*s(i,j,k)
gvs = gvs - mvs*s(i,j,k)
if (k <= kln) then
hmn = hmn + dz(k)
hvmn = hvmn + hvn
fvmn = fvmn + fvn
gvmn = gvmn + gvn
else
hbn = hbn + dz(k)
hvbn = hvbn + hvn
fvbn = fvbn + fvn
gvbn = gvbn + gvn
end if
if (k <= kls) then
hms = hms + dz(k)
hvms = hvms + hvs
fvms = fvms + fvs
gvms = gvms + gvs
else
hbs = hbs + dz(k)
hvbs = hvbs + hvs
fvbs = fvbs + fvs
gvbs = gvbs + gvs
end if
! ... w(T)z = (wT)z - T(w)z
! dwt_dz = fwt-fwb
! tdw_dz = t(i,j,k)*(mwt-mwb)
! wdt_dz = dwt_dz - tdw_dz
! advz1(i,j) = advz1(i,j) - wdt_dz
fwt = 0.5*(mwt+abs(mwt))*t(i,j,k)
& + 0.5*(mwt-abs(mwt))*t(i,j,km1)
fwb = 0.5*(mwb+abs(mwb))*t(i,j,kp1)
& + 0.5*(mwb-abs(mwb))*t(i,j,k)
fwt = fwt - mwt*t(i,j,k)
fwb = fwb - mwb*t(i,j,k)
gwt = 0.5*(mwt+abs(mwt))*s(i,j,k)
& + 0.5*(mwt-abs(mwt))*s(i,j,km1)
gwb = 0.5*(mwb+abs(mwb))*s(i,j,kp1)
& + 0.5*(mwb-abs(mwb))*s(i,j,k)
gwt = gwt - mwt*s(i,j,k)
gwb = gwb - mwb*s(i,j,k)
advz1(i,j) = advz1(i,j) - (fwt-fwb)
advz2(i,j) = advz1(i,j) - (gwt-gwb)
hm = hm + dz(k)
end do
if (hme > 0.0) fume = fume/hme
if (hmw > 0.0) fumw = fumw/hmw
if (hmn > 0.0) fvmn = fvmn/hmn
if (hms > 0.0) fvms = fvms/hms
if (hbe > 0.0) fube = fube/hbe
if (hbw > 0.0) fubw = fubw/hbw
if (hbn > 0.0) fvbn = fvbn/hbn
if (hbs > 0.0) fvbs = fvbs/hbs
if (hme > 0.0) gume = gume/hme
if (hmw > 0.0) gumw = gumw/hmw
if (hmn > 0.0) gvmn = gvmn/hmn
if (hms > 0.0) gvms = gvms/hms
if (hbe > 0.0) gube = gube/hbe
if (hbw > 0.0) gubw = gubw/hbw
if (hbn > 0.0) gvbn = gvbn/hbn
if (hbs > 0.0) gvbs = gvbs/hbs
advx1(i,j) = - (fume-fumw)/(dxt(i,j)*dyt(i,j))*dt
advy1(i,j) = - (fvmn-fvms)/(dxt(i,j)*dyt(i,j))*dt
subx1(i,j) = - (fube-fubw)/(dxt(i,j)*dyt(i,j))*dt
suby1(i,j) = - (fvbn-fvbs)/(dxt(i,j)*dyt(i,j))*dt
advz1(i,j) = advz1(i,j)/(dxt(i,j)*dyt(i,j)*hm)*dt
advx2(i,j) = - (gume-gumw)/(dxt(i,j)*dyt(i,j))*dt
advy2(i,j) = - (gvmn-gvms)/(dxt(i,j)*dyt(i,j))*dt
subx2(i,j) = - (gube-gubw)/(dxt(i,j)*dyt(i,j))*dt
suby2(i,j) = - (gvbn-gvbs)/(dxt(i,j)*dyt(i,j))*dt
advz2(i,j) = advz2(i,j)/(dxt(i,j)*dyt(i,j)*hm)*dt
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
wet = undef
wes = undef
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef) cycle
wet(i,j) = 0.0
wes(i,j) = 0.0
if (hm0(i,j) > hm1(i,j)) then
tmld = 0
smld = 0
hm = 0
do k = 1, imld1(i,j)
tmld = tmld + t1(i,j,k)*dz(k)
smld = smld + s1(i,j,k)*dz(k)
hm = hm + dz(k)
end do
tmld = tmld/hm
smld = smld/hm
tsub = 0
ssub = 0
hm = 0
do k = imld1(i,j)+1, imld1(i,j)+1 !imld(i,j)
tsub = tsub + t1(i,j,k)*dz(k)
ssub = ssub + s1(i,j,k)*dz(k)
hm = hm + dz(k)
end do
tsub = tsub/hm
ssub = ssub/hm
wet(i,j) = wet(i,j) -(hm0(i,j)-hm1(i,j))/dt2/hm0(i,j)
& *(tmld-tsub)*dt
wes(i,j) = wes(i,j) -(hm0(i,j)-hm1(i,j))/dt2/hm0(i,j)
& *(smld-ssub)*dt
end if
if (hm2(i,j) > hm0(i,j)) then
tmld = 0
smld = 0
hm = 0
do k = 1, imld(i,j)
tmld = tmld + t(i,j,k)*dz(k)
smld = smld + s(i,j,k)*dz(k)
hm = hm + dz(k)
end do
tmld = tmld/hm
smld = smld/hm
tsub = 0
ssub = 0
hm = 0
do k = imld(i,j)+1, imld(i,j)+1 !imld2(i,j)
tsub = tsub + t(i,j,k)*dz(k)
ssub = ssub + s(i,j,k)*dz(k)
hm = hm + dz(k)
end do
tsub = tsub/hm
ssub = ssub/hm
wet(i,j) = wet(i,j) -(hm2(i,j)-hm0(i,j))/dt2/hm0(i,j)
& *(tmld-tsub)*dt
wes(i,j) = wes(i,j) -(hm2(i,j)-hm0(i,j))/dt2/hm0(i,j)
& *(smld-ssub)*dt
end if
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
term1 = 0
term2 = 0
jrec = mon
do n = 1, nvar
kb = nlev
if (n >= nvar) kb = 1
buf = undef
do k = 1, kb
call read_nc_data(tfile,trim(var(n)),
& nlon,nlat,buf(:,:,k),k,jrec,undef)
end do
do j = 1, nlat
do i = 1, nlon
do k = 1, imld(i,j)
if (buf(i,j,k) <= undef+1) exit
term1(i,j,n) = term1(i,j,n) + buf(i,j,k)*dz(k)*dt
end do
end do
end do
kb = nlev
if (n >= nvar-1) kb = 1
buf = undef
do k = 1, kb
call read_nc_data(sfile,trim(var(n)),
& nlon,nlat,buf(:,:,k),k,jrec,undef)
end do
do j = 1, nlat
do i = 1, nlon
do k = 1, imld(i,j)
if (buf(i,j,k) <= undef+1) exit
term2(i,j,n) = term2(i,j,n) + buf(i,j,k)*dz(k)*dt
end do
end do
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
tent = undef
sent = undef
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
tent(i,j) = 0
sent(i,j) = 0
hm = 0
do k = 1, imld(i,j)
tent(i,j) = tent(i,j) + (t2(i,j,k)-t1(i,j,k))/dt2*dz(k)
sent(i,j) = sent(i,j) + (s2(i,j,k)-s1(i,j,k))/dt2*dz(k)
hm = hm + dz(k)
end do
if (hm > 0) then
tent(i,j) = tent(i,j)/hm
sent(i,j) = sent(i,j)/hm
else
stop 'error'
tent(i,j) = undef
sent(i,j) = undef
end if
end do
end do
term1(:,:,nvar+1) = undef
term1(:,:,nvar+2) = undef
term1(:,:,nvar+3) = undef
term2(:,:,nvar+1) = undef
term2(:,:,nvar+2) = undef
term2(:,:,nvar+3) = undef
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
if (imld(i,j) == 0) cycle
term1(i,j,nvar+1) = 0
term1(i,j,nvar+2) = 0
term1(i,j,nvar+3) = 0
term2(i,j,nvar+1) = 0
term2(i,j,nvar+2) = 0
term2(i,j,nvar+3) = 0
do n = 1, nvar
term1(i,j,nvar+1) = term1(i,j,nvar+1) + term1(i,j,n)
term2(i,j,nvar+1) = term2(i,j,nvar+1) + term2(i,j,n)
end do
hm = 0
do k = 1, imld(i,j)
hm = hm + dz(k)
end do
do n = 1, nvar+1
term1(i,j,n) = term1(i,j,n)/hm
term2(i,j,n) = term2(i,j,n)/hm
end do
term1(i,j,nvar+2) = (tm2(i,j) - tm1(i,j))/dt2
term2(i,j,nvar+2) = (sm2(i,j) - sm1(i,j))/dt2
term1(i,j,nvar+3) = (tm2(i,j) - tm1(i,j))/dt2
& - tent(i,j)
term2(i,j,nvar+3) = (sm2(i,j) - sm1(i,j))/dt2
& - sent(i,j)
term1(i,j,nvar+2) = term1(i,j,nvar+2)*dt
term2(i,j,nvar+2) = term2(i,j,nvar+2)*dt
term1(i,j,nvar+3) = term1(i,j,nvar+3)*dt
term2(i,j,nvar+3) = term2(i,j,nvar+3)*dt
end do
end do
! ... define output
write(ofile,'(a,i4.4,i2.2,a)')
& 'dat/mld.',iyr,mon,'.bin'
open (90,file=trim(ofile),form='unformatted',
& access='direct',recl=nlon*nlat*4)
irec = 0
irec = irec + 1
write(90,rec=irec) tm
irec = irec + 1
write(90,rec=irec) sm
irec = irec + 1
write(90,rec=irec) mld
irec = irec + 1
write(90,rec=irec) ild
irec = irec + 1
write(90,rec=irec) blt
irec = irec + 1
write(90,rec=irec) d20
irec = irec + 1
write(90,rec=irec) um
irec = irec + 1
write(90,rec=irec) vm
do n = 1, nvar+3
irec = irec + 1
write(90,rec=irec) term1(:,:,n)
end do
do n = 1, nvar+3
irec = irec + 1
write(90,rec=irec) term2(:,:,n)
end do
irec = irec + 1
write(90,rec=irec) ((advx1(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((advy1(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((advz1(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((subx1(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((suby1(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((wet(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((advx2(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((advy2(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((advz2(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((subx2(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((suby2(i,j),i=1,nlon),j=1,nlat)
irec = irec + 1
write(90,rec=irec) ((wes(i,j),i=1,nlon),j=1,nlat)
do k = 1, nlev
irec = irec + 1
write(90,rec=irec) ((t(i,j,k),i=1,nlon),j=1,nlat)
end do
do k = 1, nlev
irec = irec + 1
write(90,rec=irec) ((s(i,j,k),i=1,nlon),j=1,nlat)
end do
do k = 1, nlev
irec = irec + 1
write(90,rec=irec) ((r(i,j,k),i=1,nlon),j=1,nlat)
end do
end do
end do
stop
end
subroutine mixed_layer(nlon,nlat,nlev,zw,dz,t,s,r,
& imld,mld,ild,ivd,blt,tm,sm,undef)
! ======================================================================
!
! computes mixed layer depth
!
! ======================================================================
implicit none
integer :: nlon, nlat, nlev
real :: zw(0:nlev)
real :: dz(nlev)
real :: t(nlon,nlat,nlev)
real :: s(nlon,nlat,nlev)
real :: r(nlon,nlat,nlev)
real :: tm(nlon,nlat)
real :: sm(nlon,nlat)
integer :: imld(nlon,nlat)
real :: mld(nlon,nlat)
real :: ild(nlon,nlat)
real :: ivd(nlon,nlat)
real :: blt(nlon,nlat)
integer :: ild2(nlon,nlat)
real :: undef
real(4) :: hm
real(8) :: del_r0 = 0.125
real(8) :: del_t0 = 0.2
real(8) :: del_r
real(8) :: del_t
real :: t1, t2, r1, r2, s1, s2
integer :: i, j, k, kb, kv
! ... set
imld = 0
mld = undef
ild = undef
ivd = undef
blt = undef
ild2 = 0
! ... density
do j = 1, nlat
do i = 1, nlon
do k = 1, nlev
if (t(i,j,k) <= undef+1) then
r(i,j,k) = undef
else
call dens(t(i,j,k),s(i,j,k),r(i,j,k))
end if
end do
end do
end do
! ... mixed layer depth
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
t1 = t(i,j,1)
t2 = t(i,j,1)-del_t0
s1 = s(i,j,1)
s2 = s(i,j,1)
call dens(t1,s1,r1)
call dens(t2,s2,r2)
del_r0 = r2-r1
kb = 0
do k = 1, nlev
if (r(i,j,k) <= undef+1) exit
del_r = r(i,j,k) - r(i,j,1)
if (del_r > del_r0) exit
kb = k
end do
imld(i,j) = kb
mld(i,j) = zw(kb)
end do
end do
! ... isothermal layer depth
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
kb = 0
do k = 1, nlev
if (t(i,j,k) <= undef+1) exit
del_t = t(i,j,1) - t(i,j,k)
if (del_t > del_t0) exit
kb = k
end do
ild(i,j) = zw(kb)
ild2(i,j) = kb
end do
end do
! ... check temperature inversion
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
kv = 0
do k = kb, 1, -1
if (t(i,j,k) <= undef+1) exit
del_t = t(i,j,k) - t(i,j,1)
if (del_t <= 0.0) cycle
kv = k
exit
end do
ivd(i,j) = zw(kv)
if (kv == 0) ivd(i,j) = undef
end do
end do
! ... reset mld
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
if (mld(i,j) > ild(i,j)) then
mld(i,j) = ild(i,j)
imld(i,j) = ild2(i,j)
end if
end do
end do
! ... test : fixed mld
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
kb = 0
do k = 1, 7
if (t(i,j,k) < undef+1) exit
kb = k
end do
! imld(i,j) = kb
! mld(i,j) = zw(kb)
end do
end do
! ... barrier layer thickness
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
blt(i,j) = ild(i,j) - mld(i,j)
end do
end do
! ... tempeature and salinity in the mixed layer
tm = undef
sm = undef
do j = 1, nlat
do i = 1, nlon
if (t(i,j,1) <= undef+1) cycle
tm(i,j) = 0
sm(i,j) = 0
hm = 0
do k = 1, imld(i,j)
if (t(i,j,k) <= undef+1) cycle
if (s(i,j,k) <= undef+1) cycle
tm(i,j) = tm(i,j) + t(i,j,k)*dz(k)
sm(i,j) = sm(i,j) + s(i,j,k)*dz(k)
hm = hm + dz(k)
end do
if (hm > 0) then
tm(i,j) = tm(i,j)/hm
sm(i,j) = sm(i,j)/hm
else
stop 'error in tempeature and salinity in the mixed layer'
tm(i,j) = undef
sm(i,j) = undef
end if
end do
end do
return
end
subroutine dens(t,s,rho)
! ======================================================================
!
! computes density based on unescos formula as given in Gill
!
! ======================================================================
implicit none
real :: t
real :: s
real :: rho
! ... local
real(8) :: rhor
real(8) :: tr1, tr2, tr3, tr4, tr5
real(8) :: sr1, sr2
tr1 = t
tr2 = tr1*tr1
tr3 = tr1*tr2
tr4 = tr1*tr3
tr5 = tr1*tr4
sr1 = s
sr2 = sr1*sr1
rhor = 999.842594
rhor = rhor + ( 6.793952e-2*tr1
& - 9.095290e-3*tr2
& + 1.001685e-4*tr3
& - 1.120083e-6*tr4
& + 6.536332e-9*tr5)
rhor = rhor + ( 0.824493
& - 4.0899e-3 *tr1
& + 7.6438e-5 *tr2
& - 8.2467e-7 *tr3
& + 5.3875e-9 *tr4)*sr1
rhor = rhor + (-5.72466e-3
& + 1.0227e-4 *tr1
& - 1.6546e-6 *tr2)*abs(sr1)**1.5
rhor = rhor + ( 4.8314e-4 )*sr2
rho = rhor - 1000.0
! rho = rho - 1000.0
return
end
! ======================================================================
!
! ======================================================================
subroutine convect(nz,t,s,dz,mask)
implicit none
!---------------------------------------------------------------------
! convective adjustment for gravitationally unstable columns
! of water
!---------------------------------------------------------------------
integer :: nz
real :: t(nz), s(nz)
real :: zw(nz), dz(nz)
real :: mask(nz)
real :: rho(nz)
real :: zm
integer :: n, k, ncon
ncon = nz-1
do n = 1, ncon
do k = 1, 1
call density(t(k),s(k),rho(k))
end do
do k = 1, nz-1
if (mask(k+1) == 0) exit
call density(t(k+1),s(k+1),rho(k+1))
if (rho(k) > rho(k+1)) then
zm = dz(k) + dz(k+1)
t(k) = (t(k)*dz(k) + t(k+1)*dz(k+1))/zm
t(k+1) = t(k)
s(k) = (s(k)*dz(k) + s(k+1)*dz(k+1))/zm
s(k+1) = s(k)
call density(t(k+1),s(k+1),rho(k+1))
end if
end do
end do
return
end
! ======================================================================
subroutine density(t,s,rho)
! ======================================================================
!
! computes density based on unescos formula as given in Gill
!
! ======================================================================
implicit none
real(4) :: t, s, rho
real(8) :: rhor
real(8) :: tr1, tr2, tr3, tr4, tr5
real(8) :: sr1, sr2
tr1 = t
tr2 = tr1*tr1
tr3 = tr1*tr2
tr4 = tr1*tr3
tr5 = tr1*tr4
sr1 = s
sr2 = sr1*sr1
! rhor = 999.842594
rhor = 0.0
rhor = rhor + ( 6.793952e-2*tr1
& - 9.095290e-3*tr2
& + 1.001685e-4*tr3
& - 1.120083e-6*tr4
& + 6.536332e-9*tr5)
rhor = rhor + ( 0.824493
& - 4.0899e-3 *tr1
& + 7.6438e-5 *tr2
& - 8.2467e-7 *tr3
& + 5.3875e-9 *tr4)*sr1
rhor = rhor + (-5.72466e-3
& + 1.0227e-4 *tr1
& - 1.6546e-6 *tr2)*sr1**1.5
rhor = rhor + ( 4.8314e-4 )*sr2
rho = rhor
return
end
EOF
set FC = "ifort"
#set NCLIB = "$home/local/netcdf/*.o -L/usr/local/netcdf-3.6.2/lib -lnetcdf"
set NCLIB = "../netcdf/*.o -L/usr/local/lib -lnetcdf"
set OPT = "-assume byterecl -ftrapuv -fpe0 -traceback -CB"
set DBG =
#set DBG = "-g -fbounds-check -o Wuninitialized -fbacktrace -ffpe-trap=invalid,zero,overflow"
#set DBG = "-g -fbounds-check -Wall -fbacktrace -ffpe-trap=invalid,zero,overflow"
$FC $DBG $OPT a.f $NCLIB && ./a.out
#gfortran $DBG a.f $NCLIB && ./a.out
#/bin/rm a.f a.out