Commit a707c8b1 authored by HOURDIN Christophe's avatar HOURDIN Christophe
Browse files

Minor adds

parent c25fa193
! $Id$
!
!=========================================================================
! ROMS_AGRIF is a branch of ROMS developped at IRD and INRIA, in France.
! The two other branches, from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al), are under MIT/X style license.
! ROMS_AGRIF specific routines (nesting) are under CeCILL-C license.
!
! ROMS_AGRIF website : http://www.romsagrif.org
!=========================================================================
!
subroutine biology_tile(Istr,Iend,Jstr,Jend)
!------------------------------------------------------------------
!
! ROUTINE biology_pisces : PISCES MODEL
! *************************************
!
!
! PURPOSE.
! --------
! *ROMS_PISCES ECOSYSTEM MODEL FOR THE WHOLE OCEAN
! THIS ROUTINE COMPUTES INTERACTIONS
! BETWEEN THE DIFFERENT COMPARTMENTS OF THE
! MODEL
!----------------------------------------------------------------
USE sms_pisces
USE trcsms_pisces
C implicit none
INTEGER Istr,Iend,Jstr,Jend
INTEGER i, j, jk, k, jn
REAL zdiag(GLOBAL_2D_ARRAY,0:N)
#include "ocean2pisces.h90"
CALL ocean_2_pisces( Istr,Iend,Jstr,Jend)
! modif SPOUS ASAP
#ifdef AGRIF
if (AGRIF_Root()) then
DO jk = KRANGE
DO j = JRANGE
DO i = IRANGE
tsn(i,j,K,jp_sal) =
& MAX(0. , tsn(i,j,K,jp_sal) )
ENDDO
ENDDO
ENDDO
endif
#endif
! modif SPOUS ASAP
DO jn = 1, jptra
DO jk = KRANGE
DO j = JRANGE
DO i = IRANGE ! masked grid volume
trb(i,j,K,jn)
! MODIF ASAP SPOUS
!! & = MAX( 0., trb(i,j,K,jn) * 1.e-6 )
& = MAX( 1.e-14, trb(i,j,K,jn) * 1.e-6 )
! MODIF ASAP SPOUS
ENDDO
ENDDO
ENDDO
END DO
DO jk = KRANGE
DO j = JRANGE
DO i = IRANGE ! masked grid volume
trb(i,j,K,jpno3)
& = trb(i,j,K,jpno3) / rno3
trb(i,j,K,jpnh4)
& = trb(i,j,K,jpnh4) / rno3
trb(i,j,K,jppo4)
& = trb(i,j,K,jppo4) / po4r
ENDDO
ENDDO
ENDDO
CALL trc_sms_pisces( iic )
DO jn = 1, jptra
DO jk = KRANGE
DO j = JRANGE
DO i = IRANGE ! masked grid volume
trb(i,j,K,jn)
& = trb(i,j,K,jn) * 1.e6
ENDDO
ENDDO
ENDDO
END DO
DO jk = KRANGE
DO j = JRANGE
DO i = IRANGE ! masked grid volume
trb(i,j,K,jpno3)
& = trb(i,j,K,jpno3) * rno3
trb(i,j,K,jpnh4)
& = trb(i,j,K,jpnh4) * rno3
trb(i,j,K,jppo4)
& = trb(i,j,K,jppo4) * po4r
ENDDO
ENDDO
ENDDO
RETURN
END
! $Id: hmix_coef.F 1618 2014-12-18 14:39:51Z rblod $
!
!======================================================================
! CROCO is a branch of ROMS developped at IRD and INRIA, in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#if defined VIS_COEF_3D || (!defined SOLVE3D && defined UV_VIS_SMAGO)
subroutine hvisc_coef (tile)
implicit none
integer tile, trd, omp_get_thread_num
# include "param.h"
# include "private_scratch.h"
# include "compute_tile_bounds.h"
trd=omp_get_thread_num()
call hvisc_coef_tile (Istr,Iend,Jstr,Jend, A3d(1,1,trd))
return
end
!
subroutine hvisc_coef_tile (Istr,Iend,Jstr,Jend, defrate)
!
!==================================================================
!
! Compute horizontal viscosity
!
! Patrick Marchesiello and Pierrick Penven, IRD 2007
!
!==================================================================
!
implicit none
# include "param.h"
# include "grid.h"
# include "mixing.h"
# include "ocean3d.h"
# include "ocean2d.h"
# include "scalars.h"
integer Istr,Iend,Jstr,Jend, i,j,k,itrc,
& imin,imax,jmin,jmax, kp
real horcon, cff,cff0,cff1,cff2,cff3,cff4,cff5
# ifdef UV_VIS_SMAGO_3D
real Ric, Pr
parameter (horcon=0.05, Ric=0.25, Pr=0.3)
# else
parameter (horcon=0.025)
# endif
# ifdef SWASH
real dZdt,dZdt0,Bbr
# endif
real defrate(PRIVATE_2D_SCRATCH_ARRAY,1:N)
!
# include "compute_auxiliary_bounds.h"
!
# ifndef EW_PERIODIC
if (WESTERN_EDGE) then
imin=istr
else
imin=istr-1
endif
if (EASTERN_EDGE) then
imax=iend
else
imax=iend+1
endif
# else
imin=istr-1
imax=iend+1
# endif
# ifndef NS_PERIODIC
if (SOUTHERN_EDGE) then
jmin=jstr
else
jmin=jstr-1
endif
if (NORTHERN_EDGE) then
jmax=jend
else
jmax=jend+1
endif
# else
jmin=jstr-1
jmax=jend+1
# endif
!
# if defined UV_VIS2 && defined GLS_MIXING_3D
!---------------------------------------------------------------
! GLS horizontal viscosity at RHO points
!
do j=jmin,jmax
do i=imin,imax
cff0=om_r(i,j)*on_r(i,j)
cff1=0.01*cff0/dt
do k=1,N
cff2=0.5*(Akv(i,j,k)+Akv(i,j,k-1))
visc3d_r(i,j,k)=min(cff1,cff2)
enddo
enddo
enddo
# elif defined UV_VIS2 && defined UV_VIS_SMAGO
!---------------------------------------------------------------
! Smagorinsky viscosity at RHO points
!
! A = Cs L^2 D
!
! with deformation rate D:
!
! D = sqrt[ 2 DijDij ] -- Dij = 0.5*(dui/dxj + duj/dxi) --
! = sqrt[ 2 dudx^2 + 2 dvdy^2 + (dvdx+dudy)^2 --> 2D
! + 2 dwdz^2 + (dudz+dwdx)^2 + (dvdz+dwdy)^2 ] --> 3D
!
! L Length-scale: horizontal Lh = sqrt(dx*dy) --> visc3d
! vertical Lv = Hz --> Akv
!
! For LES problems, Lilly's buoyancy correction is applied based
! on the ratio of gradient and critical Ri (mixing goes to 0 for
! Rig > Ric)
!---------------------------------------------------------------
!
! Horizontal deformation rate (squared & halved) at RHO point:
!
! dudx^2 + dvdy^2 + 0.5(dvdx+dudy)^2
!
do k=1,N
do j=jmin,jmax
do i=imin,imax
defrate(i,j,k)=
& ((u(i+1,j,k,nrhs)-u(i,j,k,nrhs))*pm(i,j))**2
& +((v(i,j+1,k,nrhs)-v(i,j,k,nrhs))*pn(i,j))**2
& +0.125*( pn(i,j)*
& (u(i,j+1,k,nrhs)+u(i+1,j+1,k,nrhs)
& -u(i,j-1,k,nrhs)-u(i+1,j-1,k,nrhs))
& +pm(i,j)*
& (v(i+1,j,k,nrhs)+v(i+1,j+1,k,nrhs)
& -v(i-1,j,k,nrhs)-v(i-1,j+1,k,nrhs)) )**2
enddo
enddo
enddo
!
# ifdef UV_VIS_SMAGO_3D
!
! Add vertical deformation rate
!
! dwdz^2 + 0.5*(dudz+dwdx)^2 + 0.5*(dvdz+dwdy)^2
!
do k=2,N-1
do j=jmin,jmax
do i=imin,imax
defrate(i,j,k)=defrate(i,j,k)
# ifdef NBQ
& +((wz(i,j,k,nrhs)-wz(i,j,k-1,nrhs))/Hzr(i,j,k))**2
# endif
& +0.125*( (u(i,j,k+1,nrhs)+u(i+1,j,k+1,nrhs)
& -u(i,j,k-1,nrhs)-u(i+1,j,k-1,nrhs))
& /(z_r(i,j,k+1)-z_r(i,j,k-1))
# ifdef NBQ
& +pm(i,j)*
& (wz(i+1,j,k,nrhs)+wz(i+1,j,k+1,nrhs)
& -wz(i-1,j,k,nrhs)-wz(i-1,j,k+1,nrhs))
# endif
& )**2
& +0.125*( (v(i,j,k+1,nrhs)+v(i,j+1,k+1,nrhs)
& -v(i,j,k-1,nrhs)-v(i,j+1,k-1,nrhs))
& /(z_r(i,j,k+1)-z_r(i,j,k-1))
# ifdef NBQ
& +pn(i,j)*
& (wz(i,j+1,k,nrhs)+wz(i,j+1,k+1,nrhs)
& -wz(i,j-1,k,nrhs)-wz(i,j-1,k+1,nrhs))
# endif
& )**2
enddo
enddo
enddo
do j=jmin,jmax
do i=imin,imax
defrate(i,j,N)=max(1.5*defrate(i,j,N-1)-0.5*defrate(i,j,N-2),
& 0.0)
enddo
enddo
# endif
!
! Compute lateral Smagorinsky viscosity at RHO points
!
! Lh = sqrt(dx*dy) --> visc3d
!
do j=jmin,jmax
do i=imin,imax
cff=horcon*om_r(i,j)*on_r(i,j)
do k=1,N
cff1=2.*defrate(i,j,k) ! D²
# ifdef UV_VIS_SMAGO_3D
! Buoyancy correction (Lilly, 1962)
cff2=sqrt(max(0.,1.-(bvf(i,j,k)/max(cff1,1.e-12))/Ric))
visc3d_r(i,j,k)=cff*sqrt(cff1)*cff2
# else
visc3d_r(i,j,k)=cff*sqrt(cff1) ! Cs*dxdy*D²
# endif
enddo
enddo
enddo
# ifdef UV_VIS_SMAGO_3D
!
! Compute vertical Smagorinsky viscosity at W points
!
! Lv = Hz --> Akv
!
# ifdef SWASH0
do j=jmin,jmax
do i=imin,imax
cff=z_w(i,j,N)-z_w(i,j,0)
dZdt=abs(zeta(i,j,knew)-zeta(i,j,kstp))/dtfast
dZdt0=0.01*sqrt(g*cff)
Bbr=max(0,INT(min(1.,dZdt/dZdt0-5.)))
cff0=om_r(i,j)*on_r(i,j)
cff1=0.01*cff0
cff2=0.01*cff0/dt
do k=1,N
cff3=cff1*sqrt(2.*defrate(i,j,k))
cff3=cff3*(1.+10.*Bbr*cff*dZdt)
visc3d_r(i,j,k)=min(cff2,cff3)
enddo
enddo
enddo
# endif
do k=1,N
kp=min(k+1,N)
do j=jmin,jmax
do i=imin,imax
Akv(i,j,k)=0.5*(visc3d_r(i,j,k)+visc3d_r(i,j,kp)) ! Cs*dz²*D²
& *pm(i,j)*pn(i,j)
& *Hz(i,j,k)*Hz(i,j,k)
# if !defined LMD_MIXING && !defined BVF_MIXING \
&& !defined GLS_MIXING
& + Akv_bak
# endif
!
! Apply SGS Prandtl number for tracers
!
Akt(i,j,k,itemp)=Pr*Akv(i,j,k)
# ifdef SALINITY
Akt(i,j,k,isalt)=Akt(i,j,k,itemp)
# endif
enddo
enddo
enddo
# endif /* UV_VIS_SMAGO_3D */
!
!------------------------------------------------------------
! Boundary conditions
!------------------------------------------------------------
!
# ifndef EW_PERIODIC
if (WESTERN_EDGE) then
do k=1,N
do j=JstrV-1,Jend
visc3d_r(Istr-1,j,k)=visc3d_r(Istr,j,k)
enddo
enddo
endif
if (EASTERN_EDGE) then
do k=1,N
do j=JstrV-1,Jend
visc3d_r(Iend+1,j,k)=visc3d_r(Iend,j,k)
enddo
enddo
endif
# endif /* !EW_PERIODIC */
# ifndef NS_PERIODIC
if (SOUTHERN_EDGE) then
do k=1,N
do i=IstrU-1,Iend
visc3d_r(i,Jstr-1,k)=visc3d_r(i,Jstr,k)
enddo
enddo
endif
if (NORTHERN_EDGE) then
do k=1,N
do i=IstrU-1,Iend
visc3d_r(i,Jend+1,k)=visc3d_r(i,Jend,k)
enddo
enddo
endif
# endif /* !NS_PERIODIC */
!
! Corners
!
# if !defined EW_PERIODIC && !defined NS_PERIODIC
if (SOUTHERN_EDGE .and. WESTERN_EDGE) then
do k=1,N
visc3d_r(Istr-1,Jstr-1,k)=0.5*
& ( visc3d_r(Istr,Jstr-1,k)
& +visc3d_r(Istr-1,Jstr,k))
enddo
endif
if (SOUTHERN_EDGE .and. EASTERN_EDGE) then
do k=1,N
visc3d_r(Iend+1,Jstr-1,k)=0.5*
& (visc3d_r(Iend,Jstr-1,k)
& +visc3d_r(Iend+1,Jstr,k))
enddo
endif
if (NORTHERN_EDGE .and. WESTERN_EDGE) then
do k=1,N
visc3d_r(Istr-1,Jend+1,k)=0.5*
& ( visc3d_r(Istr,Jend+1,k)
& +visc3d_r(Istr-1,Jend,k))
enddo
endif
if (NORTHERN_EDGE .and. EASTERN_EDGE) then
do k=1,N
visc3d_r(Iend+1,Jend+1,k)=0.5*
& ( visc3d_r(Iend,Jend+1,k)
& +visc3d_r(Iend+1,Jend,k))
enddo
endif
# endif /* !EW_PERIODIC && !NS_PERIODIC */
!
# if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
! call exchange_r3d_tile (Istr,Iend,Jstr,Jend,visc3d_r)
# endif
!
! Viscosity at PSI points
!
do k=1,N
do j=Jstr,JendR
do i=Istr,IendR
visc3d_p(i,j,k)=0.25*
& ( visc3d_r(i,j ,k)+visc3d_r(i-1,j ,k)
& +visc3d_r(i,j-1,k)+visc3d_r(i-1,j-1,k))
enddo
enddo
enddo
!
# if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
call exchange_p3d_tile(Istr,Iend,Jstr,Jend,visc3d_p)
# endif
# ifdef UV_VIS_SMAGO_3D
!------------------------------------------------------
! Lateral boundary conditions for vertical diffusivity
! in case of 3D Smagorinsky SGS model
!------------------------------------------------------
!
# define k0 1
# ifndef EW_PERIODIC
if (WESTERN_EDGE) then
do j=jstr,jend
do k=k0,N
Akv(istr-1,j,k)=Akv(istr,j,k)
Akt(istr-1,j,k,itemp)=Akt(istr,j,k,itemp)
# ifdef SALINITY
Akt(istr-1,j,k,isalt)=Akt(istr,j,k,isalt)
# endif
enddo
enddo
endif
if (EASTERN_EDGE) then
do j=jstr,jend
do k=k0,N
Akv(iend+1,j,k)=Akv(iend,j,k)
Akt(iend+1,j,k,itemp)=Akt(iend,j,k,itemp)
# ifdef SALINITY
Akt(iend+1,j,k,isalt)=Akt(iend,j,k,isalt)
# endif
enddo
enddo
endif
# endif
# ifndef NS_PERIODIC
if (SOUTHERN_EDGE) then
do i=istr,iend
do k=k0,N
Akv(i,jstr-1,k)=Akv(i,jstr,k)
Akt(i,jstr-1,k,itemp)=Akt(i,jstr,k,itemp)
# ifdef SALINITY
Akt(i,jstr-1,k,isalt)=Akt(i,jstr,k,isalt)
# endif
enddo
enddo
endif
if (NORTHERN_EDGE) then
do i=istr,iend
do k=k0,N
Akv(i,jend+1,k)=Akv(i,jend,k)
Akt(i,jend+1,k,itemp)=Akt(i,jend,k,itemp)
# ifdef SALINITY
Akt(i,jend+1,k,isalt)=Akt(i,jend,k,isalt)
# endif
enddo
enddo
endif
# ifndef EW_PERIODIC
if (WESTERN_EDGE .and. SOUTHERN_EDGE) then
do k=k0,N
Akv(istr-1,jstr-1,k)=Akv(istr,jstr,k)
Akt(istr-1,jstr-1,k,itemp)=Akt(istr,jstr,k,itemp)
# ifdef SALINITY
Akt(istr-1,jstr-1,k,isalt)=Akt(istr,jstr,k,isalt)
# endif
enddo
endif
if (WESTERN_EDGE .and. NORTHERN_EDGE) then
do k=k0,N
Akv(istr-1,jend+1,k)=Akv(istr,jend,k)
Akt(istr-1,jend+1,k,itemp)=Akt(istr,jend,k,itemp)
# ifdef SALINITY
Akt(istr-1,jend+1,k,isalt)=Akt(istr,jend,k,isalt)
# endif
enddo
endif
if (EASTERN_EDGE .and. SOUTHERN_EDGE) then
do k=k0,N
Akv(iend+1,jstr-1,k)=Akv(iend,jstr,k)
Akt(iend+1,jstr-1,k,itemp)=Akt(iend,jstr,k,itemp)
# ifdef SALINITY
Akt(iend+1,jstr-1,k,isalt)=Akt(iend,jstr,k,isalt)
# endif
enddo
endif
if (EASTERN_EDGE .and. NORTHERN_EDGE) then
do k=k0,N
Akv(iend+1,jend+1,k)=Akv(iend,jend,k)
Akt(iend+1,jend+1,k,itemp)=Akt(iend,jend,k,itemp)
# ifdef SALINITY
Akt(iend+1,jend+1,k,isalt)=Akt(iend,jend,k,isalt)
# endif
enddo
endif
# endif
# endif
!# if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
! call exchange_w3d_tile (istr,iend,jstr,jend, Akv)
! call exchange_w3d_tile (istr,iend,jstr,jend,
! & Akt(START_2D_ARRAY,0,itemp))
!# ifdef SALINITY
! call exchange_w3d_tile (istr,iend,jstr,jend,
! & Akt(START_2D_ARRAY,0,isalt))
!# endif
!# endif
# endif /* UV_VIS_SMAGO_3D */
# endif /* UV_VIS_SMAGO */
#else
subroutine hvisc_coef_empty
#endif
return
end
!
!===================================================================
! DIFFUSIVITY
!===================================================================
!
#ifdef DIF_COEF_3D
subroutine hdiff_coef (tile)
implicit none
integer tile, trd, omp_get_thread_num
# include "param.h"
# include "private_scratch.h"
# include "compute_tile_bounds.h"
trd=omp_get_thread_num()
call hdiff_coef_tile (Istr,Iend,Jstr,Jend,A2d(1,1,trd))
return
end
!
subroutine hdiff_coef_tile (Istr,Iend,Jstr,Jend,grd_scale)
!
!==================================================================
!
! Compute horizontal diffusivity coefficients
!
! Patrick Marchesiello and Pierrick Penven, IRD 2007
!
!==================================================================
!
implicit none
integer Istr,Iend,Jstr,Jend, i,j