diff --git a/.gitignore b/.gitignore index 8293d27e..17fb9116 100644 --- a/.gitignore +++ b/.gitignore @@ -604,3 +604,7 @@ hycom/hycom_ALL/hycom_2.2.72_ALL/topo/src_2.2.72/topo_smooth_skip hycom/hycom_ALL/hycom_2.2.72_ALL/topo/src_2.2.72/topo_subset hycom/hycom_ALL/hycom_2.2.72_ALL/topo/src_2.2.72/topo_zcells hycom/hycom_ALL/hycom_2.2.72_ALL/topo/src_2.2.72/topo_zthin + +hycom/MSCPROGS/src/Hyc2proj/TMP/* +hycom/MSCPROGS/src/IceDriftNC/TMP/* + diff --git a/hycom/MSCPROGS/src/Average/README.md b/hycom/MSCPROGS/src/Average/README.md new file mode 100644 index 00000000..10084c8a --- /dev/null +++ b/hycom/MSCPROGS/src/Average/README.md @@ -0,0 +1,52 @@ +## hycom_mean + +```hycom_mean``` is F90 version of ```hycom/hycom_ALL/hycom_2.2.72_ALL/meanstd/src/hycom_mean``` for making an ensemble average of hycom archive files. Note that other options, such as meanstd or meansq, available in the original ```hycom_mean``` are not available at this moment. + +### Step1. Prepare configuration file ```mean_hycom.in```: + +``` + 50 'kk ' = number of layers involved + 0 'meansq' = form meansq, rather than mean (0=F,1=T) + 2 'narchs' = number of archives to read (==0 to end input) +file_mem001.a +file_mem002.a + 0 'narchs' = number of archives to read (==0 to end input) +file_mean +``` + +where ```narchs``` is number of hycom ab files to be averaged, + +``` +file_mem001.a +file_mem002.a +``` + +are a list of input hycom a files, + +``` +file_mean +``` + +specifies a name of mean ab files: ```file_mean.a```,```file_mean.b```, generated by ```hycom_mean```. + +### Step 2. Run script at command line: + +```bash +hycom_mean < mean_hycom.in +``` + +### Notes + +Upon change of BGC variables list on hycom b file, you need to modify variable registraion in ```mod_mean.F90```. + +In order to use ```hycom_mean``` with HYCOM-CICE setup without ECOSMO, you need to turn off ECOSMO option set in ```hycom_mean.F90``` by changing: + +``` +lecosmo = .true. +``` + +to + +``` +lecosmo = .false. +``` diff --git a/hycom/MSCPROGS/src/Average/bigrid.F90 b/hycom/MSCPROGS/src/Average/bigrid.F90 new file mode 100644 index 00000000..b3520b19 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/bigrid.F90 @@ -0,0 +1,72 @@ +subroutine bigrid(depth) + use mod_mean ! HYCOM mean array interface + use mod_za ! HYCOM array I/O interface + implicit none + + integer :: i, j + real, dimension(0:ii,0:jj) :: depth + + ! --- set land/sea masks for irregular basin in c-grid configuration + ! --- q,u,v,p are vorticity, u-velocity, v-velocity, and mass points, resp. + ! --- 'depth' = basin depth array, zero values indicate land + + integer, allocatable, dimension (:,:) :: ip0 + + allocate( ip0(0:ii,0:jj) ) + + write (lp,'(a,6i6)') 'bigrid called with ii,jj,ii1,jj1 =', ii,jj,ii1,jj1 + + ! --- mass points are defined where water depth is greater than zero + do i=0,ii + do j=0,jj + if (depth(i,j).gt.0.) then + ip0(i,j)=1 + else + ip0(i,j)=0 + endif + enddo + enddo + + do i=1,ii + do j=1,jj + ip(i,j)=ip0(i,j) + enddo + enddo + + ! --- u,v points are located halfway between any 2 adjoining mass points + do j=1,jj + do i=1,ii + if (ip0(i-1,j).gt.0.and.ip0(i,j).gt.0) then + iu(i,j)=1 + else + iu(i,j)=0 + endif + if (ip0(i,j-1).gt.0.and.ip0(i,j).gt.0) then + iv(i,j)=1 + else + iv(i,j)=0 + endif + enddo + enddo + + do j=1,jj + do i=1,ii + iq(i,j)=0 + + ! --- 'interior' q points require water on all 4 sides. + if (min0(ip(i,j),ip(i-1,j),ip(i,j-1),ip(i-1,j-1)).gt.0) then + iq(i,j)=1 + endif + + ! --- 'promontory' q points require water on 3 + ! --- (or at least 2 diametrically opposed) sides + if ((ip(i,j).gt.0.and.ip(i-1,j-1).gt.0).or.(ip(i-1,j).gt.0.and.ip(i,j-1).gt.0)) then + iq(i,j)=1 + endif + enddo + enddo + + deallocate( ip0 ) + + return +end subroutine bigrid diff --git a/hycom/MSCPROGS/src/Average/blkin.F90 b/hycom/MSCPROGS/src/Average/blkin.F90 new file mode 100644 index 00000000..51c9f6a7 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/blkin.F90 @@ -0,0 +1,366 @@ +module mod_blkin + use mod_za ! HYCOM array I/O interface + implicit none + + contains + + subroutine blkind(dvar,cvar,cfmt) + implicit none + + real(kind=8) :: dvar + character(len=*) :: cfmt + character(len=6) :: cvar + + ! integer ::lp + ! common/linepr/lp + + ! read in one d.p. value from stdin + + character(len=6) :: cvarin + + read(*,*) dvar,cvarin + write(lp,cfmt) cvarin,dvar + call flush(lp) + + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkind - input ',cvarin,' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + end subroutine blkind + + subroutine blkinr(rvar,cvar,cfmt) + implicit none + + real(kind=4) :: rvar + character(len=*) :: cfmt + character(len=6) :: cvar + + !integer :: lp + !common/linepr/lp + + !read in one real value from stdin + + character(len=6) :: cvarin + + read(*,*) rvar,cvarin + write(lp,cfmt) cvarin,rvar + call flush(lp) + + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkinr - input ',cvarin,' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + end subroutine blkinr + + subroutine blkinr2(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2) + implicit none + + real(kind=4) :: rvar + integer :: nvar + character(len=*) :: cfmt1, cfmt2 + character(len=6) :: cvar1,cvar2 + + ! integer :: lp + ! common/linepr/lp + + !read in one real value from stdin, + !identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) + + character*6 cvarin + + read(*,*) rvar,cvarin + + if (cvar1.eq.cvarin) then + nvar = 1 + write(lp,cfmt1) cvarin,rvar + call flush(lp) + elseif (cvar2.eq.cvarin) then + nvar = 2 + write(lp,cfmt2) cvarin,rvar + call flush(lp) + else + write(lp,*) + write(lp,*) 'error in blkinr2 - input ',cvarin,' but should be ',cvar1,' or ',cvar2 + write(lp,*) + call flush(lp) + stop + endif + return + end subroutine blkinr2 + + subroutine blkinr3(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3) + implicit none + + real(kind=4) :: rvar + integer :: nvar + character(len=*) :: cfmt1,cfmt2,cfmt3 + character(len=6) :: cvar1,cvar2,cvar3 + + ! integer :: lp + ! common/linepr/lp + + ! read in one of three possible real values from stdin, + ! identified as either cvar1 (return nvar=1) or + ! cvar2 (return nvar=2) or + ! cvar3 (return nvar=3) + + call blkinr9(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3, & + 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & + 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & + 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & + 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & + 'XXXXXX','("blkinr: ",a6," =",f11.4," ")', & + 'XXXXXX','("blkinr: ",a6," =",f11.4," ")') + return + end subroutine blkinr3 + + subroutine blkinr9(rvar,nvar, & + cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3, & + cvar4,cfmt4,cvar5,cfmt5,cvar6,cfmt6, & + cvar7,cfmt7,cvar8,cfmt8,cvar9,cfmt9) + implicit none + + real(kind=4) :: rvar + integer :: nvar + character(len=6) :: cvar1,cvar2,cvar3,cvar4,cvar5,cvar6,cvar7,cvar8,cvar9 + character(len=*) :: cfmt1,cfmt2,cfmt3,cfmt4,cfmt5,cfmt6,cfmt7,cfmt8,cfmt9 + + ! integer :: lp + ! common/linepr/lp + + ! read in one of nine possible real values from stdin, + ! identified as either cvar1 (return nvar=1) or + ! cvar2 (return nvar=2) or + ! ... + ! cvar9 (return nvar=9) + + character*6 cvarin + + read(*,*) rvar,cvarin + + if (cvar1.eq.cvarin) then + nvar = 1 + write(lp,cfmt1) cvarin,rvar + call flush(lp) + elseif (cvar2.eq.cvarin) then + nvar = 2 + write(lp,cfmt2) cvarin,rvar + call flush(lp) + elseif (cvar3.eq.cvarin) then + nvar = 3 + write(lp,cfmt3) cvarin,rvar + call flush(lp) + elseif (cvar4.eq.cvarin) then + nvar = 4 + write(lp,cfmt4) cvarin,rvar + call flush(lp) + elseif (cvar5.eq.cvarin) then + nvar = 5 + write(lp,cfmt5) cvarin,rvar + call flush(lp) + elseif (cvar6.eq.cvarin) then + nvar = 6 + write(lp,cfmt6) cvarin,rvar + call flush(lp) + elseif (cvar7.eq.cvarin) then + nvar = 7 + write(lp,cfmt7) cvarin,rvar + call flush(lp) + elseif (cvar8.eq.cvarin) then + nvar = 8 + write(lp,cfmt8) cvarin,rvar + call flush(lp) + elseif (cvar9.eq.cvarin) then + nvar = 9 + write(lp,cfmt9) cvarin,rvar + call flush(lp) + else + write(lp,*) + write(lp,*) 'error in blkinr9 - input ',cvarin,' but should be:' + write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or' + write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or' + write(lp,*) cvar7,' or ',cvar8,' or ',cvar9 + write(lp,*) + call flush(lp) + stop + endif + return + end subroutine blkinr9 + + subroutine blkini(ivar,cvar) + implicit none + + integer :: ivar + character(len=6) :: cvar + + ! integer :: lp + ! common/linepr/lp + + ! read in one integer value from stdin + + character(len=6) :: cvarin + + read(*,*) ivar,cvarin + write(lp,6000) cvarin,ivar + call flush(lp) + + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkini - input ',cvarin,' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkini: ',a6,' =',i6) + end subroutine blkini + + subroutine blkini2(ivar,nvar,cvar1,cvar2) + implicit none + + integer :: ivar,nvar + character(len=6) :: cvar1,cvar2 + + ! integer :: lp + ! common/linepr/lp + + ! read in one integer value from stdin + ! identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) + + character(len=6) :: cvarin + + read(*,*) ivar,cvarin + write(lp,6000) cvarin,ivar + call flush(lp) + + if (cvarin.eq.cvar1) then + nvar = 1 + elseif (cvarin.eq.cvar2) then + nvar = 2 + else + write(lp,*) + write(lp,*) 'error in blkini2 - input ',cvarin,' but should be ',cvar1,' or ',cvar2 + write(lp,*) + call flush(lp) + stop + endif + return +6000 format('blkini: ',a6,' =',i6) + end subroutine blkini2 + + subroutine blkini3(ivar,nvar,cvar1,cvar2,cvar3) + implicit none + + integer ::ivar,nvar + character(len=6) :: cvar1,cvar2,cvar3 + + ! integer :: lp + ! common/linepr/lp + + ! read in one of three possible integer values from stdin + ! identified as either cvar1 (return nvar=1) or + ! cvar2 (return nvar=2) or + ! cvar3 (return nvar=3) + + call blkini9(ivar,nvar,cvar1,cvar2,cvar3, & + 'XXXXXX','XXXXXX','XXXXXX', & + 'XXXXXX','XXXXXX','XXXXXX') + return + end subroutine blkini3 + + subroutine blkini9(ivar,nvar,cvar1,cvar2,cvar3, & + cvar4,cvar5,cvar6, & + cvar7,cvar8,cvar9) + implicit none + + integer :: ivar,nvar + character(len=6) :: cvar1,cvar2,cvar3,cvar4,cvar5,cvar6,cvar7,cvar8,cvar9 + + ! integer :: lp + ! common/linepr/lp + + ! read in one of nine possible integer values from stdin + ! identified as either cvar1 (return nvar=1) or + ! cvar2 (return nvar=2) or + ! cvar3 (return nvar=3) or + ! ... + ! cvar9 (return nvar=9) + + character(len=6) :: cvarin + + read(*,*) ivar,cvarin + write(lp,6000) cvarin,ivar + call flush(lp) + + if (cvarin.eq.cvar1) then + nvar = 1 + elseif (cvarin.eq.cvar2) then + nvar = 2 + elseif (cvarin.eq.cvar3) then + nvar = 3 + elseif (cvar4.eq.cvarin) then + nvar = 4 + elseif (cvar5.eq.cvarin) then + nvar = 5 + elseif (cvar6.eq.cvarin) then + nvar = 6 + elseif (cvar7.eq.cvarin) then + nvar = 7 + elseif (cvar8.eq.cvarin) then + nvar = 8 + elseif (cvar9.eq.cvarin) then + nvar = 9 + else + write(lp,*) + write(lp,*) 'error in blkini9 - input ',cvarin,' but should be:' + write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or' + write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or' + write(lp,*) cvar7,' or ',cvar8,' or ',cvar9 + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkini: ',a6,' =',i6) + end subroutine blkini9 + + subroutine blkinl(lvar,cvar) + implicit none + + logical :: lvar + character(len=6) :: cvar + + ! integer :: lp + ! common/linepr/lp + + ! read in one logical value from stdin + ! due to a SGI bug for logical I/O: read in an integer 0=F,1=T + + character*6 cvarin + integer ivar + + read(*,*) ivar,cvarin + lvar = ivar .ne. 0 + write(lp,6000) cvarin,lvar + call flush(lp) + + if (cvar.ne.cvarin) then + write(lp,*) + write(lp,*) 'error in blkinr - input ',cvarin,' but should be ',cvar + write(lp,*) + call flush(lp) + stop + endif + return + 6000 format('blkinl: ',a6,' =',l6) + end subroutine blkinl + +end module mod_blkin diff --git a/hycom/MSCPROGS/src/Average/extrct.F90 b/hycom/MSCPROGS/src/Average/extrct.F90 new file mode 100644 index 00000000..8af1cfa5 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/extrct.F90 @@ -0,0 +1,20 @@ +subroutine extrct(work,n,m,io,jo,array,no,mo) + implicit none + + integer :: n,m,io,jo,no,mo + real :: work(n,m),array(no,mo) +! +! --- array = work(io:io+no-1,jo:jo+mo-1) +! +! --- this version c y c l i c in i + + integer :: i,j + + do j=1,min(mo,m-jo+1) + do i=1,no + array(i,j)=work(mod(io+i-2,n)+1,jo+j-1) + enddo + enddo + + return +end subroutine extrct diff --git a/hycom/MSCPROGS/src/Average/getdat.F90 b/hycom/MSCPROGS/src/Average/getdat.F90 new file mode 100644 index 00000000..c6a62fc9 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/getdat.F90 @@ -0,0 +1,931 @@ +subroutine getdat(flnm,time,iweight,mntype,iexpt,yrflag,thbase) + use mod_mean ! HYCOM mean array interface + use mod_za ! HYCOM array I/O interface + + character(len=*) :: flnm + real*8, dimension(3) :: time + real :: thbase + !logical, intent(in) :: icegln,trcout + integer :: iweight,iexpt,mntype,yrflag + + ! --- read model fields and extract portion of global fields. + ! --- HYCOM 2.0 array I/O archive or mean/mnsq archive file. + + character(len=80) :: cline + character(len=79), dimension(5) :: preambl + character(len=6) :: cvarin + real :: hminb,hmaxb + real :: thet + integer :: i,j,k,iversn,l,layer,sigver + integer :: ios,ktr,ntr + logical :: nodens + + real, allocatable :: work(:,:) + + data ni/14/ + + allocate( work(idm,jdm) ) + + l = len_trim(flnm) + open (unit=ni,file=flnm(1:l-2)//'.b',form='formatted',status='old',action='read') + call zaiopf(flnm(1:l-2)//'.a','old', ni) + + read( ni,'(a80/a80/a80/a80)') ctitle + write(lp,'(a80/a80/a80/a80)') ctitle + + read( ni,*) iversn,cvarin + write(lp,*) cvarin,' = ',iversn + if (cvarin.ne.'iversn') then + write(lp,*) + write(lp,*) 'error in getdat - input ',cvarin,' but should be iversn' + write(lp,*) + call flush(lp) + stop + endif + + read( ni,*) iexpt,cvarin + write(lp,*) cvarin,' = ',iexpt + if (cvarin.ne.'iexpt ') then + write(lp,*) + write(lp,*) 'error in getdat - input ',cvarin,' but should be iexpt ' + write(lp,*) + call flush(lp) + stop + endif + + read( ni,*) yrflag,cvarin + write(lp,*) cvarin,' = ',yrflag + if (cvarin.ne.'yrflag') then + write(lp,*) + write(lp,*) 'error in getdat - input ',cvarin,' but should be yrflag' + write(lp,*) + call flush(lp) + stop + endif + + read( ni,*) idmtst,cvarin + write(lp,*) cvarin,' = ',idmtst + if (cvarin.ne.'idm ') then + write(lp,*) + write(lp,*) 'error in getdat - input ',cvarin,' but should be idm ' + write(lp,*) + call flush(lp) + stop + endif + + read( ni,*) jdmtst,cvarin + write(lp,*) cvarin,' = ',jdmtst + if (cvarin.ne.'jdm ') then + write(lp,*) + write(lp,*) 'error in getdat - input ',cvarin,' but should be jdm ' + write(lp,*) + call flush(lp) + stop + endif + + if (idmtst.ne.idm .or. jdmtst.ne.jdm) then + write(lp,*) + write(lp,*) 'error in getdat - input idm,jdm',' not consistent with parameters' + write(lp,*) 'idm,jdm = ',idm, jdm, ' (REGION.h)' + write(lp,*) 'idm,jdm = ',idmtst,jdmtst,' (input)' + write(lp,*) + call flush(lp) + stop + endif + + read( ni,'(a)') cline + write(lp,'(a)') cline(1:len_trim(cline)) + + if (cline(24:28).eq.'model') then + iweight = 1 ! standard archive + mntype = 0 + elseif (cline(24:28).eq.' mean') then + iweight = 0 ! mean archive, see below + mntype = 1 + elseif (cline(24:28).eq.' mnsq') then + iweight = 0 ! mean archive, see below + mntype = 2 + else + write(lp,*) + write(lp,*) 'error in getdat - wrong archive type.' + write(lp,*) 'only "model", " mean", or " mnsq" allowed.' + write(lp,*) + call flush(lp) + stop + endif + + !-------------------------------- + ! 'montg1 ' > montg(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(1),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,montg,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'montg ' + + ! --- detect version 2.2 normal archive files + + nodens = layer.ne.0 + if (nodens) then + sigver = layer + thbase = thet + endif + !-------------------------------- + ! 'srfhgt ' > srfht(:,:) + !-------------------------------- + read(ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(2),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,srfht,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'srfht ' + + if (lsteric) then + !-------------------------------- + ! 'steric ' > steric(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work,ni,hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,steric,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'steric ' + endif + !-------------------------------- + ! 'oneta ' > oneta(:,:) + !-------------------------------- + if (loneta) then + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,oneta,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'oneta ' + !-------------------------------- + ! 'mnoneta ' > onetaw(:,:) + !-------------------------------- + if (mntype.eq.2) then ! mnsq + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work,ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,onetaw,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'mnoneta ' + if (cline(1:8).ne.'mnoneta ') then + write(lp,*) + write(lp,*) 'error in getdat - must have mnoneta in mnsq archive' + write(lp,*) + call flush(lp) + call xcstop('error - must have mnoneta in mnsq archive') + stop + endif + else + onetaw(:,:) = oneta(:,:) + endif !mnsq:else + else + oneta(:,:) = 1.0 + onetaw(:,:) = 1.0 + endif + !-------------------------------- + ! 'surflx ' > surflx(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,surflx,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'surflx ' + !-------------------------------- + ! 'wtrflx ' > wtrflx(:,:) + !-------------------------------- + if (lwtrflx) then + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,wtrflx,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'wtrflx ' + else + wtrflx(:,:) = 0.0 + endif + !-------------------------------- + ! 'salflx ' > salflx(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,salflx,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'salflx ' + !-------------------------------- + ! 'bl_dpth ' > dpbl(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,dpbl,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'dpbl ' + !-------------------------------- + ! 'mix_dpth' > dpmixl(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,dpmixl,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'dpmixl ' + + if (.not. nodens) then + !-------------------------------- + ! 'tmix ' > tmix(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,tmix,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'tmix ' + !-------------------------------- + ! 'smix ' > smix(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,smix,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'smix ' + !-------------------------------- + ! 'thmix ' > thmix(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + thbase = thet + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,thmix,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'thmix ' + !-------------------------------- + ! 'umix ' > umix(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,umix,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'umix ' + !-------------------------------- + ! 'vmix ' > vmix(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,vmix,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'vmix ' + endif !.not. nodens + + !-------------------------------- + ! 'kemix ' > kemix(:,:) + !-------------------------------- + if (iweight.eq.0) then ! mean archive + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,kemix,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'kemix ' + endif + + ! --- is there ice? + if (icegln) then + !-------------------------------- + ! 'covice ' > covice(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,covice,ii,jj) + !-------------------------------- + ! 'thcice ' > thkice(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,thkice,ii,jj) + !-------------------------------- + ! 'temice ' > temice(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,temice,ii,jj) + else + covice(:,:) = 0.0 + thkice(:,:) = 0.0 + temice(:,:) = 0.0 + endif + !-------------------------------- + ! 'u_btrop ' > ubaro(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,ubaro,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'ubaro ' + !-------------------------------- + ! 'v_btrop ' > vbaro(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,vbaro,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'vbaro ' + !-------------------------------- + ! 'kebtrop ' > kebaro(:,:) + !-------------------------------- + if (iweight.eq.0) then ! mean archive + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,kebaro,ii,jj) + write(lp,'("input ",a," into ",a)') cline(1:8),'kebaro ' + endif + + do k=1,kk + !-------------------------------- + ! 'u-vel ' > u(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,u(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'u ',k + if (cline(1:8).ne.'u-vel. ') then + write(lp,*) + write(lp,*) 'error in getdat - layer ',k,' does not exist (kk= ',kk,')' + write(lp,*) + call flush(lp) + stop + endif + !-------------------------------- + ! 'v-vel ' > v(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .true.) + call extrct(work,idm,jdm,iorign,jorign,v(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'v ',k + !-------------------------------- + ! 'k.e. ' > ke(:,:,k) + !-------------------------------- + if (iweight.eq.0) then ! mean archive + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,ke(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'k.e. ',k + endif + !-------------------------------- + ! 'thknss ' > dp(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,dp(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'dp ',k + !-------------------------------- + ! 'mnthknss' > dw(:,:,k) + !-------------------------------- + if (mntype.eq.2) then ! mnsq + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,dw(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'dw ',k + if (cline(1:8).ne.'mnthknss') then + write(lp,*) + write(lp,*) 'error in getdat - must have mnthknss in mnsq archive' + write(lp,*) + call flush(lp) + stop + endif + else + dw(:,:,k) = dp(:,:,k) + endif + !-------------------------------- + ! 'temp ' > temp(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,temp(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'temp ',k + !-------------------------------- + ! 'salin ' > saln(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,saln(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'saln ',k + !-------------------------------- + ! 'density ' > th3d(:,:,k) + !-------------------------------- + if (.not. nodens) then + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,th3d(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'th3d ',k + else + call th3d_p(temp(1,1,k),saln(1,1,k),th3d(1,1,k),ii,jj,sigver,thbase) + write(lp,'(" ",a8,"calculate ",a,i3)') " ",'th3d ',k + if (k.eq.1) then + tmix(:,:) = temp(:,:,1) + smix(:,:) = saln(:,:,1) + thmix(:,:) = th3d(:,:,1) + umix(:,:) = u(:,:,1) + vmix(:,:) = v(:,:,1) + write(lp,'("copy ",a," into ",a)') 'temp.1 ','tmix ' + write(lp,'("copy ",a," into ",a)') 'saln.1 ','smix ' + write(lp,'("copy ",a," into ",a)') 'th3d.1 ','thmix ' + write(lp,'("copy ",a," into ",a)') ' u.1 ','umix ' + write(lp,'("copy ",a," into ",a)') ' v.1 ','vmix ' + endif !k==1 + endif !.not.nodens:else + !-------------------------------- + ! 'tracer ' > tracer(:,:,k,ktr) + !-------------------------------- + if (trcout) then + do ktr= 1,ntracr + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,tracer(1,1,k,ktr),ii,jj) + write(lp,'("input ",a," into ",a,2i3)') cline(1:8),'tracer ',k,ktr + enddo !ktr + endif !trcout + !-------------------------------- + ! 'viscty ' > visc(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,visc(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'visc ',k + !-------------------------------- + ! 't-diff ' > tdff(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,tdff(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'tdff ',k + !-------------------------------- + ! 's-diff ' > sdff(:,:,k) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,sdff(1,1,k),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'sdff ',k + + write(lp,'(a,f9.5)') 'finished reading data for layer',thet + call flush(lp) + theta(k)=thet + + enddo !k-loop + + if (bgcout) then + !-------------------------------- + ! > tracer_bgc_3d + !-------------------------------- + do k=1,kk + do ktr= 1,ntracr_bgc_3d + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,tracer_bgc_3d(1,1,k,ktr),ii,jj) + write(lp,'("input ",a," into ",a,2i3)') cline(1:8),'tracer_bgc_3d',k,ktr + enddo !ktr + enddo !k-loop + !-------------------------------- + ! > tracer_bgc_2d + !-------------------------------- + do ktr= 1,ntracr_bgc_2d + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,tracer_bgc_2d(1,1,ktr),ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'tracer_bgc_2d' + enddo !ktr + endif ! bgcout + + if (mntype.eq.1.and.icegln) then + !-------------------------------- + ! 'si_u ' > si_u(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,si_u,ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'si_u' + !-------------------------------- + ! 'si_v ' > si_v(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,si_v,ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'si_v' + !-------------------------------- + ! 'surtx ' > surtx(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,surtx,ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'surtx' + !-------------------------------- + ! 'surty ' > surty(:,:) + !-------------------------------- + read (ni,'(a)',end=6) cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) nstep,time(3),layer,thet,hminb,hmaxb + call getfld(work, ni, hminb,hmaxb, .false.) + call extrct(work,idm,jdm,iorign,jorign,surty,ii,jj) + write(lp,'("input ",a," into ",a,i3)') cline(1:8),'surty' + endif + + if (iweight.eq.0) then + iweight = nstep ! mean archive + endif + + close( unit=ni) + call zaiocl(ni) + + deallocate( work ) + + return + + ! --- unexpected end of file +6 continue + write (lp,*) '***** unexpected end of archive file *****' + call flush(lp) + stop '(e-o-f)' + +end subroutine getdat + +subroutine getdepth(dpthfil) + use mod_mean ! HYCOM mean array interface + use mod_za + + character :: dpthfil*(*) + + ! --- acquire basin depths + + character :: cline*80 + character(len=79), dimension(5) :: preambl + character(len=6) :: cvarin + real :: hminb,hmaxb + integer :: i,j,iversn,l + + real, allocatable :: work(:,:) + + allocate( work(idm,jdm) ) + + open(unit=9,file=dpthfil(1:len_trim(dpthfil))//'.b',form='formatted',status='old',action='read') + read(9, '(a79)') preambl + write(lp,'(a79)') preambl + read (9, '(a)') cline + write(lp,'(a)') cline(1:len_trim(cline)) + i = index(cline,'=') + read (cline(i+1:),*) hminb,hmaxb + close(unit=9) + + call zaiopf(dpthfil(1:len_trim(dpthfil))//'.a','old', 9) + call getfld(work, 9, hminb,hmaxb, .true.) + call zaiocl(9) + write(lp,'("read ",a," into ",a)') preambl(1)(1:8),'depths ' + do j= 1,jj + do i= 1,ii + depths(i,j) = work(i,j) + enddo + enddo + do i= 0,ii + depths(i,0) = 0.0 + enddo + do j= 0,jj + depths(0,j) = depths(ii,j) ! assumed periodic + enddo + + deallocate( work ) + + return + +end subroutine getdepth + +subroutine getfld(work, iunit, hminb,hmaxb, lzero) + use mod_za ! HYCOM array I/O interface + + ! --- read a single array + + logical :: lzero + integer :: iunit + real, dimension(idm,jdm) :: work + integer, dimension(idm,jdm) :: mask + real :: hminb,hmaxb + real :: hmina,hmaxa + + call zaiord(work,mask,.false., hmina,hmaxa, iunit) + + if (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or. abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4 ) then + write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)') & + ,'error - .a and .b files not consistent:' & + ,'.a,.b min = ',hmina,hminb,hmina-hminb & + ,'.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb + call flush(lp) + stop + endif + + if (lzero) then + do j= 1,jdm + do i= 1,idm + if (work(i,j).gt.2.0**99) then + work(i,j) = 0.0 + endif + enddo + enddo + endif + + return +end subroutine getfld + +subroutine th3d_p(temp,saln,th3d,no,mo,sigver,thbase) + implicit none + + integer no,mo,sigver + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the appropriate equation of state. + + if (sigver.eq.1) then + call th3d_p1(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.2) then + call th3d_p2(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.3) then + call th3d_p3(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.4) then + call th3d_p4(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.5) then + call th3d_p5(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.6) then + call th3d_p6(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.7) then + call th3d_p7(temp,saln,th3d,no,mo,thbase) + elseif (sigver.eq.8) then + call th3d_p8(temp,saln,th3d,no,mo,thbase) + else !unknown + th3d(:,:) = 0.0 + endif + return + +end subroutine th3d_p + +subroutine th3d_p1(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA0_7term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + return + +end subroutine th3d_p1 + +subroutine th3d_p2(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA2_7term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + return +end subroutine th3d_p2 + +subroutine th3d_p3(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA0_9term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + return +end subroutine th3d_p3 + +subroutine th3d_p4(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA2_9term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + return +end subroutine th3d_p4 + +subroutine th3d_p5(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA0_17term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + return + +end subroutine th3d_p5 + +subroutine th3d_p6(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA2_17term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + return +end subroutine th3d_p6 + +subroutine th3d_p7(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA0_12term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + + return +end subroutine th3d_p7 + +subroutine th3d_p8(temp,saln,th3d,no,mo,thbase) + implicit none + + integer :: no,mo + real, dimension(no,mo) :: temp,saln,th3d + real :: thbase + + ! --- calculate density using the equation of state. + + integer :: i,j + + include './include/stmt_fns_SIGMA2_12term.h90' + + do j= 1,mo + do i= 1,no + th3d(i,j) = sig(r8(temp(i,j)),r8(saln(i,j))) - thbase + enddo !i + enddo !j + + return +end subroutine th3d_p8 diff --git a/hycom/MSCPROGS/src/Average/hycom_mean b/hycom/MSCPROGS/src/Average/hycom_mean new file mode 100755 index 00000000..74c10cd6 Binary files /dev/null and b/hycom/MSCPROGS/src/Average/hycom_mean differ diff --git a/hycom/MSCPROGS/src/Average/hycom_mean.F90 b/hycom/MSCPROGS/src/Average/hycom_mean.F90 new file mode 100644 index 00000000..ee5e1445 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/hycom_mean.F90 @@ -0,0 +1,135 @@ +program hycom_mean + ! + ! F90 version of hycom_All/hycom_mean + ! + use mod_mean ! HYCOM mean array interface + use mod_za ! HYCOM array I/O interface + use mod_blkin + implicit none + + ! --- Form the mean (or mean-squared) of a sequence of HYCOM archive files. + ! --- Layered means weighted by layer thickness. + + character(len=81) :: label + character(len=18) :: text + character(len=80) :: flnm + !logical :: meansq,single,trcout,icegln,hisurf + integer :: mntype,iweight,i,j,n + integer :: narchs,iarch + integer :: iexpt,jexpt,kkin,yrflag,kpalet,mxlflg + real :: thbase + real(kind=8) :: time_min,time_max,time_ave,time(3) + + call xcspmd() !define idm,jdm + call zaiost() + !lp=6 + + iexpt = 0 + ii = idm + jj = jdm + iorign = 1 + jorign = 1 + + ! --- 'ntracr' = number of tracers to include in the mean + ! --- optional, default is 0 + ! --- 'kk ' = number of layers involved (1 for surface only means) + + call blkini2(i,j, 'ntracr','kk ') !read ntracr or kk as integer + if (j.eq.1) then !'ntracr' + ntracr = i + call blkini(kk, 'kk ') + else !kk + ntracr = 0 + kk = i + endif + + hisurf = .true. + trcout = ntracr.gt.0 + icegln = .true. + + ! --- array allocation and initialiation + + if (ntracr_bgc_3d.gt.0 .or. ntracr_bgc_2d.gt.0) then + bgcout = .true. + endif + + call mean_alloc + + ! --- land masks. + + call getdepth('regional.depth') + call bigrid(depths) + call mean_depths + + ! --- read and sum a sequence of archive files. + + time_min = huge(time_min) + time_max = -huge(time_max) + time_ave = 0.0 + nstep_m = 0 + + ! --- 'single' = treat input mean archive as a single sample + ! --- optional, default is meansq + ! --- 'meansq' = form meansq (rather than mean) + + call blkini2(i,j,'single','meansq') !read single or meansq as integer + if (j.eq.1) then !'single' + single = i .ne. 0 !0=F + call blkinl(meansq,'meansq') + if (meansq .and. .not.single) then + write (lp,'(a)') 'error: meansq=T requires single=T' + call flush(lp) + stop + endif !error + else !'meannsq' + meansq = i .ne. 0 !0=F + single = meansq + endif + + do ! loop until input narchs==0 + ! --- 'narchs' = number of archives to read (==0 to end input) + call blkini(narchs,'narchs') + if (narchs.eq.0) then + exit + endif + do iarch= 1,narchs + read (*,'(a)') flnm + write (lp,'(2a)') ' input file: ',trim(flnm) + call getdat(flnm,time,iweight,mntype,iexpt,yrflag,thbase) + if (iweight.eq.1) then + call mean_velocity ! calculate full velocity and KE + endif + if (single .and. mntype.eq.1) then + iweight = 1 !treat mean as a standard archive + endif + if (iweight.eq.1 .and. meansq) then + call mean_addsq(iweight) + else + call mean_add(iweight) + endif + time_min = min( time_min, time(1) ) + time_max = max( time_max, time(2) ) + time_ave = time_ave + time(3) * iweight + nstep_m = nstep_m + nstep * iweight + enddo + enddo + + ! --- output the mean or mean square + + call mean_end + time_ave = time_ave/nmean + nstep_m = nstep_m/nmean + + read (*,'(a)') flnm + write (lp,'(2a)') 'output file: ',trim(flnm) + call flush(lp) + if (meansq) then + mntype = 2 + else + mntype = 1 + endif + jexpt=iexpt + call putdat(flnm,time_min,time_max,time_ave,mntype,iexpt,jexpt,yrflag,thbase) + +end program hycom_mean + diff --git a/hycom/MSCPROGS/src/Average/include/README.stmt_fns b/hycom/MSCPROGS/src/Average/include/README.stmt_fns new file mode 100644 index 00000000..58fdc950 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/README.stmt_fns @@ -0,0 +1,27 @@ +Starting with 2.2, the HYCOM archive files no longer include density. This +is because T&S define density via the equation of state. However, this means +that the archive file must include information about the actual equation of +state. So "k dens" for montg1 (the first record in the archive file) are +actually sigver and thbase, where sigver is a small integer identifying the +equation of state (taken from a parameter defined in stmt_fns.h). + +So far we have (odd for sigma-0, even for sigma-2, >=40 for sigma4): + +sigver = 1 --- 7-term sigma-0 +sigver = 2 --- 7-term sigma-2 +sigver = 3 --- 9-term sigma-0 +sigver = 4 --- 9-term sigma-2 +sigver = 5 --- 17-term sigma-0 +sigver = 6 --- 17-term sigma-2 +sigver = 46 --- 17-term sigma-4 +sigver = 7 --- 12-term sigma-0 +sigver = 8 --- 12-term sigma-2 +sigver = 48 --- 12-term sigma-4 + +The 17-term and 12-term versions are rational function approximations, and +are equal to the associated 25-term and 18-term sigloc evaluated at the +fixed reference pressure (0 or 2000 or 4000 dbar). Note that the 17-term +version does not have a closed form expression for tofsig or sofsig, and +so can't be used if potential density is a prognostic variable. In +diagnostic programs, the 17-term tofsig and sofsig are calculated using +Newton iterations from a 12-term 1st guess. diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_12term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_12term.h90 new file mode 100644 index 00000000..ff74daa6 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_12term.h90 @@ -0,0 +1,155 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver= 7 !12-term sigma-0 +!sig2& sigver= 8 !12-term sigma-2 +!sig4& sigver=48 !12-term sigma-4 + + real(kind=8) :: sig,dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sig_n,sig_d,sig_q, dsigdt_n,dsigdt_d, dsigds_n,dsigds_d + real(kind=8) :: tofsig_a,tofsig_b,tofsig_c + real(kind=8) :: sofsig_a,sofsig_b,sofsig_c + real(kind=8) :: sigloc_n,sigloc_d,sigloc_q,dsiglocdt_n,dsiglocdt_d, dsiglocds_n,dsiglocds_d + real(kind=8) :: r,s,t,pdb,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + aone =1.0, & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- REFERENCE? + +! --- coefficients for 18-term rational function sigloc(). + real(kind=8), parameter :: & + c001=-1.4627567840659594d-01, & !num. constant coefficent + c002= 6.4247392832635697d-02, & !num. T coefficent + c003= 8.1213979591704621d-01, & !num. S coefficent + c004=-8.1321489441909698d-03, & !num. T^2 coefficent + c005= 4.5199845091090296d-03, & !num. T S coefficent + c006= 4.6347888132781394d-04, & !num. S^2 coefficent + c007= 5.0879498675039621d-03, & !num. P coefficent + c008= 1.6333913018305079d-05, & !num. P T coefficent + c009= 4.3899924880543972d-06 !num. P S coefficent + real(kind=8), parameter :: & + c011= 1.0000000000000000d+00, & !den. constant coefficent + c012= 1.0316374535350838d-02, & !den. T coefficent + c013= 8.9521792365142522d-04, & !den. S coefficent + c014=-2.8438341552142710d-05, & !den. T^2 coefficent + c015=-1.1887778959461776d-05, & !den. T S coefficent + c016=-4.0163964812921489d-06, & !den. S^2 coefficent + c017= 1.1995545126831476d-05, & !den. P coefficent + c018= 5.5234008384648383d-08, & !den. P T coefficent + c019= 8.4310335919950873d-09 !den. P S coefficent + real(kind=8), parameter :: & + c004x2=c004*2.d0, & !for dsigdt and dsiglocdt + c014x2=c014*2.d0, & !for dsigdt and dsiglocdt + c006x2=c006*2.d0, & !for dsigds and dsiglocds + c016x2=c016*2.d0 !for dsigds and dsiglocds + real(kind=8), parameter :: & + sqrmin=0.d0, & !sqrt arg can't be negative + sofmin=0.d0 !salinity can't be negative +! --- reference pressure. + real(kind=8), parameter :: prs2pdb=1.d-4 !Pascals to dbar + real(kind=8), parameter :: rpdb=0.d0 !reference pressure in dbar, sigma-0 +!sig2 real(kind=8), parameter :: rpdb=2000.d0 !reference pressure in dbar, sigma-2 +!sig4 real(kind=8), parameter :: rpdb=4000.d0 !reference pressure in dbar, sigma-4 +! --- coefficients for 12-term rational function sig() at rpdb. + real(kind=8), parameter :: & + c101=c001+rpdb*c007, & !num. constant coefficent + c102=c002+rpdb*c008, & !num. T coefficent + c103=c003+rpdb*c009 !num. S coefficent + real(kind=8), parameter :: & + c111=c011+rpdb*c017, & !num. constant coefficent + c112=c012+rpdb*c018, & !num. T coefficent + c113=c013+rpdb*c019 !num. S coefficent + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma at rpdb (dbar) as a function of temp (deg c) and salinity (psu) + + sig_n(t,s) = c101+(c102+c004*t+c005*s)*t + & + (c103+ c006*s)*s + sig_d(t,s) = c111+(c112+c014*t+c015*s)*t + & + (c113 +c016*s)*s + sig_q(t,s) = aone/sig_d(t,s) + sig( t,s) = sig_n(t,s)*sig_q(t,s) + +! --- d(sig)/dt + dsigdt_n(t,s) = c102+c004x2*t+c005*s + dsigdt_d(t,s) = c112+c014x2*t+c015*s + dsigdt( t,s) = (dsigdt_n(t,s)- & + dsigdt_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- d(sig)/ds + dsigds_n(t,s) = c103+c005*t+c006x2*s + dsigds_d(t,s) = c113+c015*t+c016x2*s + dsigds( t,s) = (dsigds_n(t,s)- & + dsigds_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a quadrati! polynominal root of a*t**2+b*t+c=0 + tofsig_a(r,s)=( c004 - & + r* c014 ) !quadratic coefficient + tofsig_b(r,s)=( (c102+ c005*s) - & + r*(c112+ c015*s) ) !linear coefficient + tofsig_c(r,s)=( (c101+(c103+c006*s)*s) - & + r*(c111+(c113+c016*s)*s) ) !constant coefficient + tofsig(r,s)=( -tofsig_b(r,s) & + -sqrt(max(sqrmin, & + tofsig_b(r,s)**2 - & + 4.0*tofsig_a(r,s)*tofsig_c(r,s))) ) / & + (2.0*tofsig_a(r,s)) + +! --- salinity (psu) as a function of sigma and temperature (deg c) +! --- find a quadratic polynominal root of a*s**2+b*s+c=0 + sofsig_a(r,t)=( c006 - & + r* c016 ) !quadratic coefficient + sofsig_b(r,t)=( (c103+ c005*t) - & + r*(c113+ c015*t) ) !linear coefficient + sofsig_c(r,t)=( (c101+(c102+c004*t)*t) - & + r*(c111+(c112+c014*t)*t) ) !constant coefficient + sofsig(r,s)=max(sofmin, & + ( -sofsig_b(r,s) & + +sqrt(max(sqrmin, & + sofsig_b(r,s)**2 - & + 4.0*sofsig_a(r,s)*sofsig_c(r,s))) ) / & + (2.0*sofsig_a(r,s)) ) + +! --- locally referenced sigma, using the 18-term equation of state. +! --- t: potential temperature; s: psu; prs: pressure + + sigloc_n(t,s,pdb) = c001+(c002+c004*t+c005*s)*t + & + (c003+ c006*s)*s + & + (c007+c008*t+c009*s)*pdb + sigloc_d(t,s,pdb) = c011+(c012+c014*t+c015*s)*t + & + (c013 +c016*s)*s + & + (c017+c018*t+c019*s)*pdb + sigloc_q(t,s,pdb) = aone/sigloc_d(t,s,pdb) + sigloc( t,s,prs)=sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/dt + dsiglocdt_n(t,s,pdb) = c002+c004x2*t+c005*s+c008*pdb + dsiglocdt_d(t,s,pdb) = c012+c014x2*t+c015*s+c018*pdb + dsiglocdt( t,s,prs)=(dsiglocdt_n(t,s,prs*prs2pdb)- & + dsiglocdt_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/ds + dsiglocds_n(t,s,pdb) = c003+c005*t+c006x2*s+c009*pdb + dsiglocds_d(t,s,pdb) = c013+c015*t+c016x2*s+c019*pdb + dsiglocds( t,s,prs)=(dsiglocds_n(t,s,prs*prs2pdb)- & + dsiglocds_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_17term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_17term.h90 new file mode 100644 index 00000000..ea406c3b --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_17term.h90 @@ -0,0 +1,208 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver= 5 !17-term sigma-0 +!sig2& sigver= 6 !17-term sigma-2 +!sig4& sigver=46 !17-term sigma-4 + + real(kind=8) :: sig,dsigdt,dsigds,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sig_n,sig_d,sig_q, dsigdt_n,dsigdt_d, dsigds_n,dsigds_d + real(kind=8) :: sigloc_n,sigloc_d,sigloc_q,dsiglocdt_n,dsiglocdt_d, dsiglocds_n,dsiglocds_d + real(kind=8) :: r,s,t,pdb,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + aone =1.0, & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- Jackett, McDougall, Feistel, Wright and Griffies (2006), +! --- Algorithms for Density, Potential Temperature, Conservative +! --- Temperature, and the Freezing Temperature of Seawater, JAOT + +! --- coefficients for 25-term rational function sigloc(). + real(kind=8), parameter :: & + c001= 9.9984085444849347d+02, & !num. constant coefficent + c002= 7.3471625860981584d+00, & !num. T coefficent + c003=-5.3211231792841769d-02, & !num. T^2 coefficent + c004= 3.6492439109814549d-04, & !num. T^3 coefficent + c005= 2.5880571023991390d+00, & !num. S coefficent + c006= 6.7168282786692355d-03, & !num. T S coefficent + c007= 1.9203202055760151d-03, & !num. S^2 coefficent + c008= 1.0000000000000000d+00, & !den. constant coefficent + c009= 7.2815210113327091d-03, & !den. T coefficent + c010=-4.4787265461983921d-05, & !den. T^2 coefficent + c011= 3.3851002965802430d-07, & !den. T^3 coefficent + c012= 1.3651202389758572d-10, & !den. T^4 coefficent + c013= 1.7632126669040377d-03, & !den. S coefficent + c014= 8.8066583251206474d-06, & !den. T S coefficent + c015= 1.8832689434804897d-10, & !den. T^3S coefficent + c016= 5.7463776745432097d-06, & !den. T S^1.5 coefficent + c017= 1.4716275472242334d-09 !den. T^3S^1.5 coefficent + real(kind=8), parameter :: & + c018= 1.1798263740430364d-02, & !num. P coefficent + c019= 9.8920219266399117d-08, & !num. P T^2 coefficent + c020= 4.6996642771754730d-06, & !num. P S coefficent + c021= 2.5862187075154352d-08, & !num. P^2 coefficent + c022= 3.2921414007960662d-12, & !num. P^2T^2 coefficent + c023= 6.7103246285651894d-06, & !den. P coefficent + c024= 2.4461698007024582d-17, & !den. P^2T^3 coefficent + c025= 9.1534417604289062d-18 !den. P^3T coefficent +! --- additional coefficients for dsiglocdt(). + real(kind=8), parameter :: & + c031= 7.3471625860981580d+00, & !num. constant coefficent + c032=-1.0642246358568354d-01, & !num. T coefficent + c033= 1.0947731732944364d-03, & !num. T^2 coefficent + c034= 6.7168282786692355d-03, & !num. S coefficent + c035= 7.2815210113327090d-03, & !den. constant coefficent + c036=-8.9574530923967840d-05, & !den. T coefficent + c037= 1.0155300889740728d-06, & !den. T^2 coefficent + c038= 5.4604809559034290d-10, & !den. T^3 coefficent + c039=-8.8066583251206470d-06, & !den. S coefficent + c040= 5.6498068304414700d-10, & !den. T^2S coefficent + c041= 2.9432550944484670d-09, & !den. T S^1.5 coefficent + c042= 1.9784043853279823d-07, & !num. P T coefficent + c043= 6.5842828015921320d-12, & !num. P^2T coefficent + c044= 7.3385094021073750d-17, & !den. P^2T^2 coefficent + c045= 9.1534417604289060d-18 !den. P^3 coefficent +! --- additional coefficients for dsiglocds(). + real(kind=8), parameter :: & + c051= 2.5880571023991390d+00, & !num. constant coefficent + c052= 6.7168282786692355d-03, & !num. T coefficent + c053= 3.8406404111520300d-03, & !num. S coefficent + c054= 1.7632126669040377d-03, & !den. constant coefficent + c055=-8.8066583251206470d-06, & !den. T coefficent + c056= 1.8832689434804897d-10, & !den. T^3 coefficent + c057= 8.6195665118148150d-06, & !den. S^0.5 coefficent + c058= 2.2074413208363504d-09, & !den. T^2S^0.5 coefficent + c059= 4.6996642771754730d-06 !num. P coefficent + + real(kind=8), parameter :: sqrmin=0.d0 !sqrt arg can't be negative +! --- reference pressure. + real(kind=8), parameter :: prs2pdb=1.d-4 !Pascals to dbar + real(kind=8), parameter :: rpdb=0.d0 !reference pressure in dbar, sigma0 +!sig2 real(kind=8), parameter :: rpdb=2000.d0 !reference pressure in dbar, sigma2 +!sig4 real(kind=8), parameter :: rpdb=4000.d0 !reference pressure in dbar, sigma4 +! --- coefficients for 17-term rational function sig() at rpdb. + real(kind=8), parameter :: & + c101=c001+(c018-c021*rpdb)*rpdb, & !num. constant coefficent + c103=c003+(c019-c022*rpdb)*rpdb, & !num. T^2 coefficent + c105=c005+c020*rpdb, & !num. S coefficent + c108=c008+c023*rpdb, & !den. constant coefficent + c109=c009-c025*rpdb**3, & !den. T coefficent + c111=c011-c024*rpdb**2 !den. T^3 coefficent +! --- additional coefficients for dsigdt(). + real(kind=8), parameter :: & + c132=c032+(c042-c043*rpdb)*rpdb, & !num. T coefficent + c135=c035-c045*rpdb**3, & !den. constant coefficent + c137=c037-c044*rpdb**2 !den. T^2 coefficent +! --- additional coefficients for dsigds(). + real(kind=8), parameter :: & + c151=c051+c059*rpdb !num. constant coefficent + real(kind=8), parameter :: & + rhoref=1.d3 !rhoref=qthref kg/m^3 + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma at rpdb (dbar) as a function of pot.temp (deg c) and salinity (psu) + + sig_n(t,s) = c101 + t*(c002+t*(c103+t*c004)) + & + s*(c105-t*c006+s*c007) + sig_d(t,s) = c108 + t*(c109+t*(c010+t*(c111+t*c012))) + & + s*(c013-t*(c014+t*t*c015) + & + sqrt(max(sqrmin,s))*(c016+t*t*c017)) + sig_q(t,s) = aone/sig_d(t,s) + sig( t,s) = sig_n(t,s)*sig_q(t,s) - rhoref + +! --- d(sig)/dt + dsigdt_n(t,s) = c031 + t*(c132+t*c033) - & + s* c034 + dsigdt_d(t,s) = c135 + t*(c036+t*(c137+t*c038)) + & + s*(c039-t*t*c040+ & + sqrt(max(sqrmin,s))*t*c041 ) + dsigdt( t,s) = (dsigdt_n(t,s)- & + dsigdt_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- d(sig)/ds +! --- additional coefficients for dsigds(). + dsigds_n(t,s) = c151 - t*c052 + s*c053 + dsigds_d(t,s) = c054 + t*(c055-t*t*c056) + & + sqrt(max(sqrmin,s))*(c057+t*t*c058) + dsigds( t,s) = (dsigds_n(t,s)- & + dsigds_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- locally referenced sigma, using the 25-term equation of state. +! --- t: potential temperature (degC); s: salinity (psu); prs: pressure (dbar) + sigloc_n(t,s,pdb) = c001 + & + t*(c002 + & + t*(c003 + & + t* c004 )) + & + s*(c005 - & + t* c006 + & + s* c007 ) + & + pdb*(c018 + & + t*t* c019 + & + s* c020 - & + pdb*(c021 + & + t*t* c022 )) + sigloc_d(t,s,pdb) = c008 + & + t*(c009 + & + t*(c010 + & + t*(c011 + & + t* c012 ))) + & + s*(c013 - & + t*(c014 + & + t*t* c015 ) + & + sqrt(max(sqrmin,s))*(c016 + & + t*t* c017 )) + & + pdb*(c023 - & + pdb*t*(t*t* c024 + & + pdb* c025 )) + sigloc_q(t,s,pdb) = aone/sigloc_d(t,s,pdb) + sigloc(t,s,prs)=sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) - rhoref + +! --- d(sig)/dt + dsiglocdt_n(t,s,pdb) = c031 + & + t*(c032 + & + t* c033 ) - & + s* c034 + & + pdb*t*(c042 - & + pdb* c043 ) + dsiglocdt_d(t,s,pdb) = c035 + & + t*(c036 + & + t*(c037 + & + t* c038 )) + & + s*(c039 - & + t*t* c040 + & + sqrt(max(sqrmin,s))*t* c041 ) - & + pdb*pdb*(t*t* c044 + & + pdb* c045 ) + dsiglocdt(t,s,prs)=(dsiglocdt_n(t,s,prs*prs2pdb)- & + dsiglocdt_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/ds + dsiglocds_n(t,s,pdb) = c051 - & + t* c052 + & + s* c053 + & + pdb* c059 + dsiglocds_d(t,s,pdb) = c054 + & + t*(c055 - & + t*t* c056 ) + & + sqrt(max(sqrmin,s))*(c057 + & + t*t* c058 ) + dsiglocds(t,s,prs)=(dsiglocds_n(t,s,prs*prs2pdb)- & + dsiglocds_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_7term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_7term.h90 new file mode 100644 index 00000000..4e5a5446 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_7term.h90 @@ -0,0 +1,113 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver=1 !7-term sigma-0 +!sig2& sigver=2 !7-term sigma-2 + + real(kind=8) :: sig,dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: a0,a1,a2,cubr,cubq,cuban,cubrl,cubim + real(kind=8) :: c1l,c2l,c3l,c4l,c5l,c6l,c7l + real(kind=8) :: r,s,t,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- coefficients for sigma-0 (based on Brydon & Sun fit) + real(kind=8), parameter :: & + c1=-1.36471E-01, & !const. coefficent + c2= 4.68181E-02, & !T coefficent + c3= 8.07004E-01, & ! S coefficent + c4=-7.45353E-03, & !T^2 coefficent + c5=-2.94418E-03, & !T S coefficent + c6= 3.43570E-05, & !T^3 coefficent + rc6= 1.0/c6, & + c7= 3.48658E-05 !T^2S coefficent + +! --- coefficients for sigma-2 (based on Brydon & Sun fit) +!sig2 real(kind=8), parameter :: +!sig2& c1= 9.77093E+00, !const. coefficent +!sig2& c2=-2.26493E-02, !T coefficent +!sig2& c3= 7.89879E-01, ! S coefficent +!sig2& c4=-6.43205E-03, !T^2 coefficent +!sig2& c5=-2.62983E-03, !T S coefficent +!sig2& c6= 2.75835E-05, !T^3 coefficent +!sig2& rc6= 1.0/c6, +!sig2& c7= 3.15235E-05 !T^2S coefficent + +! --- HYCOM pressure to bar, for locally referenced equations + real(kind=8), parameter :: prs2pb=1.e-5 !Pascals to bar + +! --- sub-coefficients for locally referenced sigma +! --- a fit towards Jackett & McDougall (1995) + real(kind=8), parameter, dimension(7) :: & + alphap = (/ -0.1364705627213484 , 0.04681812123458564, & + 0.80700383913187 ,-0.007453530323180844, & + -0.002944183249153631 , 0.00003435702568990446, & + 0.0000348657661057688 /) & + ,betap = (/ 0.05064226654169138 ,-0.0003571087848996894, & + -0.0000876148051892879, 5.252431910751829e-6, & + 1.579762259448864e-6 ,-3.466867400295792e-8, & + -1.687643078774232e-8 /) & + ,gammap = (/ -5.526396144304812e-6 , 4.885838128243163e-8, & + 9.96026931578033e-9 ,-7.251389796582352e-10, & + -3.987360250058777e-11, 4.006307891935698e-12, & + 8.26367520608008e-13 /) + +! --- auxiliary statements for finding root of cubic polynomial + a0(s,r)=(c1+c3*s-r)*rc6 !constant coefficient + a1(s) =(c2+c5*s )*rc6 !linear coefficient + a2(s) =(c4+c7*s )*rc6 !quadratic coefficient + !cubic coefficient is c6*rc6=1.0 + cubq(s)=a3rd*a1(s)-(a3rd*a2(s))**2 + cubr(r,s)=a3rd*(0.5d0*a1(s)*a2(s)-1.5d0*a0(s,r))-(a3rd*a2(s))**3 +! --- if q**3+r**2>0, water is too dense to yield real root at given +! --- salinitiy. setting q**3+r**2=0 in that case is equivalent to +! --- lowering sigma until a double real root is obtained. + cuban(r,s)=a3rd*atan2(sqrt(max(0.d0,-(cubq(s)**3+cubr(r,s)**2))), & + cubr(r,s)) + cubrl(r,s)=sqrt(-cubq(s))*cos(cuban(r,s)) + cubim(r,s)=sqrt(-cubq(s))*sin(cuban(r,s)) + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma-theta as a function of temp (deg c) and salinity (psu) +! --- (friedrich-levitus, polynomial fit that is cubic in T and linear in S) + + sig(t,s)=(c1+c3*s+t*(c2+c5*s+t*(c4+c7*s+c6*t))) + +! --- d(sig)/dt + dsigdt(t,s)=(c2+c5*s+2.0*t*(c4+c7*s+1.5*c6*t)) + +! --- d(sig)/ds + dsigds(t,s)=(c3+t*(c5+t*c7)) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a cubic polynominal root of t**3+a2*t**2+a1*t+a0=0 + tofsig(r,s)=-cubrl(r,s)+sqrt(3.0)*cubim(r,s)-a3rd*a2(s) + +! --- salinity (psu) as a function of sigma and temperature (deg c) + sofsig(r,t)=(r-c1-t*(c2+t*(c4+c6*t)))/(c3+t*(c5+c7*t)) + +! --- locally referenced sigma, a fit towards Jackett & McDougall (1995) +! --- t: potential temperature; s: psu; prs: pressure + c1l(prs)=alphap(1)+prs2pb*prs*(betap(1)+prs2pb*prs*gammap(1)) + c2l(prs)=alphap(2)+prs2pb*prs*(betap(2)+prs2pb*prs*gammap(2)) + c3l(prs)=alphap(3)+prs2pb*prs*(betap(3)+prs2pb*prs*gammap(3)) + c4l(prs)=alphap(4)+prs2pb*prs*(betap(4)+prs2pb*prs*gammap(4)) + c5l(prs)=alphap(5)+prs2pb*prs*(betap(5)+prs2pb*prs*gammap(5)) + c6l(prs)=alphap(6)+prs2pb*prs*(betap(6)+prs2pb*prs*gammap(6)) + c7l(prs)=alphap(7)+prs2pb*prs*(betap(7)+prs2pb*prs*gammap(7)) + sigloc(t,s,prs)=c1l(prs)+c3l(prs)*s+ & + t*(c2l(prs)+c5l(prs)*s+t*(c4l(prs)+c7l(prs)*s+c6l(prs)*t)) + dsiglocdt(t,s,prs)=(c2l(prs)+c5l(prs)*s+ & + 2.0*t*(c4l(prs)+c7l(prs)*s+1.5*c6l(prs)*t)) + dsiglocds(t,s,prs)=(c3l(prs)+t*(c5l(prs)+t*c7l(prs))) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_9term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_9term.h90 new file mode 100644 index 00000000..26923ff9 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA0_9term.h90 @@ -0,0 +1,129 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver=3 !9-term sigma-0 +!sig2& sigver=4 !9-term sigma-2 + + real(kind=8) :: sig,dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sofsig_a,sofsig_b,sofsig_c + real(kind=8) :: a0,a1,a2,cubr,cubq,cuban,cubrl,cubim + real(kind=8) :: c1l,c2l,c3l,c4l,c5l,c6l,c7l + real(kind=8) :: r,s,t,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- sigma-theta as a function of temp (deg c) and salinity (psu) +! --- (9-term polynomial fit to T:[-2:30],S:[18:38]) + +! --- coefficients for sigma-0. + real(kind=8), parameter :: & + c1=-4.311829E-02, & !const. coefficent + c2= 5.429948E-02, & !T coefficent + c3= 8.011774E-01, & ! S coefficent + c4=-7.641336E-03, & !T^2 coefficent + c5=-3.258442E-03, & !T S coefficent + c6= 3.757643E-05, & !T^3 coefficent + rc6=1.0/c6, & + c7= 3.630361E-05, & !T^2S coefficent + c8= 8.675546E-05, & ! S^2 coefficent + c9= 3.995086E-06 !T S^2 coefficent +! --- coefficients for sigma-2. +!sig2 real(kind=8), parameter :: +!sig2& c1= 9.903308E+00, !const. coefficent +!sig2& c2=-1.618075E-02, !T coefficent +!sig2& c3= 7.819166E-01, ! S coefficent +!sig2& c4=-6.593939E-03, !T^2 coefficent +!sig2& c5=-2.896464E-03, !T S coefficent +!sig2& c6= 3.038697E-05, !T^3 coefficent +!sig2& rc6= 1.0/c6, +!sig2& c7= 3.266933E-05, !T^2S coefficent +!sig2& c8= 1.180109E-04, ! S^2 coefficent +!sig2& c9= 3.399511E-06 !T S^2 coefficent + +! --- HYCOM pressure to bar, for locally referenced equations + real(kind=8), parameter :: prs2pb=1.e-5 !Pascals to bar + +! --- sub-coefficients for locally referenced sigma +! --- a fit towards Jackett & McDougall (1995) + real(kind=8), parameter, dimension(7) :: & + alphap = (/ -0.1364705627213484 , 0.04681812123458564, & + 0.80700383913187 ,-0.007453530323180844, & + -0.002944183249153631 , 0.00003435702568990446, & + 0.0000348657661057688 /) & + ,betap = (/ 0.05064226654169138 ,-0.0003571087848996894, & + -0.0000876148051892879, 5.252431910751829e-6, & + 1.579762259448864e-6 ,-3.466867400295792e-8, & + -1.687643078774232e-8 /) & + ,gammap = (/ -5.526396144304812e-6 , 4.885838128243163e-8, & + 9.96026931578033e-9 ,-7.251389796582352e-10, & + -3.987360250058777e-11, 4.006307891935698e-12, & + 8.26367520608008e-13 /) + +! --- auxiliary statements for finding root of cubic polynomial + a0(s,r)=(c1+s*(c3+s*c8)-r)*rc6 !constant coefficient + a1(s) =(c2+s*(c5+s*c9) )*rc6 !linear coefficient + a2(s) =(c4+s* c7 )*rc6 !quadratic coefficient + !cubic coefficient is c6*rc6=1.0 + cubq(s) =a3rd* a1(s) -(a3rd*a2(s))**2 + cubr(r,s)=a3rd*(0.5*a1(s)*a2(s)-1.5*a0(s,r))-(a3rd*a2(s))**3 +! --- if q**3+r**2>0, water is too dense to yield real root at given +! --- salinitiy. setting q**3+r**2=0 in that case is equivalent to +! --- lowering sigma until a double real root is obtained. + cuban(r,s)=a3rd*atan2(sqrt(max(0.d0,-(cubq(s)**3+cubr(r,s)**2))), & + cubr(r,s)) + cubrl(r,s)=sqrt(-cubq(s))*cos(cuban(r,s)) + cubim(r,s)=sqrt(-cubq(s))*sin(cuban(r,s)) + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma-theta as a function of temp (deg c) and salinity (psu) +! --- (polynomial fit that is cubic in T and quadratic in S) + + sig(t,s)=(c1+s*(c3+s* c8)+ & + t*(c2+s*(c5+s*c9)+t*(c4+s*c7+t*c6))) + +! --- d(sig)/dt + dsigdt(t,s)=(c2+s*(c5+s*c9)+2.0*t*(c4+s*c7+1.5*t*c6)) + +! --- d(sig)/ds + dsigds(t,s)=(c3+t*(c5+t*c7)+2.0*s* c8) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a cubic polynominal root of t**3+a2*t**2+a1*t+a0=0 + tofsig(r,s)=-cubrl(r,s)+sqrt(3.0)*cubim(r,s)-a3rd*a2(s) + +! --- salinity (psu) as a function of sigma and temperature (deg c) +! --- find a quadratic polynominal root of a*s**2+b*s+c=0 + sofsig_a(r,t)=(c8+t* c9) !quadratic coefficient + sofsig_b(r,t)=(c3+t*(c5+t* c7)) !linear coefficient + sofsig_c(r,t)=(c1+t*(c2+t*(c4+t*c6))-r) !constant coefficient + sofsig(r,t)=(2.d0*sofsig_c(r,t))/ & + (-sofsig_b(r,t) & + -sign(sqrt(sofsig_b(r,t)**2- & + 4.d0*sofsig_a(r,t)*sofsig_c(r,t)), & + sofsig_b(r,t))) + +! --- locally referenced sigma, a fit towards Jackett & McDougall (1995) +! --- t: potential temperature; s: psu; prs: pressure + c1l(prs)=alphap(1)+prs2pb*prs*(betap(1)+prs2pb*prs*gammap(1)) + c2l(prs)=alphap(2)+prs2pb*prs*(betap(2)+prs2pb*prs*gammap(2)) + c3l(prs)=alphap(3)+prs2pb*prs*(betap(3)+prs2pb*prs*gammap(3)) + c4l(prs)=alphap(4)+prs2pb*prs*(betap(4)+prs2pb*prs*gammap(4)) + c5l(prs)=alphap(5)+prs2pb*prs*(betap(5)+prs2pb*prs*gammap(5)) + c6l(prs)=alphap(6)+prs2pb*prs*(betap(6)+prs2pb*prs*gammap(6)) + c7l(prs)=alphap(7)+prs2pb*prs*(betap(7)+prs2pb*prs*gammap(7)) + sigloc(t,s,prs)=c1l(prs)+c3l(prs)*s+ & + t*(c2l(prs)+c5l(prs)*s+t*(c4l(prs)+c7l(prs)*s+c6l(prs)*t)) + dsiglocdt(t,s,prs)=(c2l(prs)+c5l(prs)*s+ & + 2.0*t*(c4l(prs)+c7l(prs)*s+1.5*c6l(prs)*t)) + dsiglocds(t,s,prs)=(c3l(prs)+t*(c5l(prs)+t*c7l(prs))) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_12term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_12term.h90 new file mode 100644 index 00000000..12998503 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_12term.h90 @@ -0,0 +1,156 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver= 8 !12-term sigma-2 +!sig0& sigver= 7 !12-term sigma-0 +!sig4& sigver=48 !12-term sigma-4 + + real(kind=8) :: sig + real(kind=8) :: dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sig_n,sig_d,sig_q, dsigdt_n,dsigdt_d, dsigds_n,dsigds_d + real(kind=8) :: tofsig_a,tofsig_b,tofsig_c + real(kind=8) :: sofsig_a,sofsig_b,sofsig_c + real(kind=8) :: sigloc_n,sigloc_d,sigloc_q,dsiglocdt_n,dsiglocdt_d, dsiglocds_n,dsiglocds_d + real(kind=8) :: r,s,t,pdb,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + aone =1.0, & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- REFERENCE? + +! --- coefficients for 18-term rational function sigloc(). + real(kind=8), parameter :: & + c001=-1.4627567840659594d-01, & !num. constant coefficent + c002= 6.4247392832635697d-02, & !num. T coefficent + c003= 8.1213979591704621d-01, & !num. S coefficent + c004=-8.1321489441909698d-03, & !num. T^2 coefficent + c005= 4.5199845091090296d-03, & !num. T S coefficent + c006= 4.6347888132781394d-04, & !num. S^2 coefficent + c007= 5.0879498675039621d-03, & !num. P coefficent + c008= 1.6333913018305079d-05, & !num. P T coefficent + c009= 4.3899924880543972d-06 !num. P S coefficent + real(kind=8), parameter :: & + c011= 1.0000000000000000d+00, & !den. constant coefficent + c012= 1.0316374535350838d-02, & !den. T coefficent + c013= 8.9521792365142522d-04, & !den. S coefficent + c014=-2.8438341552142710d-05, & !den. T^2 coefficent + c015=-1.1887778959461776d-05, & !den. T S coefficent + c016=-4.0163964812921489d-06, & !den. S^2 coefficent + c017= 1.1995545126831476d-05, & !den. P coefficent + c018= 5.5234008384648383d-08, & !den. P T coefficent + c019= 8.4310335919950873d-09 !den. P S coefficent + real(kind=8), parameter :: & + c004x2=c004*2.d0, & !for dsigdt and dsiglocdt + c014x2=c014*2.d0, & !for dsigdt and dsiglocdt + c006x2=c006*2.d0, & !for dsigds and dsiglocds + c016x2=c016*2.d0 !for dsigds and dsiglocds + real(kind=8), parameter :: & + sqrmin=0.d0, & !sqrt arg can't be negative + sofmin=0.d0 !salinity can't be negative +! --- reference pressure. + real(kind=8), parameter :: prs2pdb=1.d-4 !Pascals to dbar +!sig0 real(kind=8), parameter :: rpdb=0.d0 !reference pressure in dbar, sigma-0 + real(kind=8), parameter :: rpdb=2000.d0 !reference pressure in dbar, sigma-2 +!sig4 real(kind=8), parameter :: rpdb=4000.d0 !reference pressure in dbar, sigma-4 +! --- coefficients for 12-term rational function sig() at rpdb. + real(kind=8), parameter :: & + c101=c001+rpdb*c007, & !num. constant coefficent + c102=c002+rpdb*c008, & !num. T coefficent + c103=c003+rpdb*c009 !num. S coefficent + real(kind=8), parameter :: & + c111=c011+rpdb*c017, & !num. constant coefficent + c112=c012+rpdb*c018, & !num. T coefficent + c113=c013+rpdb*c019 !num. S coefficent + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma at rpdb (dbar) as a function of temp (deg c) and salinity (psu) + + sig_n(t,s) = c101+(c102+c004*t+c005*s)*t + & + (c103+ c006*s)*s + sig_d(t,s) = c111+(c112+c014*t+c015*s)*t + & + (c113 +c016*s)*s + sig_q(t,s) = aone/sig_d(t,s) + sig( t,s) = sig_n(t,s)*sig_q(t,s) + +! --- d(sig)/dt + dsigdt_n(t,s) = c102+c004x2*t+c005*s + dsigdt_d(t,s) = c112+c014x2*t+c015*s + dsigdt( t,s) = (dsigdt_n(t,s)- & + dsigdt_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- d(sig)/ds + dsigds_n(t,s) = c103+c005*t+c006x2*s + dsigds_d(t,s) = c113+c015*t+c016x2*s + dsigds( t,s) = (dsigds_n(t,s)- & + dsigds_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a quadratic polynominal root of a*t**2+b*t+c=0 + tofsig_a(r,s)=( c004 - & + r* c014 ) !quadratic coefficient + tofsig_b(r,s)=( (c102+ c005*s) - & + r*(c112+ c015*s) ) !linear coefficient + tofsig_c(r,s)=( (c101+(c103+c006*s)*s) - & + r*(c111+(c113+c016*s)*s) ) !constant coefficient + tofsig(r,s)=( -tofsig_b(r,s) & + -sqrt(max(sqrmin, & + tofsig_b(r,s)**2 - & + 4.0*tofsig_a(r,s)*tofsig_c(r,s))) ) / & + (2.0*tofsig_a(r,s)) + +! --- salinity (psu) as a function of sigma and temperature (deg c) +! --- find a quadratic polynominal root of a*s**2+b*s+c=0 + sofsig_a(r,t)=( c006 - & + r* c016 ) !quadratic coefficient + sofsig_b(r,t)=( (c103+ c005*t) - & + r*(c113+ c015*t) ) !linear coefficient + sofsig_c(r,t)=( (c101+(c102+c004*t)*t) - & + r*(c111+(c112+c014*t)*t) ) !constant coefficient + sofsig(r,s)=max(sofmin, & + ( -sofsig_b(r,s) & + +sqrt(max(sqrmin, & + sofsig_b(r,s)**2 - & + 4.0*sofsig_a(r,s)*sofsig_c(r,s))) ) / & + (2.0*sofsig_a(r,s)) ) + +! --- locally referenced sigma, using the 18-term equation of state. +! --- t: potential temperature; s: psu; prs: pressure + + sigloc_n(t,s,pdb) = c001+(c002+c004*t+c005*s)*t + & + (c003+ c006*s)*s + & + (c007+c008*t+c009*s)*pdb + sigloc_d(t,s,pdb) = c011+(c012+c014*t+c015*s)*t + & + (c013 +c016*s)*s + & + (c017+c018*t+c019*s)*pdb + sigloc_q(t,s,pdb) = aone/sigloc_d(t,s,pdb) + sigloc( t,s,prs)=sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/dt + dsiglocdt_n(t,s,pdb) = c002+c004x2*t+c005*s+c008*pdb + dsiglocdt_d(t,s,pdb) = c012+c014x2*t+c015*s+c018*pdb + dsiglocdt( t,s,prs)=(dsiglocdt_n(t,s,prs*prs2pdb)- & + dsiglocdt_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/ds + dsiglocds_n(t,s,pdb) = c003+c005*t+c006x2*s+c009*pdb + dsiglocds_d(t,s,pdb) = c013+c015*t+c016x2*s+c019*pdb + dsiglocds( t,s,prs)=(dsiglocds_n(t,s,prs*prs2pdb)- & + dsiglocds_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_17term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_17term.h90 new file mode 100644 index 00000000..7c1d2281 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_17term.h90 @@ -0,0 +1,208 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver= 6 !17-term sigma-2 +!sig0& sigver= 5 !17-term sigma-0 +!sig4& sigver=46 !17-term sigma-4 + + real(kind=8) :: sig,dsigdt,dsigds,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sig_n,sig_d,sig_q, dsigdt_n,dsigdt_d, dsigds_n,dsigds_d + real(kind=8) :: sigloc_n,sigloc_d,sigloc_q,dsiglocdt_n,dsiglocdt_d, dsiglocds_n,dsiglocds_d + real(kind=8) :: r,s,t,pdb,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + aone =1.0, & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- Jackett, McDougall, Feistel, Wright and Griffies (2006), +! --- Algorithms for Density, Potential Temperature, Conservative +! --- Temperature, and the Freezing Temperature of Seawater, JAOT + +! --- coefficients for 25-term rational function sigloc(). + real(kind=8), parameter :: & + c001= 9.9984085444849347d+02, & !num. constant coefficent + c002= 7.3471625860981584d+00, & !num. T coefficent + c003=-5.3211231792841769d-02, & !num. T^2 coefficent + c004= 3.6492439109814549d-04, & !num. T^3 coefficent + c005= 2.5880571023991390d+00, & !num. S coefficent + c006= 6.7168282786692355d-03, & !num. T S coefficent + c007= 1.9203202055760151d-03, & !num. S^2 coefficent + c008= 1.0000000000000000d+00, & !den. constant coefficent + c009= 7.2815210113327091d-03, & !den. T coefficent + c010=-4.4787265461983921d-05, & !den. T^2 coefficent + c011= 3.3851002965802430d-07, & !den. T^3 coefficent + c012= 1.3651202389758572d-10, & !den. T^4 coefficent + c013= 1.7632126669040377d-03, & !den. S coefficent + c014= 8.8066583251206474d-06, & !den. T S coefficent + c015= 1.8832689434804897d-10, & !den. T^3S coefficent + c016= 5.7463776745432097d-06, & !den. T S^1.5 coefficent + c017= 1.4716275472242334d-09 !den. T^3S^1.5 coefficent + real(kind=8), parameter :: & + c018= 1.1798263740430364d-02, & !num. P coefficent + c019= 9.8920219266399117d-08, & !num. P T^2 coefficent + c020= 4.6996642771754730d-06, & !num. P S coefficent + c021= 2.5862187075154352d-08, & !num. P^2 coefficent + c022= 3.2921414007960662d-12, & !num. P^2T^2 coefficent + c023= 6.7103246285651894d-06, & !den. P coefficent + c024= 2.4461698007024582d-17, & !den. P^2T^3 coefficent + c025= 9.1534417604289062d-18 !den. P^3T coefficent +! --- additional coefficients for dsiglocdt(). + real(kind=8), parameter :: & + c031= 7.3471625860981580d+00, & !num. constant coefficent + c032=-1.0642246358568354d-01, & !num. T coefficent + c033= 1.0947731732944364d-03, & !num. T^2 coefficent + c034= 6.7168282786692355d-03, & !num. S coefficent + c035= 7.2815210113327090d-03, & !den. constant coefficent + c036=-8.9574530923967840d-05, & !den. T coefficent + c037= 1.0155300889740728d-06, & !den. T^2 coefficent + c038= 5.4604809559034290d-10, & !den. T^3 coefficent + c039=-8.8066583251206470d-06, & !den. S coefficent + c040= 5.6498068304414700d-10, & !den. T^2S coefficent + c041= 2.9432550944484670d-09, & !den. T S^1.5 coefficent + c042= 1.9784043853279823d-07, & !num. P T coefficent + c043= 6.5842828015921320d-12, & !num. P^2T coefficent + c044= 7.3385094021073750d-17, & !den. P^2T^2 coefficent + c045= 9.1534417604289060d-18 !den. P^3 coefficent +! --- additional coefficients for dsiglocds(). + real(kind=8), parameter :: & + c051= 2.5880571023991390d+00, & !num. constant coefficent + c052= 6.7168282786692355d-03, & !num. T coefficent + c053= 3.8406404111520300d-03, & !num. S coefficent + c054= 1.7632126669040377d-03, & !den. constant coefficent + c055=-8.8066583251206470d-06, & !den. T coefficent + c056= 1.8832689434804897d-10, & !den. T^3 coefficent + c057= 8.6195665118148150d-06, & !den. S^0.5 coefficent + c058= 2.2074413208363504d-09, & !den. T^2S^0.5 coefficent + c059= 4.6996642771754730d-06 !num. P coefficent + + real(kind=8), parameter :: sqrmin=0.d0 !sqrt arg can't be negative +! --- reference pressure. + real(kind=8), parameter :: prs2pdb=1.d-4 !Pascals to dbar +!sig0 real(kind=8), parameter :: rpdb=0.d0 !reference pressure in dbar, sigma0 + real(kind=8), parameter :: rpdb=2000.d0 !reference pressure in dbar, sigma2 +!sig4 real(kind=8), parameter :: rpdb=4000.d0 !reference pressure in dbar, sigma4 +! --- coefficients for 17-term rational function sig() at rpdb. + real(kind=8), parameter :: & + c101=c001+(c018-c021*rpdb)*rpdb, & !num. constant coefficent + c103=c003+(c019-c022*rpdb)*rpdb, & !num. T^2 coefficent + c105=c005+c020*rpdb, & !num. S coefficent + c108=c008+c023*rpdb, & !den. constant coefficent + c109=c009-c025*rpdb**3, & !den. T coefficent + c111=c011-c024*rpdb**2 !den. T^3 coefficent +! --- additional coefficients for dsigdt(). + real(kind=8), parameter :: & + c132=c032+(c042-c043*rpdb)*rpdb, & !num. T coefficent + c135=c035-c045*rpdb**3, & !den. constant coefficent + c137=c037-c044*rpdb**2 !den. T^2 coefficent +! --- additional coefficients for dsigds(). + real(kind=8), parameter :: & + c151=c051+c059*rpdb !num. constant coefficent + real, parameter :: & + rhoref=1.d3 !rhoref=qthref kg/m^3 + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma at rpdb (dbar) as a function of pot.temp (deg c) and salinity (psu) + + sig_n(t,s) = c101 + t*(c002+t*(c103+t*c004)) + & + s*(c105-t*c006+s*c007) + sig_d(t,s) = c108 + t*(c109+t*(c010+t*(c111+t*c012))) + & + s*(c013-t*(c014+t*t*c015) + & + sqrt(max(sqrmin,s))*(c016+t*t*c017)) + sig_q(t,s) = aone/sig_d(t,s) + sig( t,s) = sig_n(t,s)*sig_q(t,s) - rhoref + +! --- d(sig)/dt + dsigdt_n(t,s) = c031 + t*(c132+t*c033) - & + s* c034 + dsigdt_d(t,s) = c135 + t*(c036+t*(c137+t*c038)) + & + s*(c039-t*t*c040+ & + sqrt(max(sqrmin,s))*t*c041 ) + dsigdt( t,s) = (dsigdt_n(t,s)- & + dsigdt_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- d(sig)/ds +! --- additional coefficients for dsigds(). + dsigds_n(t,s) = c151 - t*c052 + s*c053 + dsigds_d(t,s) = c054 + t*(c055-t*t*c056) + & + sqrt(max(sqrmin,s))*(c057+t*t*c058) + dsigds( t,s) = (dsigds_n(t,s)- & + dsigds_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- locally referenced sigma, using the 25-term equation of state. +! --- t: potential temperature (degC); s: salinity (psu); prs: pressure (dbar) + sigloc_n(t,s,pdb) = c001 + & + t*(c002 + & + t*(c003 + & + t* c004 )) + & + s*(c005 - & + t* c006 + & + s* c007 ) + & + pdb*(c018 + & + t*t* c019 + & + s* c020 - & + pdb*(c021 + & + t*t* c022 )) + sigloc_d(t,s,pdb) = c008 + & + t*(c009 + & + t*(c010 + & + t*(c011 + & + t* c012 ))) + & + s*(c013 - & + t*(c014 + & + t*t* c015 ) + & + sqrt(max(sqrmin,s))*(c016 + & + t*t* c017 )) + & + pdb*(c023 - & + pdb*t*(t*t* c024 + & + pdb* c025 )) + sigloc_q(t,s,pdb) = aone/sigloc_d(t,s,pdb) + sigloc(t,s,prs)=sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) - rhoref + +! --- d(sig)/dt + dsiglocdt_n(t,s,pdb) = c031 + & + t*(c032 + & + t* c033 ) - & + s* c034 + & + pdb*t*(c042 - & + pdb* c043 ) + dsiglocdt_d(t,s,pdb) = c035 + & + t*(c036 + & + t*(c037 + & + t* c038 )) + & + s*(c039 - & + t*t* c040 + & + sqrt(max(sqrmin,s))*t* c041 ) - & + pdb*pdb*(t*t* c044 + & + pdb* c045 ) + dsiglocdt(t,s,prs)=(dsiglocdt_n(t,s,prs*prs2pdb)- & + dsiglocdt_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/ds + dsiglocds_n(t,s,pdb) = c051 - & + t* c052 + & + s* c053 + & + pdb* c059 + dsiglocds_d(t,s,pdb) = c054 + & + t*(c055 - & + t*t* c056 ) + & + sqrt(max(sqrmin,s))*(c057 + & + t*t* c058 ) + dsiglocds(t,s,prs)=(dsiglocds_n(t,s,prs*prs2pdb)- & + dsiglocds_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_7term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_7term.h90 new file mode 100644 index 00000000..45938c24 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_7term.h90 @@ -0,0 +1,116 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver=2 !7-term sigma-2 +!sig0& sigver=1 !7-term sigma-0 + + real*8 sig,dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + + real*8 r8 + real*4 r4 + + real*8 a0,a1,a2,cubr,cubq,cuban,cubrl,cubim + real c1l,c2l,c3l,c4l,c5l,c6l,c7l + + real*8 r,s,t,prs + + real*8, parameter :: & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- coefficients for sigma-0 (based on Brydon & Sun fit) +!sig0 real*8, parameter :: & +!sig0& c1=-1.36471E-01, & !const. coefficent +!sig0& c2= 4.68181E-02, & !T coefficent +!sig0& c3= 8.07004E-01, & ! S coefficent +!sig0& c4=-7.45353E-03, & !T^2 coefficent +!sig0& c5=-2.94418E-03, & !T S coefficent +!sig0& c6= 3.43570E-05, & !T^3 coefficent +!sig0& rc6= 1.0/c6, & +!sig0& c7= 3.48658E-05 !T^2S coefficent + +! --- coefficients for sigma-2 (based on Brydon & Sun fit) + real*8, parameter :: & + c1= 9.77093E+00, & !const. coefficent + c2=-2.26493E-02, & !T coefficent + c3= 7.89879E-01, & ! S coefficent + c4=-6.43205E-03, & !T^2 coefficent + c5=-2.62983E-03, & !T S coefficent + c6= 2.75835E-05, & !T^3 coefficent + rc6= 1.0/c6, & + c7= 3.15235E-05 !T^2S coefficent + +! --- HYCOM pressure to bar, for locally referenced equations + real*8, parameter :: prs2pb=1.e-5 !Pascals to bar + +! --- sub-coefficients for locally referenced sigma +! --- a fit towards Jackett & McDougall (1995) + real*8, parameter, dimension(7) :: & + alphap = (/ -0.1364705627213484 , 0.04681812123458564, & + 0.80700383913187 ,-0.007453530323180844, & + -0.002944183249153631 , 0.00003435702568990446, & + 0.0000348657661057688 /) & + ,betap = (/ 0.05064226654169138 ,-0.0003571087848996894, & + -0.0000876148051892879, 5.252431910751829e-6, & + 1.579762259448864e-6 ,-3.466867400295792e-8, & + -1.687643078774232e-8 /) & + ,gammap = (/ -5.526396144304812e-6 , 4.885838128243163e-8, & + 9.96026931578033e-9 ,-7.251389796582352e-10, & + -3.987360250058777e-11, 4.006307891935698e-12, & + 8.26367520608008e-13 /) + +! --- auxiliary statements for finding root of cubic polynomial + a0(s,r)=(c1+c3*s-r)*rc6 !constant coefficient + a1(s) =(c2+c5*s )*rc6 !linear coefficient + a2(s) =(c4+c7*s )*rc6 !quadratic coefficient + !cubic coefficient is c6*rc6=1.0 + cubq(s)=a3rd*a1(s)-(a3rd*a2(s))**2 + cubr(r,s)=a3rd*(0.5d0*a1(s)*a2(s)-1.5d0*a0(s,r))-(a3rd*a2(s))**3 +! --- if q**3+r**2>0, water is too dense to yield real root at given +! --- salinitiy. setting q**3+r**2=0 in that case is equivalent to +! --- lowering sigma until a double real root is obtained. + cuban(r,s)=a3rd*atan2(sqrt(max(0.d0,-(cubq(s)**3+cubr(r,s)**2))), & + cubr(r,s)) + cubrl(r,s)=sqrt(-cubq(s))*cos(cuban(r,s)) + cubim(r,s)=sqrt(-cubq(s))*sin(cuban(r,s)) + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma-theta as a function of temp (deg c) and salinity (psu) +! --- (friedrich-levitus, polynomial fit that is cubic in T and linear in S) + + sig(t,s)=(c1+c3*s+t*(c2+c5*s+t*(c4+c7*s+c6*t))) + +! --- d(sig)/dt + dsigdt(t,s)=(c2+c5*s+2.0*t*(c4+c7*s+1.5*c6*t)) + +! --- d(sig)/ds + dsigds(t,s)=(c3+t*(c5+t*c7)) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a cubic polynominal root of t**3+a2*t**2+a1*t+a0=0 + tofsig(r,s)=-cubrl(r,s)+sqrt(3.0)*cubim(r,s)-a3rd*a2(s) + +! --- salinity (psu) as a function of sigma and temperature (deg c) + sofsig(r,t)=(r-c1-t*(c2+t*(c4+c6*t)))/(c3+t*(c5+c7*t)) + +! --- locally referenced sigma, a fit towards Jackett & McDougall (1995) +! --- t: potential temperature; s: psu; prs: pressure + c1l(prs)=alphap(1)+prs2pb*prs*(betap(1)+prs2pb*prs*gammap(1)) + c2l(prs)=alphap(2)+prs2pb*prs*(betap(2)+prs2pb*prs*gammap(2)) + c3l(prs)=alphap(3)+prs2pb*prs*(betap(3)+prs2pb*prs*gammap(3)) + c4l(prs)=alphap(4)+prs2pb*prs*(betap(4)+prs2pb*prs*gammap(4)) + c5l(prs)=alphap(5)+prs2pb*prs*(betap(5)+prs2pb*prs*gammap(5)) + c6l(prs)=alphap(6)+prs2pb*prs*(betap(6)+prs2pb*prs*gammap(6)) + c7l(prs)=alphap(7)+prs2pb*prs*(betap(7)+prs2pb*prs*gammap(7)) + sigloc(t,s,prs)=c1l(prs)+c3l(prs)*s+ & + t*(c2l(prs)+c5l(prs)*s+t*(c4l(prs)+c7l(prs)*s+c6l(prs)*t)) + dsiglocdt(t,s,prs)=(c2l(prs)+c5l(prs)*s+ & + 2.0*t*(c4l(prs)+c7l(prs)*s+1.5*c6l(prs)*t)) + dsiglocds(t,s,prs)=(c3l(prs)+t*(c5l(prs)+t*c7l(prs))) + +! --- auxiliary statement for real to real*8 conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_9term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_9term.h90 new file mode 100644 index 00000000..7e86d518 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA2_9term.h90 @@ -0,0 +1,129 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver=4 !9-term sigma-2 +!sig0& sigver=3 !9-term sigma-0 + + real(kind=8) :: sig,dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sofsig_a,sofsig_b,sofsig_c + real(kind=8) :: a0,a1,a2,cubr,cubq,cuban,cubrl,cubim + real(kind=8) :: c1l,c2l,c3l,c4l,c5l,c6l,c7l + real(kind=8) :: r,s,t,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- sigma-theta as a function of temp (deg c) and salinity (psu) +! --- (9-term polynomial fit to T:[-2:30],S:[18:38]) + +! --- coefficients for sigma-0. +!sig0 real(kind=8), parameter :: +!sig0& c1=-4.311829E-02, !const. coefficent +!sig0& c2= 5.429948E-02, !T coefficent +!sig0& c3= 8.011774E-01, ! S coefficent +!sig0& c4=-7.641336E-03, !T^2 coefficent +!sig0& c5=-3.258442E-03, !T S coefficent +!sig0& c6= 3.757643E-05, !T^3 coefficent +!sig0& rc6=1.0/c6, +!sig0& c7= 3.630361E-05, !T^2S coefficent +!sig0& c8= 8.675546E-05, ! S^2 coefficent +!sig0& c9= 3.995086E-06 !T S^2 coefficent +! --- coefficients for sigma-2. + real(kind=8), parameter :: & + c1= 9.903308E+00, & !const. coefficent + c2=-1.618075E-02, & !T coefficent + c3= 7.819166E-01, & ! S coefficent + c4=-6.593939E-03, & !T^2 coefficent + c5=-2.896464E-03, & !T S coefficent + c6= 3.038697E-05, & !T^3 coefficent + rc6= 1.0/c6, & + c7= 3.266933E-05, & !T^2S coefficent + c8= 1.180109E-04, & ! S^2 coefficent + c9= 3.399511E-06 !T S^2 coefficent + +! --- HYCOM pressure to bar, for locally referenced equations + real(kind=8), parameter :: prs2pb=1.e-5 !Pascals to bar + +! --- sub-coefficients for locally referenced sigma +! --- a fit towards Jackett & McDougall (1995) + real(kind=8), parameter, dimension(7) :: & + alphap = (/ -0.1364705627213484 , 0.04681812123458564, & + 0.80700383913187 ,-0.007453530323180844, & + -0.002944183249153631 , 0.00003435702568990446, & + 0.0000348657661057688 /) & + ,betap = (/ 0.05064226654169138 ,-0.0003571087848996894, & + -0.0000876148051892879, 5.252431910751829e-6, & + 1.579762259448864e-6 ,-3.466867400295792e-8, & + -1.687643078774232e-8 /) & + ,gammap = (/ -5.526396144304812e-6 , 4.885838128243163e-8, & + 9.96026931578033e-9 ,-7.251389796582352e-10, & + -3.987360250058777e-11, 4.006307891935698e-12, & + 8.26367520608008e-13 /) + +! --- auxiliary statements for finding root of cubic polynomial + a0(s,r)=(c1+s*(c3+s*c8)-r)*rc6 !constant coefficient + a1(s) =(c2+s*(c5+s*c9) )*rc6 !linear coefficient + a2(s) =(c4+s* c7 )*rc6 !quadratic coefficient + !cubic coefficient is c6*rc6=1.0 + cubq(s) =a3rd* a1(s) -(a3rd*a2(s))**2 + cubr(r,s)=a3rd*(0.5*a1(s)*a2(s)-1.5*a0(s,r))-(a3rd*a2(s))**3 +! --- if q**3+r**2>0, water is too dense to yield real root at given +! --- salinitiy. setting q**3+r**2=0 in that case is equivalent to +! --- lowering sigma until a double real root is obtained. + cuban(r,s)=a3rd*atan2(sqrt(max(0.d0,-(cubq(s)**3+cubr(r,s)**2))), & + cubr(r,s)) + cubrl(r,s)=sqrt(-cubq(s))*cos(cuban(r,s)) + cubim(r,s)=sqrt(-cubq(s))*sin(cuban(r,s)) + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma-theta as a function of temp (deg c) and salinity (psu) +! --- (polynomial fit that is cubic in T and quadratic in S) + + sig(t,s)=(c1+s*(c3+s* c8)+ & + t*(c2+s*(c5+s*c9)+t*(c4+s*c7+t*c6))) + +! --- d(sig)/dt + dsigdt(t,s)=(c2+s*(c5+s*c9)+2.0*t*(c4+s*c7+1.5*t*c6)) + +! --- d(sig)/ds + dsigds(t,s)=(c3+t*(c5+t*c7)+2.0*s* c8) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a cubic polynominal root of t**3+a2*t**2+a1*t+a0=0 + tofsig(r,s)=-cubrl(r,s)+sqrt(3.0)*cubim(r,s)-a3rd*a2(s) + +! --- salinity (psu) as a function of sigma and temperature (deg c) +! --- find a quadratic polynominal root of a*s**2+b*s+c=0 + sofsig_a(r,t)=(c8+t* c9) !quadratic coefficient + sofsig_b(r,t)=(c3+t*(c5+t* c7)) !linear coefficient + sofsig_c(r,t)=(c1+t*(c2+t*(c4+t*c6))-r) !constant coefficient + sofsig(r,t)=(2.d0*sofsig_c(r,t))/ & + (-sofsig_b(r,t) & + -sign(sqrt(sofsig_b(r,t)**2- & + 4.d0*sofsig_a(r,t)*sofsig_c(r,t)), & + sofsig_b(r,t))) + +! --- locally referenced sigma, a fit towards Jackett & McDougall (1995) +! --- t: potential temperature; s: psu; prs: pressure + c1l(prs)=alphap(1)+prs2pb*prs*(betap(1)+prs2pb*prs*gammap(1)) + c2l(prs)=alphap(2)+prs2pb*prs*(betap(2)+prs2pb*prs*gammap(2)) + c3l(prs)=alphap(3)+prs2pb*prs*(betap(3)+prs2pb*prs*gammap(3)) + c4l(prs)=alphap(4)+prs2pb*prs*(betap(4)+prs2pb*prs*gammap(4)) + c5l(prs)=alphap(5)+prs2pb*prs*(betap(5)+prs2pb*prs*gammap(5)) + c6l(prs)=alphap(6)+prs2pb*prs*(betap(6)+prs2pb*prs*gammap(6)) + c7l(prs)=alphap(7)+prs2pb*prs*(betap(7)+prs2pb*prs*gammap(7)) + sigloc(t,s,prs)=c1l(prs)+c3l(prs)*s+ & + t*(c2l(prs)+c5l(prs)*s+t*(c4l(prs)+c7l(prs)*s+c6l(prs)*t)) + dsiglocdt(t,s,prs)=(c2l(prs)+c5l(prs)*s+ & + 2.0*t*(c4l(prs)+c7l(prs)*s+1.5*c6l(prs)*t)) + dsiglocds(t,s,prs)=(c3l(prs)+t*(c5l(prs)+t*c7l(prs))) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA4_12term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA4_12term.h90 new file mode 100644 index 00000000..a1333df0 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA4_12term.h90 @@ -0,0 +1,155 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver=48 !12-term sigma-4 +!sig2& sigver= 8 !12-term sigma-2 +!sig0& sigver= 7 !12-term sigma-0 + + real(kind=8) :: sig,dsigdt,dsigds,tofsig,sofsig,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sig_n,sig_d,sig_q, dsigdt_n,dsigdt_d, dsigds_n,dsigds_d + real(kind=8) :: tofsig_a,tofsig_b,tofsig_c + real(kind=8) :: sofsig_a,sofsig_b,sofsig_c + real(kind=8) :: sigloc_n,sigloc_d,sigloc_q,dsiglocdt_n,dsiglocdt_d, dsiglocds_n,dsiglocds_d + real(kind=8) :: r,s,t,pdb,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + aone =1.0, & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- REFERENCE? + +! --- coefficients for 18-term rational function sigloc(). + real(kind=8), parameter :: & + c001=-1.4627567840659594d-01, & !num. constant coefficent + c002= 6.4247392832635697d-02, & !num. T coefficent + c003= 8.1213979591704621d-01, & !num. S coefficent + c004=-8.1321489441909698d-03, & !num. T^2 coefficent + c005= 4.5199845091090296d-03, & !num. T S coefficent + c006= 4.6347888132781394d-04, & !num. S^2 coefficent + c007= 5.0879498675039621d-03, & !num. P coefficent + c008= 1.6333913018305079d-05, & !num. P T coefficent + c009= 4.3899924880543972d-06 !num. P S coefficent + real(kind=8), parameter :: & + c011= 1.0000000000000000d+00, & !den. constant coefficent + c012= 1.0316374535350838d-02, & !den. T coefficent + c013= 8.9521792365142522d-04, & !den. S coefficent + c014=-2.8438341552142710d-05, & !den. T^2 coefficent + c015=-1.1887778959461776d-05, & !den. T S coefficent + c016=-4.0163964812921489d-06, & !den. S^2 coefficent + c017= 1.1995545126831476d-05, & !den. P coefficent + c018= 5.5234008384648383d-08, & !den. P T coefficent + c019= 8.4310335919950873d-09 !den. P S coefficent + real(kind=8), parameter :: & + c004x2=c004*2.d0, & !for dsigdt and dsiglocdt + c014x2=c014*2.d0, & !for dsigdt and dsiglocdt + c006x2=c006*2.d0, & !for dsigds and dsiglocds + c016x2=c016*2.d0 !for dsigds and dsiglocds + real(kind=8), parameter :: & + sqrmin=0.d0, & !sqrt arg can't be negative + sofmin=0.d0 !salinity can't be negative +! --- reference pressure. + real(kind=8), parameter :: prs2pdb=1.d-4 !Pascals to dbar +!sig0 real(kind=8), parameter :: rpdb=0.d0 !reference pressure in dbar, sigma-0 +!sig2 real(kind=8), parameter :: rpdb=2000.d0 !reference pressure in dbar, sigma-2 + real(kind=8), parameter :: rpdb=4000.d0 !reference pressure in dbar, sigma-4 +! --- coefficients for 12-term rational function sig() at rpdb. + real(kind=8), parameter :: & + c101=c001+rpdb*c007, & !num. constant coefficent + c102=c002+rpdb*c008, & !num. T coefficent + c103=c003+rpdb*c009 !num. S coefficent + real(kind=8), parameter :: & + c111=c011+rpdb*c017, & !num. constant coefficent + c112=c012+rpdb*c018, & !num. T coefficent + c113=c013+rpdb*c019 !num. S coefficent + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma at rpdb (dbar) as a function of temp (deg c) and salinity (psu) + + sig_n(t,s) = c101+(c102+c004*t+c005*s)*t + & + (c103+ c006*s)*s + sig_d(t,s) = c111+(c112+c014*t+c015*s)*t + & + (c113 +c016*s)*s + sig_q(t,s) = aone/sig_d(t,s) + sig( t,s) = sig_n(t,s)*sig_q(t,s) + +! --- d(sig)/dt + dsigdt_n(t,s) = c102+c004x2*t+c005*s + dsigdt_d(t,s) = c112+c014x2*t+c015*s + dsigdt( t,s) = (dsigdt_n(t,s)- & + dsigdt_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- d(sig)/ds + dsigds_n(t,s) = c103+c005*t+c006x2*s + dsigds_d(t,s) = c113+c015*t+c016x2*s + dsigds( t,s) = (dsigds_n(t,s)- & + dsigds_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- temp (deg c) as a function of sigma and salinity (psu) +! --- find a quadratic polynominal root of a*t**2+b*t+c=0 + tofsig_a(r,s)=( c004 - & + r* c014 ) !quadratic coefficient + tofsig_b(r,s)=( (c102+ c005*s) - & + r*(c112+ c015*s) ) !linear coefficient + tofsig_c(r,s)=( (c101+(c103+c006*s)*s) - & + r*(c111+(c113+c016*s)*s) ) !constant coefficient + tofsig(r,s)=( -tofsig_b(r,s) & + -sqrt(max(sqrmin, & + tofsig_b(r,s)**2 - & + 4.0*tofsig_a(r,s)*tofsig_c(r,s))) ) / & + (2.0*tofsig_a(r,s)) + +! --- salinity (psu) as a function of sigma and temperature (deg c) +! --- find a quadratic polynominal root of a*s**2+b*s+c=0 + sofsig_a(r,t)=( c006 - & + r* c016 ) !quadratic coefficient + sofsig_b(r,t)=( (c103+ c005*t) - & + r*(c113+ c015*t) ) !linear coefficient + sofsig_c(r,t)=( (c101+(c102+c004*t)*t) - & + r*(c111+(c112+c014*t)*t) ) !constant coefficient + sofsig(r,s)=max(sofmin, & + ( -sofsig_b(r,s) & + +sqrt(max(sqrmin, & + sofsig_b(r,s)**2 - & + 4.0*sofsig_a(r,s)*sofsig_c(r,s))) ) / & + (2.0*sofsig_a(r,s)) ) + +! --- locally referenced sigma, using the 18-term equation of state. +! --- t: potential temperature; s: psu; prs: pressure + + sigloc_n(t,s,pdb) = c001+(c002+c004*t+c005*s)*t + & + (c003+ c006*s)*s + & + (c007+c008*t+c009*s)*pdb + sigloc_d(t,s,pdb) = c011+(c012+c014*t+c015*s)*t + & + (c013 +c016*s)*s + & + (c017+c018*t+c019*s)*pdb + sigloc_q(t,s,pdb) = aone/sigloc_d(t,s,pdb) + sigloc( t,s,prs)=sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/dt + dsiglocdt_n(t,s,pdb) = c002+c004x2*t+c005*s+c008*pdb + dsiglocdt_d(t,s,pdb) = c012+c014x2*t+c015*s+c018*pdb + dsiglocdt( t,s,prs)=(dsiglocdt_n(t,s,prs*prs2pdb)- & + dsiglocdt_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/ds + dsiglocds_n(t,s,pdb) = c003+c005*t+c006x2*s+c009*pdb + dsiglocds_d(t,s,pdb) = c013+c015*t+c016x2*s+c019*pdb + dsiglocds( t,s,prs)=(dsiglocds_n(t,s,prs*prs2pdb)- & + dsiglocds_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA4_17term.h90 b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA4_17term.h90 new file mode 100644 index 00000000..7930370c --- /dev/null +++ b/hycom/MSCPROGS/src/Average/include/stmt_fns_SIGMA4_17term.h90 @@ -0,0 +1,208 @@ +!----------------------------------------------------------------------------- + integer, parameter :: & + sigver=46 !17-term sigma-4 +!sig2& sigver= 6 !17-term sigma-2 +!sig0& sigver= 5 !17-term sigma-0 + + real(kind=8) :: sig,dsigdt,dsigds,sigloc,dsiglocdt,dsiglocds + real(kind=8) :: sig_n,sig_d,sig_q, dsigdt_n,dsigdt_d, dsigds_n,dsigds_d + real(kind=8) :: sigloc_n,sigloc_d,sigloc_q,dsiglocdt_n,dsiglocdt_d, dsiglocds_n,dsiglocds_d + real(kind=8) :: r,s,t,pdb,prs + real(kind=8) :: r8 + real(kind=4) :: r4 + + real(kind=8), parameter :: & + aone =1.0, & + ahalf=1.0/2.0, & + a3rd =1.0/3.0, athird =a3rd, & + a4th =1.0/4.0, afourth=a4th + +! --- Jackett, McDougall, Feistel, Wright and Griffies (2006), +! --- Algorithms for Density, Potential Temperature, Conservative +! --- Temperature, and the Freezing Temperature of Seawater, JAOT + +! --- coefficients for 25-term rational function sigloc(). + real(kind=8), parameter :: & + c001= 9.9984085444849347d+02, & !num. constant coefficent + c002= 7.3471625860981584d+00, & !num. T coefficent + c003=-5.3211231792841769d-02, & !num. T^2 coefficent + c004= 3.6492439109814549d-04, & !num. T^3 coefficent + c005= 2.5880571023991390d+00, & !num. S coefficent + c006= 6.7168282786692355d-03, & !num. T S coefficent + c007= 1.9203202055760151d-03, & !num. S^2 coefficent + c008= 1.0000000000000000d+00, & !den. constant coefficent + c009= 7.2815210113327091d-03, & !den. T coefficent + c010=-4.4787265461983921d-05, & !den. T^2 coefficent + c011= 3.3851002965802430d-07, & !den. T^3 coefficent + c012= 1.3651202389758572d-10, & !den. T^4 coefficent + c013= 1.7632126669040377d-03, & !den. S coefficent + c014= 8.8066583251206474d-06, & !den. T S coefficent + c015= 1.8832689434804897d-10, & !den. T^3S coefficent + c016= 5.7463776745432097d-06, & !den. T S^1.5 coefficent + c017= 1.4716275472242334d-09 !den. T^3S^1.5 coefficent + real(kind=8), parameter :: & + c018= 1.1798263740430364d-02, & !num. P coefficent + c019= 9.8920219266399117d-08, & !num. P T^2 coefficent + c020= 4.6996642771754730d-06, & !num. P S coefficent + c021= 2.5862187075154352d-08, & !num. P^2 coefficent + c022= 3.2921414007960662d-12, & !num. P^2T^2 coefficent + c023= 6.7103246285651894d-06, & !den. P coefficent + c024= 2.4461698007024582d-17, & !den. P^2T^3 coefficent + c025= 9.1534417604289062d-18 !den. P^3T coefficent +! --- additional coefficients for dsiglocdt(). + real(kind=8), parameter :: & + c031= 7.3471625860981580d+00, & !num. constant coefficent + c032=-1.0642246358568354d-01, & !num. T coefficent + c033= 1.0947731732944364d-03, & !num. T^2 coefficent + c034= 6.7168282786692355d-03, & !num. S coefficent + c035= 7.2815210113327090d-03, & !den. constant coefficent + c036=-8.9574530923967840d-05, & !den. T coefficent + c037= 1.0155300889740728d-06, & !den. T^2 coefficent + c038= 5.4604809559034290d-10, & !den. T^3 coefficent + c039=-8.8066583251206470d-06, & !den. S coefficent + c040= 5.6498068304414700d-10, & !den. T^2S coefficent + c041= 2.9432550944484670d-09, & !den. T S^1.5 coefficent + c042= 1.9784043853279823d-07, & !num. P T coefficent + c043= 6.5842828015921320d-12, & !num. P^2T coefficent + c044= 7.3385094021073750d-17, & !den. P^2T^2 coefficent + c045= 9.1534417604289060d-18 !den. P^3 coefficent +! --- additional coefficients for dsiglocds(). + real(kind=8), parameter :: & + c051= 2.5880571023991390d+00, & !num. constant coefficent + c052= 6.7168282786692355d-03, & !num. T coefficent + c053= 3.8406404111520300d-03, & !num. S coefficent + c054= 1.7632126669040377d-03, & !den. constant coefficent + c055=-8.8066583251206470d-06, & !den. T coefficent + c056= 1.8832689434804897d-10, & !den. T^3 coefficent + c057= 8.6195665118148150d-06, & !den. S^0.5 coefficent + c058= 2.2074413208363504d-09, & !den. T^2S^0.5 coefficent + c059= 4.6996642771754730d-06 !num. P coefficent + + real(kind=8), parameter :: sqrmin=0.d0 !sqrt arg can't be negative +! --- reference pressure. + real(kind=8), parameter :: prs2pdb=1.d-4 !Pascals to dbar +!sig0 real(kind=8), parameter :: rpdb=0.d0 !reference pressure in dbar, sigma0 +!sig2 real(kind=8), parameter :: rpdb=2000.d0 !reference pressure in dbar, sigma2 + real(kind=8), parameter :: rpdb=4000.d0 !reference pressure in dbar, sigma4 +! --- coefficients for 17-term rational function sig() at rpdb. + real(kind=8), parameter :: & + c101=c001+(c018-c021*rpdb)*rpdb, & !num. constant coefficent + c103=c003+(c019-c022*rpdb)*rpdb, & !num. T^2 coefficent + c105=c005+c020*rpdb, & !num. S coefficent + c108=c008+c023*rpdb, & !den. constant coefficent + c109=c009-c025*rpdb**3, & !den. T coefficent + c111=c011-c024*rpdb**2 !den. T^3 coefficent +! --- additional coefficients for dsigdt(). + real(kind=8), parameter :: & + c132=c032+(c042-c043*rpdb)*rpdb, & !num. T coefficent + c135=c035-c045*rpdb**3, & !den. constant coefficent + c137=c037-c044*rpdb**2 !den. T^2 coefficent +! --- additional coefficients for dsigds(). + real(kind=8), parameter :: & + c151=c051+c059*rpdb !num. constant coefficent + real, parameter :: & + rhoref=1.d3 !rhoref=qthref kg/m^3 + +! --- ----------------- +! --- equation of state +! --- ----------------- + +! --- sigma at rpdb (dbar) as a function of pot.temp (deg c) and salinity (psu) + + sig_n(t,s) = c101 + t*(c002+t*(c103+t*c004)) + & + s*(c105-t*c006+s*c007) + sig_d(t,s) = c108 + t*(c109+t*(c010+t*(c111+t*c012))) + & + s*(c013-t*(c014+t*t*c015) + & + sqrt(max(sqrmin,s))*(c016+t*t*c017)) + sig_q(t,s) = aone/sig_d(t,s) + sig( t,s) = sig_n(t,s)*sig_q(t,s) - rhoref + +! --- d(sig)/dt + dsigdt_n(t,s) = c031 + t*(c132+t*c033) - & + s* c034 + dsigdt_d(t,s) = c135 + t*(c036+t*(c137+t*c038)) + & + s*(c039-t*t*c040+ & + sqrt(max(sqrmin,s))*t*c041 ) + dsigdt( t,s) = (dsigdt_n(t,s)- & + dsigdt_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- d(sig)/ds +! --- additional coefficients for dsigds(). + dsigds_n(t,s) = c151 - t*c052 + s*c053 + dsigds_d(t,s) = c054 + t*(c055-t*t*c056) + & + sqrt(max(sqrmin,s))*(c057+t*t*c058) + dsigds( t,s) = (dsigds_n(t,s)- & + dsigds_d(t,s)*sig_n(t,s)*sig_q(t,s))*sig_q(t,s) + +! --- locally referenced sigma, using the 25-term equation of state. +! --- t: potential temperature (degC); s: salinity (psu); prs: pressure (dbar) + sigloc_n(t,s,pdb) = c001 + & + t*(c002 + & + t*(c003 + & + t* c004 )) + & + s*(c005 - & + t* c006 + & + s* c007 ) + & + pdb*(c018 + & + t*t* c019 + & + s* c020 - & + pdb*(c021 + & + t*t* c022 )) + sigloc_d(t,s,pdb) = c008 + & + t*(c009 + & + t*(c010 + & + t*(c011 + & + t* c012 ))) + & + s*(c013 - & + t*(c014 + & + t*t* c015 ) + & + sqrt(max(sqrmin,s))*(c016 + & + t*t* c017 )) + & + pdb*(c023 - & + pdb*t*(t*t* c024 + & + pdb* c025 )) + sigloc_q(t,s,pdb) = aone/sigloc_d(t,s,pdb) + sigloc(t,s,prs)=sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) - rhoref + +! --- d(sig)/dt + dsiglocdt_n(t,s,pdb) = c031 + & + t*(c032 + & + t* c033 ) - & + s* c034 + & + pdb*t*(c042 - & + pdb* c043 ) + dsiglocdt_d(t,s,pdb) = c035 + & + t*(c036 + & + t*(c037 + & + t* c038 )) + & + s*(c039 - & + t*t* c040 + & + sqrt(max(sqrmin,s))*t* c041 ) - & + pdb*pdb*(t*t* c044 + & + pdb* c045 ) + dsiglocdt(t,s,prs)=(dsiglocdt_n(t,s,prs*prs2pdb)- & + dsiglocdt_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- d(sig)/ds + dsiglocds_n(t,s,pdb) = c051 - & + t* c052 + & + s* c053 + & + pdb* c059 + dsiglocds_d(t,s,pdb) = c054 + & + t*(c055 - & + t*t* c056 ) + & + sqrt(max(sqrmin,s))*(c057 + & + t*t* c058 ) + dsiglocds(t,s,prs)=(dsiglocds_n(t,s,prs*prs2pdb)- & + dsiglocds_d(t,s,prs*prs2pdb)* & + sigloc_n(t,s,prs*prs2pdb)* & + sigloc_q(t,s,prs*prs2pdb) ) * & + sigloc_q(t,s,prs*prs2pdb) + +! --- auxiliary statement for real to real(kind=8) conversion + r8(r4) = r4 +!----------------------------------------------------------------------------- diff --git a/hycom/MSCPROGS/src/Average/makefile b/hycom/MSCPROGS/src/Average/makefile index d75b3ce4..01ee6ea6 100755 --- a/hycom/MSCPROGS/src/Average/makefile +++ b/hycom/MSCPROGS/src/Average/makefile @@ -33,7 +33,8 @@ endif TARGET = hycave TARGET2 = ensave -all: $(TARGET) $(TARGET2) +TARGET3 = hycom_mean +all: $(TARGET) $(TARGET2) $(TARGET3) ############################################################################### @@ -48,10 +49,16 @@ OBJECTS2= p_ensave.o $(TARGET2): $(OBJECTS2) cd ./TMP ; $(LD) $(LINKFLAGS) -o ../$(TARGET2) $(OBJECTS2) $(LIBS) ############################################################################### +############################################################################### +OBJECTS3= mod_mean.o bigrid.o blkin.o extrct.o getdat.o putdat.o hycom_mean.o + +$(TARGET3): $(OBJECTS3) + cd ./TMP ; $(LD) $(LINKFLAGS) -o ../$(TARGET3) $(OBJECTS3) $(LIBS) + ############################################################################### clean: rm ./TMP/*.o TMP/*.mod TMP/*.f TMP/*.f90 $(TARGET) $(TARGET2) install : all mkdir -p ../../bin - cp $(TARGET) $(TARGET2) ../../bin + cp $(TARGET) $(TARGET2) $(TARGET3) ../../bin diff --git a/hycom/MSCPROGS/src/Average/mean_hycom.in b/hycom/MSCPROGS/src/Average/mean_hycom.in new file mode 100644 index 00000000..bd3235b3 --- /dev/null +++ b/hycom/MSCPROGS/src/Average/mean_hycom.in @@ -0,0 +1,8 @@ + 50 'kk ' = number of layers involved + 0 'meansq' = form meansq, rather than mean (0=F,1=T) + 2 'narchs' = number of archives to read (==0 to end input) +file_mem001.a +file_mem002.a + 0 'narchs' = number of archives to read (==0 to end input) +file + diff --git a/hycom/MSCPROGS/src/Average/mod_mean.F90 b/hycom/MSCPROGS/src/Average/mod_mean.F90 new file mode 100644 index 00000000..3b5bafac --- /dev/null +++ b/hycom/MSCPROGS/src/Average/mod_mean.F90 @@ -0,0 +1,1254 @@ +module mod_mean + implicit none + + ! --- HYCOM mean: array allocation and calculation interface. + + ! --- ii = 1st dimension of array (==idm) + ! --- jj = 2nd dimension of array (==jdm) + ! --- kk = number of layers (typically 1) + ! --- nmean = number of archive records in the mean + ! --- ntracr= number of tracers to form the mean of + + integer, save :: ii,ii1,ii2,iorign,jj,jj1,jj2,jorign,kk,ntracr + integer, save :: nmean,nstep,nstep_m + + ! --- lsteric = steric in input (and output) means + ! --- lwtrflx = wtrflx in input (and output) means + + logical, save :: loneta, lsteric, lwtrflx + logical, save :: meansq, single, trcout, icegln, hisurf, bgcout + + ! --- archive header + + character(len=80), save, dimension(4) :: ctitle + + ! --- arrays: + + real, save, allocatable, dimension (:,:,:,:) :: & + tracer, & + tracer_m + + real, save, allocatable, dimension (:,:,:) :: & + u, v, ke, temp, saln, th3d, dp ,visc, tdff, sdff, dw,p, & + u_m,v_m,ke_m,temp_m,saln_m,th3d_m,dp_m,visc_m,tdff_m,sdff_m,dpu_m,dpv_m,dw_m + + real, save, allocatable, dimension (:,:) :: & + depths,depthu,depthv, & + ubaro,vbaro,pbaro,kebaro, & + montg,srfht,steric,oneta,onetaw,dpbl,dpmixl, & + tmix,smix,thmix,umix,vmix,kemix, & + surflx,salflx,wtrflx,covice,thkice,temice, & + ubaro_m,vbaro_m,pbaro_m,kebaro_m, & + montg_m,srfht_m,steric_m,dpbl_m,dpmixl_m, & + tmix_m,smix_m,thmix_m,umix_m,vmix_m,kemix_m, & + surflx_m,salflx_m,wtrflx_m,covice_m,thkice_m,temice_m, & + oneta_m,onetaw_m,onetaw_u,onetaw_v, & + si_u, si_um, si_v, si_vm, & + surtx, surtxm, surty, surtym + ! --- add two sea-ice variables for drift + ! --- add x-/y surface stress + + real, save, allocatable, dimension (:) :: & + theta + + integer, save, allocatable, dimension (:,:) :: & + ip,iq,iu,iv + + ! --- ECOSMO + + integer, parameter :: ntracr_bgc_2d = 10, ntracr_bgc_3d = 45 + + character(len=8), save, allocatable, dimension(:) :: & + nvar_bgc_2d + + character(len=8), save, allocatable, dimension(:) :: & + nvar_bgc_3d + + real, save, allocatable, dimension (:,:,:,:) :: & + tracer_bgc_3d, & + tracer_bgc_3dm + + real, save, allocatable, dimension (:,:,:) :: & + tracer_bgc_2d, & + tracer_bgc_2dm + + ! --- module subroutines + + contains + + subroutine mean_alloc + implicit none + real, parameter :: spval = 2.0**100 + + ! --- initialize allocatable arrays. + + ii1 = ii - 1 + ii2 = ii - 2 + jj1 = jj - 1 + jj2 = jj - 2 + + nmean = 0 + + lsteric = .false. !default + lwtrflx = .false. !default + loneta = .false. !default + + ! --- + + allocate( p( ii, jj,kk+1) ) + allocate( dpu_m( ii, jj,kk ) ); dpu_m = 0.0 + allocate( dpv_m( ii, jj,kk ) ); dpv_m = 0.0 + allocate( depthu( ii, jj ) ) + allocate( depthv( ii, jj ) ) + allocate( depths(0:ii,0:jj ) ) + + allocate( onetaw_u(ii,jj) ); onetaw_u = 0.0 + allocate( onetaw_v(ii,jj) ); onetaw_v = 0.0 + + allocate( ip(ii,jj) ) + allocate( iq(ii,jj) ) + allocate( iu(ii,jj) ) + allocate( iv(ii,jj) ) + + allocate( theta(kk) ) + + ! --- 2D + + allocate( montg(ii,jj), montg_m(ii,jj) ); montg_m = 0.0 + allocate( srfht(ii,jj), srfht_m(ii,jj) ); srfht_m = 0.0 + allocate( steric(ii,jj), steric_m(ii,jj) ); steric_m = 0.0 + allocate( oneta(ii,jj), oneta_m(ii,jj) ); oneta_m = 0.0 + allocate( onetaw(ii,jj), onetaw_m(ii,jj) ); onetaw_m = 0.0 + allocate( surflx(ii,jj), surflx_m(ii,jj) ); surflx_m = 0.0 + allocate( wtrflx(ii,jj), wtrflx_m(ii,jj) ); wtrflx_m = 0.0 + allocate( salflx(ii,jj), salflx_m(ii,jj) ); salflx_m = 0.0 + allocate( dpbl(ii,jj), dpbl_m(ii,jj) ); dpbl_m = 0.0 + allocate( dpmixl(ii,jj), dpmixl_m(ii,jj) ); dpmixl_m = 0.0 + allocate( tmix(ii,jj), tmix_m(ii,jj) ); tmix_m = 0.0 + allocate( smix(ii,jj), smix_m(ii,jj) ); smix_m = 0.0 + allocate( thmix(ii,jj), thmix_m(ii,jj) ); thmix_m = 0.0 + allocate( umix(ii,jj), umix_m(ii,jj) ); umix_m = 0.0 + allocate( vmix(ii,jj), vmix_m(ii,jj) ); vmix_m = 0.0 + allocate( kemix(ii,jj), kemix_m(ii,jj) ); kemix_m = 0.0 + + allocate( covice(ii,jj), covice_m(ii,jj) ); covice_m = 0.0 + allocate( thkice(ii,jj), thkice_m(ii,jj) ); thkice_m = 0.0 + allocate( temice(ii,jj), temice_m(ii,jj) ); temice_m = 0.0 + allocate( si_u(ii,jj), si_um(ii,jj) ); si_um = 0.0 + allocate( si_v(ii,jj), si_vm(ii,jj) ); si_vm = 0.0 + allocate( surtx(ii,jj), surtxm(ii,jj) ); surtxm = 0.0 + allocate( surty(ii,jj), surtym(ii,jj) ); surtym = 0.0 + + allocate( ubaro(ii,jj), ubaro_m(ii,jj) ); ubaro_m = 0.0 + allocate( vbaro(ii,jj), vbaro_m(ii,jj) ); vbaro_m = 0.0 + allocate( pbaro(ii,jj), pbaro_m(ii,jj) ); pbaro_m = 0.0 + allocate( kebaro(ii,jj), kebaro_m(ii,jj) ); kebaro_m = 0.0 + + ! 3D + + allocate( u(ii,jj,kk), u_m(ii,jj,kk) ); u_m = 0.0 + allocate( v(ii,jj,kk), v_m(ii,jj,kk) ); v_m = 0.0 + allocate( ke(ii,jj,kk), ke_m(ii,jj,kk) ); ke_m = 0.0 + allocate( dp(ii,jj,kk), dp_m(ii,jj,kk) ); dp_m = 0.0 + allocate( dw(ii,jj,kk), dw_m(ii,jj,kk) ); dw_m = 0.0 + allocate( temp(ii,jj,kk), temp_m(ii,jj,kk) ); temp_m = 0.0 + allocate( saln(ii,jj,kk), saln_m(ii,jj,kk) ); saln_m = 0.0 + allocate( th3d(ii,jj,kk), th3d_m(ii,jj,kk) ); th3d_m = 0.0 + allocate( visc(ii,jj,kk), visc_m(ii,jj,kk) ); visc_m = 0.0 + allocate( tdff(ii,jj,kk), tdff_m(ii,jj,kk) ); tdff_m = 0.0 + allocate( sdff(ii,jj,kk), sdff_m(ii,jj,kk) ); sdff_m = 0.0 + + if (trcout) then + allocate( tracer( ii,jj,kk,ntracr) ) + allocate( tracer_m(ii,jj,kk,ntracr) ); tracer_m = 0.0 + endif + + ! --- ECOSMO + + if (bgcout) then + allocate( nvar_bgc_2d(ntracr_bgc_2d) ) + allocate( nvar_bgc_3d(ntracr_bgc_3d) ) + + allocate( tracer_bgc_2d(ii,jj, ntracr_bgc_2d)) + allocate( tracer_bgc_2dm(ii,jj, ntracr_bgc_2d)); tracer_bgc_2dm(:,:,:) = 0.0 + allocate( tracer_bgc_3d(ii,jj,kk,ntracr_bgc_3d)) + allocate( tracer_bgc_3dm(ii,jj,kk,ntracr_bgc_3d)); tracer_bgc_3dm(:,:,:,:) = 0.0 + + ! register variable names (TODO: register in namelist) + + nvar_bgc_3d( 1) = 'CO2_c ' + nvar_bgc_3d( 2) = 'CO2_TA ' + nvar_bgc_3d( 3) = 'ECO_no3 ' + nvar_bgc_3d( 4) = 'ECO_nh4 ' + nvar_bgc_3d( 5) = 'ECO_pho ' + nvar_bgc_3d( 6) = 'ECO_sil ' + nvar_bgc_3d( 7) = 'ECO_oxy ' + nvar_bgc_3d( 8) = 'ECO_fla ' + nvar_bgc_3d( 9) = 'ECO_dia ' + nvar_bgc_3d(10) = 'ECO_ccl ' + nvar_bgc_3d(11) = 'ECO_cclc' + nvar_bgc_3d(12) = 'ECO_caco' + nvar_bgc_3d(13) = 'ECO_diac' + nvar_bgc_3d(14) = 'ECO_flac' + nvar_bgc_3d(15) = 'ECO_micr' + nvar_bgc_3d(16) = 'ECO_meso' + nvar_bgc_3d(17) = 'ECO_det ' + nvar_bgc_3d(18) = 'ECO_opa ' + nvar_bgc_3d(19) = 'ECO_dom ' + nvar_bgc_3d(20) = 'ECO_dsnk' + nvar_bgc_3d(21) = 'CO2_pH ' + nvar_bgc_3d(22) = 'CO2_pCO2' + nvar_bgc_3d(23) = 'CO2_Carb' + nvar_bgc_3d(24) = 'CO2_BiCa' + nvar_bgc_3d(25) = 'CO2_Carb' + nvar_bgc_3d(26) = 'CO2_Om_c' + nvar_bgc_3d(27) = 'CO2_Om_a' + nvar_bgc_3d(28) = 'ECO_prim' + nvar_bgc_3d(29) = 'ECO_secp' + nvar_bgc_3d(30) = 'ECO_netp' + nvar_bgc_3d(31) = 'ECO_parm' + nvar_bgc_3d(32) = 'ECO_Nlim' + nvar_bgc_3d(33) = 'ECO_Plim' + nvar_bgc_3d(34) = 'ECO_Slim' + nvar_bgc_3d(35) = 'ECO_Llim' + nvar_bgc_3d(36) = 'ECO_deni' + nvar_bgc_3d(37) = 'ECO_snks' + nvar_bgc_3d(38) = 'ECO_c2ch' + nvar_bgc_3d(39) = 'ECO_c2ch' + nvar_bgc_3d(40) = 'ECO_c2ch' + nvar_bgc_3d(41) = 'light_sw' + nvar_bgc_3d(42) = 'light_pa' + nvar_bgc_3d(43) = 'attenuat' + nvar_bgc_3d(44) = 'total_ch' + nvar_bgc_3d(45) = 'total_ca' + + nvar_bgc_2d( 1) = 'ECO_sed4' + nvar_bgc_2d( 2) = 'ECO_sed1' + nvar_bgc_2d( 3) = 'ECO_sed2' + nvar_bgc_2d( 4) = 'ECO_sed3' + nvar_bgc_2d( 5) = 'CO2_fair' + nvar_bgc_2d( 6) = 'CO2_wind' + nvar_bgc_2d( 7) = 'ECO_bots' + nvar_bgc_2d( 8) = 'light_pa' + nvar_bgc_2d( 9) = 'surface_' + nvar_bgc_2d(10) = 'surface_' + endif + + end subroutine mean_alloc + + subroutine mean_add(iweight) + implicit none + + integer, intent(in) :: iweight + + ! --- add an archive to the mean. + ! --- layer quantities weighted by layer thickness (i.e. by dw). + + integer :: i,im,j,jm,k,ktr + real :: s,swk + real, dimension(kk) :: sw + + nmean = nmean + iweight + + s = iweight + + p(:,:,:) = 0.0 + do j= 1,jj + do i= 1,ii + if (ip(i,j).eq.1) then + p(i,j,1) = 0.0 + do k= 1,kk + p(i,j,k+1) = p(i,j,k) + dw(i,j,k) + enddo + ! else + ! p(i,j,:) = 0.0 + endif + enddo + enddo + + do j= 1,jj + do i= 1,ii + if (iu(i,j).eq.1) then + ubaro_m(i,j) = ubaro_m(i,j) + ubaro(i,j) * s + umix_m(i,j) = umix_m(i,j) + umix(i,j) * s + + if (i.ne.1) then + im = i-1 + else + im = ii + endif + ! --- depthu is either depths(i,j) or depths(i-1,j) + if (depths(i,j).eq.depths(im,j)) then + onetaw_u(i,j) = 0.5*(onetaw(i,j)+onetaw(im,j)) + elseif (depths(i,j).eq.depthu(i,j)) then + onetaw_u(i,j) = onetaw(i,j) + else + onetaw_u(i,j) = onetaw(im,j) + endif + do k= 1,kk + swk = s*max(0.0, min(depthu(i,j),0.5*(p(i,j,k+1)+p(im,j,k+1))) & + - min(depthu(i,j),0.5*(p(i,j,k )+p(im,j,k ))) )*onetaw_u(i,j) + dpu_m(i,j,k) = dpu_m(i,j,k) + swk + u_m(i,j,k) = u_m(i,j,k) + u(i,j,k) * swk + enddo + endif !iu + + if (iv(i,j).eq.1) then + vbaro_m(i,j) = vbaro_m(i,j) + vbaro(i,j) * s + vmix_m(i,j) = vmix_m(i,j) + vmix(i,j) * s + + if (j.ne.1) then + jm = j-1 + else + jm = jj + endif + ! --- depthv is either depths(i,j) or depths(i,j-1) + if (depths(i,j).eq.depths(i,jm)) then + onetaw_v(i,j) = 0.5*(onetaw(i,j)+onetaw(i,jm)) + elseif (depths(i,j).eq.depthv(i,j)) then + onetaw_v(i,j) = onetaw(i,j) + else + onetaw_v(i,j) = onetaw(i,jm) + endif + do k= 1,kk + swk = s*max(0.0, min(depthv(i,j),0.5*(p(i,j,k+1)+p(i,jm,k+1))) & + - min(depthv(i,j),0.5*(p(i,j,k )+p(i,jm,k ))) )*onetaw_v(i,j) + dpv_m(i,j,k) = dpv_m(i,j,k) + swk + v_m(i,j,k) = v_m(i,j,k) + v(i,j,k) * swk + enddo + endif !iv + + if (ip(i,j).eq.1) then + pbaro_m(i,j) = pbaro_m(i,j) + pbaro(i,j) * s + kebaro_m(i,j) = kebaro_m(i,j) + kebaro(i,j) * s + montg_m(i,j) = montg_m(i,j) + montg(i,j) * s + srfht_m(i,j) = srfht_m(i,j) + srfht(i,j) * s + steric_m(i,j) = steric_m(i,j) + steric(i,j) * s + dpbl_m(i,j) = dpbl_m(i,j) + dpbl(i,j) * s + dpmixl_m(i,j) = dpmixl_m(i,j) + dpmixl(i,j) * s + tmix_m(i,j) = tmix_m(i,j) + tmix(i,j) * s + smix_m(i,j) = smix_m(i,j) + smix(i,j) * s + thmix_m(i,j) = thmix_m(i,j) + thmix(i,j) * s + kemix_m(i,j) = kemix_m(i,j) + kemix(i,j) * s + surflx_m(i,j) = surflx_m(i,j) + surflx(i,j) * s + salflx_m(i,j) = salflx_m(i,j) + salflx(i,j) * s + wtrflx_m(i,j) = wtrflx_m(i,j) + wtrflx(i,j) * s + if (icegln) then + covice_m(i,j) = covice_m(i,j) + covice(i,j) * s + thkice_m(i,j) = thkice_m(i,j) + thkice(i,j) * s + temice_m(i,j) = temice_m(i,j) + temice(i,j) * s + si_um(i,j) = si_um(i,j) + si_u(i,j) * s + si_vm(i,j) = si_vm(i,j) + si_v(i,j) * s + surtxm(i,j) = surtxm(i,j) + surtx(i,j) * s + surtym(i,j) = surtym(i,j) + surty(i,j) * s + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + tracer_bgc_2dm(i,j,ktr) = tracer_bgc_2dm(i,j,ktr) + tracer_bgc_2d(i,j,ktr) * s + enddo !ktr + endif + + oneta_m(i,j) = oneta_m(i,j) + oneta(i,j) * s + onetaw_m(i,j) = onetaw_m(i,j) + onetaw(i,j) * s + + sw(:) = onetaw(i,j) * dw(i,j,:) * s + + dw_m(i,j,:) = dw_m(i,j,:) + sw(:) + dp_m(i,j,:) = dp_m(i,j,:) + sw(:) + temp_m(i,j,:) = temp_m(i,j,:) + temp(i,j,:) * sw(:) + saln_m(i,j,:) = saln_m(i,j,:) + saln(i,j,:) * sw(:) + th3d_m(i,j,:) = th3d_m(i,j,:) + th3d(i,j,:) * sw(:) + ke_m(i,j,:) = ke_m(i,j,:) + ke(i,j,:) * sw(:) + visc_m(i,j,:) = visc_m(i,j,:) + visc(i,j,:) * sw(:) + tdff_m(i,j,:) = tdff_m(i,j,:) + tdff(i,j,:) * sw(:) + sdff_m(i,j,:) = sdff_m(i,j,:) + sdff(i,j,:) * sw(:) + + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,:,ktr) = tracer_m(i,j,:,ktr) + tracer(i,j,:,ktr) * sw(:) + enddo !ktr + endif + + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,:,ktr) = tracer_bgc_3dm(i,j,:,ktr) + tracer_bgc_3d(i,j,:,ktr) * sw(:) + enddo !ktr + endif + endif !ip + enddo !i + enddo !j + + ! write(6,*) 'mean_add - dp_m = ', dp_m(54, 1,1), + ! & dp( 54, 1,1) + + end subroutine mean_add + + subroutine mean_addsq(iweight) + implicit none + + integer, intent(in) :: iweight + + ! --- add an archive squared to the mean. + ! --- layer quantities weighted by layer thickness (i.e. by dw). + + integer :: i,im,j,jm,k,ktr + real :: s,swk + real, dimension(kk) :: sw + + nmean = nmean + iweight + + s = iweight + + do j= 1,jj + do i= 1,ii + if (ip(i,j).eq.1) then + p(i,j,1) = 0.0 + do k= 1,kk + p(i,j,k+1) = p(i,j,k) + dw(i,j,k) + enddo + ! else + ! p(i,j,:) = 0.0 + endif + enddo + enddo + + do j= 1,jj + do i= 1,ii + if (iu(i,j).eq.1) then + ubaro_m(i,j) = ubaro_m(i,j) + ubaro(i,j)**2 * s + umix_m(i,j) = umix_m(i,j) + umix(i,j)**2 * s + + if (i.ne.1) then + im = i-1 + else + im = ii + endif + ! --- depthu is either depths(i,j) or depths(i-1,j) + if (depths(i,j).eq.depths(im,j)) then + onetaw_u(i,j) = 0.5*(onetaw(i,j)+onetaw(im,j)) + elseif (depths(i,j).eq.depthu(i,j)) then + onetaw_u(i,j) = onetaw(i,j) + else + onetaw_u(i,j) = onetaw(im,j) + endif + do k= 1,kk + swk = s*max(0.0, min(depthu(i,j), 0.5*(p(i,j,k+1)+p(im,j,k+1))) & + - min(depthu(i,j), 0.5*(p(i,j,k )+p(im,j,k ))) )*onetaw_u(i,j) + dpu_m(i,j,k) = dpu_m(i,j,k) + swk + u_m(i,j,k) = u_m(i,j,k) + u(i,j,k)**2 * swk + enddo + endif !iu + + if (iv(i,j).eq.1) then + vbaro_m(i,j) = vbaro_m(i,j) + vbaro(i,j)**2 * s + vmix_m(i,j) = vmix_m(i,j) + vmix(i,j)**2 * s + + if (j.ne.1) then + jm = j-1 + else + jm = jj + endif + ! --- depthv is either depths(i,j) or depths(i,j-1) + if (depths(i,j).eq.depths(i,jm)) then + onetaw_v(i,j) = 0.5*(onetaw(i,j)+onetaw(i,jm)) + elseif (depths(i,j).eq.depthv(i,j)) then + onetaw_v(i,j) = onetaw(i,j) + else + onetaw_v(i,j) = onetaw(i,jm) + endif + do k= 1,kk + swk = s*max( 0.0, min(depthv(i,j),0.5*(p(i,j,k+1)+p(i,jm,k+1))) & + - min(depthv(i,j),0.5*(p(i,j,k )+p(i,jm,k ))) )*onetaw_v(i,j) + dpv_m(i,j,k) = dpv_m(i,j,k) + swk + v_m(i,j,k) = v_m(i,j,k) + v(i,j,k)**2 * swk + enddo + endif !iv + + if (ip(i,j).eq.1) then + pbaro_m(i,j) = pbaro_m(i,j) + pbaro(i,j)**2 * s + kebaro_m(i,j) = kebaro_m(i,j) + kebaro(i,j)**2 * s + montg_m(i,j) = montg_m(i,j) + montg(i,j)**2 * s + srfht_m(i,j) = srfht_m(i,j) + srfht(i,j)**2 * s + steric_m(i,j) = steric_m(i,j) + steric(i,j)**2 * s + dpbl_m(i,j) = dpbl_m(i,j) + dpbl(i,j)**2 * s + dpmixl_m(i,j) = dpmixl_m(i,j) + dpmixl(i,j)**2 * s + tmix_m(i,j) = tmix_m(i,j) + tmix(i,j)**2 * s + smix_m(i,j) = smix_m(i,j) + smix(i,j)**2 * s + thmix_m(i,j) = thmix_m(i,j) + thmix(i,j)**2 * s + kemix_m(i,j) = kemix_m(i,j) + kemix(i,j)**2 * s + surflx_m(i,j) = surflx_m(i,j) + surflx(i,j)**2 * s + salflx_m(i,j) = salflx_m(i,j) + salflx(i,j)**2 * s + wtrflx_m(i,j) = wtrflx_m(i,j) + wtrflx(i,j)**2 * s + if (icegln) then + covice_m(i,j) = covice_m(i,j) + covice(i,j)**2 * s + thkice_m(i,j) = thkice_m(i,j) + thkice(i,j)**2 * s + temice_m(i,j) = temice_m(i,j) + temice(i,j)**2 * s + si_um(i,j) = si_um(i,j) + si_u(i,j)**2 * s + si_vm(i,j) = si_vm(i,j) + si_v(i,j)**2 * s + surtxm(i,j) = surtxm(i,j) + surtx(i,j)**2 * s + surtym(i,j) = surtym(i,j) + surty(i,j)**2 * s + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + tracer_bgc_2dm(i,j,ktr) = tracer_bgc_2dm(i,j,ktr) + tracer_bgc_2d(i,j,ktr)**2 *s + enddo !ktr + endif + oneta_m(i,j) = oneta_m(i,j) + oneta(i,j)**2 * s + onetaw_m(i,j) = onetaw_m(i,j) + onetaw(i,j) * s + dp_m(i,j,:) = dp_m(i,j,:) + dp(i,j,:)**2 * s + + sw(:) = onetaw(i,j) * dw(i,j,:) * s + dw_m(i,j,:) = dw_m(i,j,:) + sw(:) + ke_m(i,j,:) = ke_m(i,j,:) + ke(i,j,:)**2 * sw(:) + temp_m(i,j,:) = temp_m(i,j,:) + temp(i,j,:)**2 * sw(:) + saln_m(i,j,:) = saln_m(i,j,:) + saln(i,j,:)**2 * sw(:) + th3d_m(i,j,:) = th3d_m(i,j,:) + th3d(i,j,:)**2 * sw(:) + visc_m(i,j,:) = visc_m(i,j,:) + visc(i,j,:)**2 * sw(:) + tdff_m(i,j,:) = tdff_m(i,j,:) + tdff(i,j,:)**2 * sw(:) + sdff_m(i,j,:) = sdff_m(i,j,:) + sdff(i,j,:)**2 * sw(:) + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,:,ktr) = tracer_m(i,j,:,ktr) + tracer(i,j,:,ktr)**2 * sw(:) + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,:,ktr) = tracer_bgc_3dm(i,j,:,ktr) + tracer_bgc_3d(i,j,:,ktr)**2 * sw(:) + enddo !ktr + endif + endif !ip + enddo + enddo + + ! write(6,*) 'mean_addsq - dp_m = ', dp_m(54, 1,1), + ! & dp( 54, 1,1)**2 + + end subroutine mean_addsq + + subroutine mean_box(a,li,lj,lk,nb,larctic,lperiod) + implicit none + + logical :: larctic,lperiod + integer :: li,lj,lk,nb + real, dimension(li,lj,lk) :: a + + ! --- form the 2*nb+1 square running average of a(:,:,k) + + real, parameter :: spval = 2.0**100 + integer :: i,iq,j,jq,k + real :: rs,qc + real, allocatable :: b(:,:),s(:),q(:) + + allocate( b(1-nb:li+nb,1-nb:lj+nb),s(1-nb:li+nb),q(1-nb:li+nb) ) + + do k= 1,lk + + do j= 1,lj + do i= 1,li + b(i,j) = a(i,j,k) + enddo + enddo + + if (.not.lperiod) then !closed domain + do j = 1,lj + do iq = 1,nb + b( 1-iq,j) = spval + b(li+iq,j) = spval + enddo !iq + enddo !j + do i = 1-nb,li+nb + do jq = 1,nb + b(i, 1-jq) = spval + b(i,lj+jq) = spval + enddo !jq + enddo !i + elseif (.not.larctic) then !lperiod only, i.e. near-global + do jq = 1,nb + do i = 1,li + b(i, 1-jq) = spval !closed bottom boundary + b(i,lj+jq) = spval !closed top boundary + enddo !i + enddo !jq + do j = 1-nb,lj+nb + do iq = 1,nb + b(-nb+iq,j) = b(li-nb+iq,j) !periodic in longitude + b( li+iq,j) = b( iq,j) !periodic in longitude + enddo !iq + enddo !j + else !global with arctic patch + do jq = 1,nb + do i = 1,li + b(i,1-jq) = spval !closed bottom boundary + enddo !i + enddo !jq + do j= lj+1,lj+nb + jq = lj-1-(j-lj) + do i= 1,li + iq = li-mod(i-1,li) + b(i,j) = b(iq,jq) !arctic patch across top boundary + enddo !i + enddo !j + do j= 1-nb,lj+nb + do iq= 1,nb + b(-nb+iq,j) = b(li-nb+iq,j) !periodic in longitude + b( li+iq,j) = b( iq,j) !periodic in longitude + enddo !iq + enddo !j + endif !domain type + + do j= 1,lj + do i= 1-nb,li+nb + rs = 0.0 + qc = 0.0 + do jq= -nb,nb + if (b(i,j+jq).ne.spval) then + rs = rs + b(i,j+jq) + qc = qc + 1.0 + endif + enddo !jq + s(i) = rs + q(i) = qc + enddo !i + do i= 1,li + if (b(i,j) .ne. spval) then + rs = 0.0 + qc = 0.0 + do iq= -nb,nb + rs = rs + s(i+iq) + qc = qc + q(i+iq) + enddo !iq + a(i,j,k) = rs/qc !qc can't be zero, since b(i,j).ne.spval + else + a(i,j,k) = spval + endif + enddo !i + enddo !j + + enddo !k + + deallocate( b,s,q ) + end subroutine mean_box + + subroutine mean_copy + implicit none + + ! --- copy archive to mean archive + + nmean = nstep + + u_m = u + v_m = v + ke_m = ke + temp_m = temp + saln_m = saln + th3d_m = th3d + dp_m = dp + dw_m = dw + visc_m = visc + tdff_m = tdff + sdff_m = sdff + + if (trcout) then + tracer_m = tracer + endif + + if (bgcout) then + tracer_bgc_3dm = tracer_bgc_3d + endif + + ubaro_m = ubaro + vbaro_m = vbaro + pbaro_m = pbaro + kebaro_m = kebaro + montg_m = montg + srfht_m = srfht + steric_m = steric + dpbl_m = dpbl + dpmixl_m = dpmixl + tmix_m = tmix + smix_m = smix + thmix_m = thmix + umix_m = umix + vmix_m = vmix + kemix_m = kemix + surflx_m = surflx + salflx_m = salflx + wtrflx_m = wtrflx + if (icegln) then + covice_m = covice + thkice_m = thkice + temice_m = temice + oneta_m = oneta + onetaw_m = onetaw + endif + if (bgcout) then + tracer_bgc_2dm = tracer_bgc_2d + endif + if (icegln) then + si_um = si_u + si_vm = si_v + surtxm = surtx + surtym = surty + endif + + ! write(6,*) 'mean_copy - dp_m = ', dp_m(54, 1,1), + ! & dp( 54, 1,1) + + end subroutine mean_copy + + subroutine mean_depths + implicit none + + ! --- calculate depthu and depthv + + integer :: i,im,j,jm + + depths(:,:) = 9806.0 * depths(:,:) ! convert to pressure units + + do j= 1,jj + do i= 1,ii + if (i.ne.1) then + im = i-1 + else + im = ii + endif + if (min(ip(i,j),ip(im,j)).eq.1) then + depthu(i,j) = min(depths(i,j),depths(im,j)) + elseif (ip(i ,j).eq.1) then + depthu(i,j) = depths(i ,j) + elseif (ip(im,j).eq.1) then + depthu(i,j) = depths(im,j) + else + depthu(i,j) = 0.0 + endif + + if (j.ne.1) then + jm = j-1 + else + jm = jj + endif + if (min(ip(i,j),ip(i,jm)).eq.1) then + depthv(i,j) = min(depths(i,j),depths(i,jm)) + elseif (ip(i,j) .eq.1) then + depthv(i,j) = depths(i,j) + elseif (ip(i,jm).eq.1) then + depthv(i,j) = depths(i,jm) + else + depthv(i,j) = 0.0 + endif + enddo + enddo + + end subroutine mean_depths + + subroutine mean_diff(nscale,nbox) + implicit none + + integer nscale,nbox + + ! --- form the difference of two archives, 1st already in _m + + real, parameter :: zero = 0.0 + + real :: q + logical :: larctic,lperiod + integer :: i,j,k,ktr + + nmean = 2 !always 2 for diff archives + q = 1.0/real(nscale) + + do j= 1,jj + do i= 1,ii + do k= 1,kk + if (iu(i,j).eq.1) then + u_m(i,j,k) = q*( u_m(i,j,k) - u(i,j,k)) + endif + if (iv(i,j).eq.1) then + v_m(i,j,k) = q*( v_m(i,j,k) - v(i,j,k)) + endif + if (ip(i,j).eq.1) then + dw_m(i,j,k) = dp(i,j,k) + dp_m(i,j,k) = q*( dp_m(i,j,k) - dp(i,j,k)) + ke_m(i,j,k) = q*( ke_m(i,j,k) - ke(i,j,k)) + temp_m(i,j,k) = q*(temp_m(i,j,k) - temp(i,j,k)) + saln_m(i,j,k) = q*(saln_m(i,j,k) - saln(i,j,k)) + th3d_m(i,j,k) = q*(th3d_m(i,j,k) - th3d(i,j,k)) + visc_m(i,j,k) = q*(visc_m(i,j,k) - visc(i,j,k)) + tdff_m(i,j,k) = q*(tdff_m(i,j,k) - tdff(i,j,k)) + sdff_m(i,j,k) = q*(sdff_m(i,j,k) - sdff(i,j,k)) + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,k,ktr) = q*(tracer_m(i,j,k,ktr) - tracer(i,j,k,ktr)) + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,k,ktr) = q*(tracer_bgc_3dm(i,j,k,ktr) - tracer_bgc_3d(i,j,k,ktr)) + enddo !ktr + endif + endif + enddo !k + + if (iu(i,j).eq.1) then + ubaro_m(i,j) = q*( ubaro_m(i,j) - ubaro(i,j)) + umix_m(i,j) = q*( umix_m(i,j) - umix(i,j)) + endif + if (iv(i,j).eq.1) then + vbaro_m(i,j) = q*( vbaro_m(i,j) - vbaro(i,j)) + vmix_m(i,j) = q*( vmix_m(i,j) - vmix(i,j)) + endif + if (ip(i,j).eq.1) then + pbaro_m(i,j) = q*( pbaro_m(i,j) - pbaro(i,j)) + kebaro_m(i,j) = q*(kebaro_m(i,j) - kebaro(i,j)) + montg_m(i,j) = q*( montg_m(i,j) - montg(i,j)) + srfht_m(i,j) = q*( srfht_m(i,j) - srfht(i,j)) + steric_m(i,j) = q*(steric_m(i,j) - steric(i,j)) + dpbl_m(i,j) = q*( dpbl_m(i,j) - dpbl(i,j)) + dpmixl_m(i,j) = q*(dpmixl_m(i,j) - dpmixl(i,j)) + tmix_m(i,j) = q*( tmix_m(i,j) - tmix(i,j)) + smix_m(i,j) = q*( smix_m(i,j) - smix(i,j)) + thmix_m(i,j) = q*( thmix_m(i,j) - thmix(i,j)) + kemix_m(i,j) = q*( kemix_m(i,j) - kemix(i,j)) + surflx_m(i,j) = q*(surflx_m(i,j) - surflx(i,j)) + salflx_m(i,j) = q*(salflx_m(i,j) - salflx(i,j)) + wtrflx_m(i,j) = q*(wtrflx_m(i,j) - wtrflx(i,j)) + if (icegln) then + covice_m(i,j) = q*(covice_m(i,j) - covice(i,j)) + thkice_m(i,j) = q*(thkice_m(i,j) - thkice(i,j)) + temice_m(i,j) = q*(temice_m(i,j) - temice(i,j)) + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + tracer_bgc_2dm(i,j,ktr) = q*(tracer_bgc_2dm(i,j,ktr) - tracer_bgc_2d(i,j,ktr)) + enddo !ktr + endif + if (icegln) then + si_um(i,j) = q*( si_um(i,j) - si_u(i,j)) + si_vm(i,j) = q*( si_vm(i,j) - si_v(i,j)) + surtxm(i,j) = q*( surtxm(i,j) - surtx(i,j)) + surtym(i,j) = q*( surtym(i,j) - surty(i,j)) + endif + endif + enddo !i + enddo !j + + if (nbox.ne.0) then + larctic = any( ip(:,jj).eq.1 ) !land on last row + lperiod = any( ip(ii,:).eq.1 ) !land on last column + + call mean_box( u_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( v_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( temp_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( saln_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( th3d_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( dp_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( ke_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( visc_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( tdff_m,ii,jj,kk,nbox,larctic,lperiod) + call mean_box( sdff_m,ii,jj,kk,nbox,larctic,lperiod) + + if (trcout) then + do ktr= 1,ntracr + call mean_box(tracer_m(1,1,1,ktr),ii,jj,kk,nbox,larctic,lperiod) + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + call mean_box(tracer_bgc_3dm(1,1,1,ktr),ii,jj,kk,nbox,larctic,lperiod) + enddo !ktr + endif + + call mean_box( ubaro_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( umix_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( vbaro_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( vmix_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( pbaro_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(kebaro_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( montg_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( srfht_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(steric_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( dpbl_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(dpmixl_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( tmix_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( smix_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( thmix_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( kemix_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(surflx_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(salflx_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(wtrflx_m,ii,jj, 1,nbox,larctic,lperiod) + if (icegln) then + call mean_box(covice_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(thkice_m,ii,jj, 1,nbox,larctic,lperiod) + call mean_box(temice_m,ii,jj, 1,nbox,larctic,lperiod) + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + call mean_box(tracer_bgc_2dm(1,1,ktr),ii,jj,1,nbox,larctic,lperiod) + enddo !ktr + endif + if (icegln) then + call mean_box( si_um,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( si_vm,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( surtxm,ii,jj, 1,nbox,larctic,lperiod) + call mean_box( surtym,ii,jj, 1,nbox,larctic,lperiod) + endif + endif !nbox + + ! write(6,*) 'mean_diff - dp_m = ', dp_m(54, 1,1), + ! & dw_m(54, 1,1), + ! & dp( 54, 1,1) + + end subroutine mean_diff + + subroutine mean_end + implicit none + +! --- reduce sum of archives to their mean. + + real, parameter :: spval = 2.0**100 + + integer :: i,im,j,jm,k,ktr + real :: s,swk + real, dimension(kk) :: sw(kk) + + s = 1.0/nmean + + do j= 1,jj + do i= 1,ii + if (iu(i,j).eq.1) then + if (i.ne.1) then + im = i-1 + else + im = ii + endif + do k= 1,kk + swk = dpu_m(i,j,k) * s + if (swk.ge.0.000001) then + swk = s/swk + u_m(i,j,k) = u_m(i,j,k) * swk + else ! project into zero thickness layers + u_m(i,j,k) = u_m(i,j,k-1) + endif + enddo + ubaro_m(i,j) = ubaro_m(i,j) * s + umix_m(i,j) = umix_m(i,j) * s + else + u_m(i,j,:) = spval + ubaro_m(i,j) = spval + umix_m(i,j) = spval + endif !iu + + if (iv(i,j).eq.1) then + if (j.ne.1) then + jm = j-1 + else + jm = jj + endif + do k= 1,kk + swk = dpv_m(i,j,k) * s + if (swk.ge.0.000001) then + swk = s/swk + v_m(i,j,k) = v_m(i,j,k) * swk + else ! project into zero thickness layers + v_m(i,j,k) = v_m(i,j,k-1) + endif + enddo + vbaro_m(i,j) = vbaro_m(i,j) * s + vmix_m(i,j) = vmix_m(i,j) * s + else + v_m(i,j,:) = spval + vbaro_m(i,j) = spval + vmix_m(i,j) = spval + endif + + if (ip(i,j).eq.1) then + pbaro_m(i,j) = pbaro_m(i,j) * s + kebaro_m(i,j) = kebaro_m(i,j) * s + montg_m(i,j) = montg_m(i,j) * s + srfht_m(i,j) = srfht_m(i,j) * s + steric_m(i,j) = steric_m(i,j) * s + dpbl_m(i,j) = dpbl_m(i,j) * s + dpmixl_m(i,j) = dpmixl_m(i,j) * s + tmix_m(i,j) = tmix_m(i,j) * s + smix_m(i,j) = smix_m(i,j) * s + thmix_m(i,j) = thmix_m(i,j) * s + kemix_m(i,j) = kemix_m(i,j) * s + surflx_m(i,j) = surflx_m(i,j) * s + salflx_m(i,j) = salflx_m(i,j) * s + wtrflx_m(i,j) = wtrflx_m(i,j) * s + if (icegln) then + covice_m(i,j) = covice_m(i,j) * s + thkice_m(i,j) = thkice_m(i,j) * s + temice_m(i,j) = temice_m(i,j) * s + endif + oneta_m(i,j) = oneta_m(i,j) * s + onetaw_m(i,j) = onetaw_m(i,j) * s + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + tracer_bgc_2dm(i,j,ktr) = tracer_bgc_2dm(i,j,ktr) * s + enddo !ktr + endif + if (icegln) then + si_um(i,j) = si_um(i,j) * s + si_vm(i,j) = si_vm(i,j) * s + surtxm(i,j) = surtxm(i,j) * s + surtym(i,j) = surtym(i,j) * s + endif + + do k= 1,kk + dw_m(i,j,k) = dw_m(i,j,k) * s + dp_m(i,j,k) = dp_m(i,j,k) * s + if (dw_m(i,j,k).ge.0.000001) then + swk = s/dw_m(i,j,k) + ke_m(i,j,k) = ke_m(i,j,k) * swk + temp_m(i,j,k) = temp_m(i,j,k) * swk + saln_m(i,j,k) = saln_m(i,j,k) * swk + th3d_m(i,j,k) = th3d_m(i,j,k) * swk + visc_m(i,j,k) = visc_m(i,j,k) * swk + tdff_m(i,j,k) = tdff_m(i,j,k) * swk + sdff_m(i,j,k) = sdff_m(i,j,k) * swk + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,k,ktr) = tracer_m(i,j,k,ktr) * swk + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,k,ktr) = tracer_bgc_3dm(i,j,k,ktr) * swk + enddo !ktr + endif + else ! project into zero thickness layers + ke_m(i,j,k) = ke_m(i,j,k-1) + temp_m(i,j,k) = temp_m(i,j,k-1) + saln_m(i,j,k) = saln_m(i,j,k-1) + th3d_m(i,j,k) = th3d_m(i,j,k-1) + visc_m(i,j,k) = visc_m(i,j,k-1) + tdff_m(i,j,k) = tdff_m(i,j,k-1) + sdff_m(i,j,k) = sdff_m(i,j,k-1) + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,k,ktr) = tracer_m(i,j,k-1,ktr) + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,k,ktr) = tracer_bgc_3dm(i,j,k-1,ktr) + enddo !ktr + endif + endif + ! --- archived dp_m is based on dp' (dp_m/oneta_m) + dp_m(i,j,k) = dp_m(i,j,k)/onetaw_m(i,j) + dw_m(i,j,k) = dw_m(i,j,k)/onetaw_m(i,j) + p(i,j,k+1) = dw_m(i,j,k) + p(i,j,k) + enddo + else + pbaro_m(i,j) = spval + kebaro_m(i,j) = spval + montg_m(i,j) = spval + srfht_m(i,j) = spval + steric_m(i,j) = spval + dpbl_m(i,j) = spval + dpmixl_m(i,j) = spval + tmix_m(i,j) = spval + smix_m(i,j) = spval + thmix_m(i,j) = spval + kemix_m(i,j) = spval + surflx_m(i,j) = spval + salflx_m(i,j) = spval + wtrflx_m(i,j) = spval + if (icegln) then + covice_m(i,j) = spval + thkice_m(i,j) = spval + temice_m(i,j) = spval + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + tracer_bgc_2dm(i,j,ktr) = spval + enddo !ktr + endif + if (icegln) then + si_um(i,j) = spval + si_vm(i,j) = spval + surtxm(i,j) = spval + surtym(i,j) = spval + endif + oneta_m(i,j) = spval + onetaw_m(i,j) = spval + + dw_m(i,j,:) = spval + dp_m(i,j,:) = spval + ke_m(i,j,:) = spval + temp_m(i,j,:) = spval + saln_m(i,j,:) = spval + th3d_m(i,j,:) = spval + visc_m(i,j,:) = spval + tdff_m(i,j,:) = spval + sdff_m(i,j,:) = spval + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,:,ktr) = spval + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,:,ktr) = spval + enddo !ktr + endif + endif + enddo + enddo + + ! write(6,*) 'mean_end - dp_m = ', dp_m(54, 1,1) + + end subroutine mean_end + + subroutine mean_std + implicit none + + ! --- form the std.dev = sqrt(mnsq-mean**2) + + real, parameter :: zero = 0.0 + integer :: i,j,k,ktr + real :: std,x + + std(x) = sqrt(max(zero,x)) + + do j= 1,jj + do i= 1,ii + do k= 1,kk + if (iu(i,j).eq.1) then + u_m(i,j,k) = std( u(i,j,k) - u_m(i,j,k)**2) + endif + if (iv(i,j).eq.1) then + v_m(i,j,k) = std( v(i,j,k) - v_m(i,j,k)**2) + endif + if (ip(i,j).eq.1) then + dw_m(i,j,k) = dp_m(i,j,k) + dp_m(i,j,k) = std( dp(i,j,k) - dp_m(i,j,k)**2) + ke_m(i,j,k) = std( ke(i,j,k) - ke_m(i,j,k)**2) + temp_m(i,j,k) = std(temp(i,j,k) - temp_m(i,j,k)**2) + saln_m(i,j,k) = std(saln(i,j,k) - saln_m(i,j,k)**2) + th3d_m(i,j,k) = std(th3d(i,j,k) - th3d_m(i,j,k)**2) + visc_m(i,j,k) = std(visc(i,j,k) - visc_m(i,j,k)**2) + tdff_m(i,j,k) = std(tdff(i,j,k) - tdff_m(i,j,k)**2) + sdff_m(i,j,k) = std(sdff(i,j,k) - sdff_m(i,j,k)**2) + if (trcout) then + do ktr= 1,ntracr + tracer_m(i,j,k,ktr) = std(tracer(i,j,k,ktr) - tracer_m(i,j,k,ktr)**2) + enddo !ktr + endif + if (bgcout) then + do ktr= 1,ntracr_bgc_3d + tracer_bgc_3dm(i,j,k,ktr) = std(tracer_bgc_3d(i,j,k,ktr) - tracer_bgc_3dm(i,j,k,ktr)**2) + enddo !ktr + endif + endif + enddo + + if (iu(i,j).eq.1) then + ubaro_m(i,j) = std( ubaro(i,j) - ubaro_m(i,j)**2) + umix_m(i,j) = std( umix(i,j) - umix_m(i,j)**2) + endif + if (iv(i,j).eq.1) then + vbaro_m(i,j) = std( vbaro(i,j) - vbaro_m(i,j)**2) + vmix_m(i,j) = std( vmix(i,j) - vmix_m(i,j)**2) + endif + if (ip(i,j).eq.1) then + pbaro_m(i,j) = std( pbaro(i,j) - pbaro_m(i,j)**2) + kebaro_m(i,j) = std(kebaro(i,j) - kebaro_m(i,j)**2) + montg_m(i,j) = std( montg(i,j) - montg_m(i,j)**2) + srfht_m(i,j) = std( srfht(i,j) - srfht_m(i,j)**2) + steric_m(i,j) = std(steric(i,j) - steric_m(i,j)**2) + dpbl_m(i,j) = std( dpbl(i,j) - dpbl_m(i,j)**2) + dpmixl_m(i,j) = std(dpmixl(i,j) - dpmixl_m(i,j)**2) + tmix_m(i,j) = std( tmix(i,j) - tmix_m(i,j)**2) + smix_m(i,j) = std( smix(i,j) - smix_m(i,j)**2) + thmix_m(i,j) = std( thmix(i,j) - thmix_m(i,j)**2) + kemix_m(i,j) = std( kemix(i,j) - kemix_m(i,j)**2) + surflx_m(i,j) = std(surflx(i,j) - surflx_m(i,j)**2) + salflx_m(i,j) = std(salflx(i,j) - salflx_m(i,j)**2) + wtrflx_m(i,j) = std(wtrflx(i,j) - wtrflx_m(i,j)**2) + if (icegln) then + covice_m(i,j) = std(covice(i,j) - covice_m(i,j)**2) + thkice_m(i,j) = std(thkice(i,j) - thkice_m(i,j)**2) + temice_m(i,j) = std(temice(i,j) - temice_m(i,j)**2) + endif + onetaw_m(i,j) = onetaw_m(i,j) + oneta_m(i,j) = std( oneta(i,j) - oneta_m(i,j)**2) + if (bgcout) then + do ktr= 1,ntracr_bgc_2d + tracer_bgc_2dm(i,j,ktr) = std(tracer_bgc_2d(i,j,ktr) - tracer_bgc_2dm(i,j,ktr)**2) + enddo !ktr + endif + if (icegln) then + si_um(i,j) = std( si_u(i,j) - si_um(i,j)**2) + si_vm(i,j) = std( si_v(i,j) - si_vm(i,j)**2) + surtxm(i,j) = std( surtx(i,j) - surtxm(i,j)**2) + surtym(i,j) = std( surty(i,j) - surtym(i,j)**2) + endif + endif + enddo + enddo + + ! write(6,*) 'mean_std - dp_m = ', dp_m(54, 1,1), + ! & dw_m(54, 1,1), + ! & dp( 54, 1,1) + + end subroutine mean_std + + subroutine mean_velocity + implicit none + + ! --- update velocity to include depth averaged component, and + ! --- calculate kinetic energy. + ! --- only called for standard archive fields. + + integer :: i,ia,ip1,j + + do j= 1,jj + do i= 1,ii + if (iu(i,j).eq.1) then + u(i,j,:) = u(i,j,:) + ubaro(i,j) + umix(i,j) = umix(i,j) + ubaro(i,j) + endif + if (iv(i,j).eq.1) then + v(i,j,:) = v(i,j,:) + vbaro(i,j) + vmix(i,j) = vmix(i,j) + vbaro(i,j) + endif + enddo + enddo + + do j= 1,jj-1 + do i= 1,ii + if (i.ne.ii) then + ip1 = i+1 + else + ip1 = 1 !global periodic region, also works for closed domains since ip(ii,:)=0 + endif + if (ip(i,j).eq.1) then + ! --- kinetic energy / mass (m**2/s**2) + ke(i,j,:) = 0.5*((0.5*( u(i,j,:) + u(ip1,j,:)))**2 + & + (0.5*( v(i,j,:) + v(i,j+1,:)))**2 ) + kemix(i,j) = 0.5*((0.5*( umix(i,j) + umix(ip1,j) ))**2 + & + (0.5*( vmix(i,j) + vmix(i,j+1) ))**2 ) + kebaro(i,j) = 0.5*((0.5*(ubaro(i,j) + ubaro(ip1,j) ))**2 + & + (0.5*(vbaro(i,j) + vbaro(i,j+1) ))**2 ) + endif + enddo + enddo +! --- arctic patch, also works for closed domains since ip(:,jj)=0 + do i= 1,ii + ia = ii-mod(i-1,ii) + if (ip(i,jj).eq.1) then + ke(i,jj,:) = ke(ia,jj-1,:) + kemix(i,jj) = kemix(ia,jj-1) + kebaro(i,jj) = kebaro(ia,jj-1) + endif + enddo !i + + ! write(6,*) 'mean_velocity - ke = ', ke(54, 1,1) + + end subroutine mean_velocity + +end module mod_mean + diff --git a/hycom/MSCPROGS/src/Average/putdat.F90 b/hycom/MSCPROGS/src/Average/putdat.F90 new file mode 100644 index 00000000..3b4968aa --- /dev/null +++ b/hycom/MSCPROGS/src/Average/putdat.F90 @@ -0,0 +1,340 @@ +subroutine putdat(flnm,time_min,time_max,time,mntype,iexpt,jexpt,yrflag,thbase) + use mod_mean ! HYCOM mean array interface + use mod_za ! HYCOM array I/O interface + + character(len=*) :: flnm + real*8 :: time_min,time_max,time + real :: thbase + integer :: mntype,iexpt,jexpt,yrflag + !logical :: icegln,trcout + + ! --- write mean or mean squared or std or diff model fields. + ! --- HYCOM 2.0 array I/O archive file. + + real :: coord,xmin,xmax + integer :: i,j,k,ktr,iversn,l,nop + + data nop/24/ + + l = len_trim(flnm) + open (unit=nop,file=flnm(1:l)//'.b',form='formatted',status='new',action='write') + call zaiopf(flnm(1:l)//'.a','new', nop) + + ! --- header. + + iversn = 20 + if (mntype.eq.4) then + write(nop,115) ctitle,iversn,jexpt,iexpt,yrflag,idm,jdm + write( lp, *) + write( lp,115) ctitle,iversn,jexpt,iexpt,yrflag,idm,jdm + else + write(nop,116) ctitle,iversn,iexpt,yrflag,idm,jdm + write( lp, *) + write( lp,116) ctitle,iversn,iexpt,yrflag,idm,jdm + endif + + 115 format(a80/a80/a80/a80/ & + ,i5,4x,'''iversn'' = hycom version number x10' / & + ,i5,4x,'''jexpt '' = 1st experiment number x10'/ & + ,i5,4x,'''iexpt '' = 2nd experiment number x10'/ & + ,i5,4x,'''yrflag'' = days in year flag' / & + ,i5,4x,'''idm '' = longitudinal array size' / & + ,i5,4x,'''jdm '' = latitudinal array size' ) + + 116 format(a80/a80/a80/a80/ & + ,i5,4x,'''iversn'' = hycom version number x10' / & + ,i5,4x,'''iexpt '' = experiment number x10' / & + ,i5,4x,'''yrflag'' = days in year flag' / & + ,i5,4x,'''idm '' = longitudinal array size' / & + ,i5,4x,'''jdm '' = latitudinal array size') + + 118 format('field no. recs ' & + ,a4,' day',' k dens min max') + + if (mntype.eq.1) then + write(nop,118) 'mean' + write( lp,118) 'mean' + elseif (mntype.eq.2) then + write(nop,118) 'mnsq' + write( lp,118) 'mnsq' + elseif (mntype.eq.3) then + write(nop,118) 'std.' + write( lp,118) 'std.' + elseif (mntype.eq.4) then + write(nop,118) 'diff' + write( lp,118) 'diff' + else + write(lp,'(/a,i5/)') 'error in putdat - illegal mntype',mntype + stop + endif + + nmean = nstep_m + + ! --- surface fields. + + coord=0.0 + + call zaiowr(montg_m,ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'montg1 ',nmean,time_min,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'montg1 ',nmean,time_min,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(srfht_m,ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'srfhgt ',nmean,time_max,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'srfhgt ',nmean,time_max,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(surflx_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'surflx ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'surflx ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(salflx_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'salflx ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'salflx ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(dpbl_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'bl_dpth ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'bl_dpth ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(dpmixl_m,ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'mix_dpth',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'mix_dpth',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(tmix_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'tmix ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'tmix ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(smix_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'smix ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'smix ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(thmix_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'thmix ',nmean,time,0,thbase,xmin,xmax + call flush(nop) + write ( lp,117) 'thmix ',nmean,time,0,thbase,xmin,xmax + call flush( lp) + + call zaiowr(umix_m,iu,.true., xmin,xmax, nop, .false.) + write (nop,117) 'umix ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'umix ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(vmix_m,iv,.true., xmin,xmax, nop, .false.) + write (nop,117) 'vmix ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'vmix ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(kemix_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'kemix ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'kemix ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + if (icegln) then + call zaiowr(covice_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'covice ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'covice ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(thkice_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'thkice ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'thkice ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(temice_m,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'temice ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'temice ',nmean,time,0,coord,xmin,xmax + call flush( lp) + endif + + ! --- depth averaged fields + + call zaiowr(ubaro_m,iu,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'u_btrop ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'u_btrop ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(vbaro_m,iv,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'v_btrop ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'v_btrop ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + call zaiowr(kebaro_m,ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'kebtrop ',nmean,time,0,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'kebtrop ',nmean,time,0,coord,xmin,xmax + call flush( lp) + + ! --- layer loop. + + do k=1, kk + coord=theta(k) + + call zaiowr(u_m(1,1,k),iu,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'u-vel. ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'u-vel. ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + call zaiowr(v_m(1,1,k),iv,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'v-vel. ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'v-vel. ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + call zaiowr(ke_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'k.e. ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'k.e. ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + if (mntype.eq.3 .or. mntype.eq.4) then + call zaiowr(dw_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'mnthknss',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'mnthknss',nmean,time,k,coord,xmin,xmax + call flush( lp) + endif + + call zaiowr(dp_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'thknss ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'thknss ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + if (mntype.eq.2) then + call zaiowr(dw_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'mnthknss',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'mnthknss',nmean,time,k,coord,xmin,xmax + call flush( lp) + endif + + call zaiowr(temp_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'temp ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'temp ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + call zaiowr(saln_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'salin ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'salin ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + call zaiowr(th3d_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'density ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'density ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + if (trcout) then + do ktr= 1,ntracr + call zaiowr(tracer_m(1,1,k,ktr),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'tracer ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'tracer ',nmean,time,k,coord,xmin,xmax + call flush( lp) + enddo !ktr + endif + + call zaiowr(visc_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 'viscty ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 'viscty ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + call zaiowr(tdff_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 't-diff ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 't-diff ',nmean,time,k,coord,xmin,xmax + call flush( lp) + + call zaiowr(sdff_m(1,1,k),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) 's-diff ',nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) 's-diff ',nmean,time,k,coord,xmin,xmax + call flush( lp) + enddo !k + + ! --- ECOSMO + + if (bgcout) then + ! --- ECOSMO 3D loop. + do k = 1, kk + coord=theta(k) + do ktr= 1,ntracr_bgc_3d + call zaiowr(tracer_bgc_3dm(1,1,k,ktr),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) nvar_bgc_3d(ktr),nmean,time,k,coord,xmin,xmax + call flush(nop) + write ( lp,117) nvar_bgc_3d(ktr),nmean,time,k,coord,xmin,xmax + call flush( lp) + enddo !ktr + enddo !k + ! --- ECOSMO 2D loop + do ktr= 1,ntracr_bgc_2d + call zaiowr(tracer_bgc_2dm(1,1,ktr),ip,.true.,xmin,xmax, nop, .false.) + write (nop,117) nvar_bgc_2d(ktr),nmean,time,0,0.0,xmin,xmax + call flush(nop) + write ( lp,117) nvar_bgc_2d(ktr),nmean,time,0,0.0,xmin,xmax + call flush( lp) + enddo !ktr + endif ! bgcout + + ! --- ice 2D + + if (icegln .and. mntype.eq.1) then + call zaiowr(si_um,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'si_u ',nmean,time,0,0.0,xmin,xmax + call flush(nop) + write ( lp,117) 'si_u ',nmean,time,0,0.0,xmin,xmax + call flush( lp) + + call zaiowr(si_vm,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'si_v ',nmean,time,0,0.0,xmin,xmax + call flush(nop) + write ( lp,117) 'si_v ',nmean,time,0,0.0,xmin,xmax + call flush( lp) + + call zaiowr(surtxm,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'surtx ',nmean,time,0,0.0,xmin,xmax + call flush(nop) + write ( lp,117) 'surtx ',nmean,time,0,0.0,xmin,xmax + call flush( lp) + + call zaiowr(surtym,ip,.true., xmin,xmax, nop, .false.) + write (nop,117) 'surty ',nmean,time,0,0.0,xmin,xmax + call flush(nop) + write ( lp,117) 'surty ',nmean,time,0,0.0,xmin,xmax + call flush( lp) + endif + + 117 format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7) + + close (unit=nop) + call zaiocl(nop) + + return + +end subroutine putdat diff --git a/hycom/MSCPROGS/src/IceDriftNC/icedrift_osisafNC b/hycom/MSCPROGS/src/IceDriftNC/icedrift_osisafNC new file mode 100755 index 00000000..276daa70 Binary files /dev/null and b/hycom/MSCPROGS/src/IceDriftNC/icedrift_osisafNC differ diff --git a/hycom/MSCPROGS/src/Nersclib/core.28559 b/hycom/MSCPROGS/src/Nersclib/core.28559 deleted file mode 100644 index 2f6c0a37..00000000 Binary files a/hycom/MSCPROGS/src/Nersclib/core.28559 and /dev/null differ