MODULE ldfdyn_smag !!====================================================================== !! *** MODULE ldftrasmag *** !! Ocean physics: variable eddy induced velocity coefficients !!====================================================================== #if defined key_dynldf_smag && defined key_dynldf_c3d !!---------------------------------------------------------------------- !! 'key_dynldf_smag' and smagorinsky diffusivity !! 'key_dynldf_c3d' 3D tracer lateral mixing coef. !!---------------------------------------------------------------------- !! ldf_eiv : compute the eddy induced velocity coefficients !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE sbc_oce ! surface boundary condition: ocean USE sbcrnf ! river runoffs USE ldfdyn_oce ! ocean tracer lateral physics USE phycst ! physical constants USE ldfslp ! iso-neutral slopes USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE iom IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC ldf_dyn_smag ! routine called by step.F90 !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z $ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- CONTAINS !!---------------------------------------------------------------------- !! *** ldfdyn_smag.F90 *** !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id: ldfdyn_c3d.h90 1581 2009-08-05 14:53:12Z smasson $ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 'key_dynldf_smag' 3D lateral eddy viscosity coefficients !!---------------------------------------------------------------------- SUBROUTINE ldf_dyn_smag( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE ldf_dyn_smag *** !! !! ** Purpose : initializations of the horizontal ocean physics !! !! ** Method : 3D eddy viscosity coef. !! M.Griffies, R.Hallberg AMS, 2000 !! for laplacian: !! Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=3-4 !! for bilaplacian: !! Bsmag=Asmag*dx*dy/8 !! D^2=(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates !! in general case du/dx ==> e2 d(u/e2)/dx; du/dy ==> e1 d(u/e1)/dy; !! dv/dx ==> e2 d(v/e2)/dx; dv/dy ==> e1 d(v/e1)/dy !! !! laplacian operator : ahm1, ahm2 defined at T- and F-points !! ahm3, ahm4 never used !! bilaplacian operator : ahm1, ahm2 never used !! : ahm3, ahm4 defined at U- and V-points !! explanation of the default is missingi !! last modified : Maria Luneva, September 2011 !!---------------------------------------------------------------------- !! * Modules used !! ahm0 here is a background viscosity !! * Arguments USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released USE wrk_nemo, ONLY: zux => wrk_2d_1 , zuy => wrk_2d_2 , zvx => wrk_2d_3 , zvy => wrk_2d_4 ! 2D workspace USE wrk_nemo, ONLY: zue1 => wrk_2d_5 , zue2 => wrk_2d_6 , zve1 => wrk_2d_7 , zve2 => wrk_2d_8 ! 2D workspac !! * local variables INTEGER :: kt ! timestep INTEGER :: ji, jj, jk ! dummy loop indices REAL (wp):: rdeltat,rdeltaf,rdeltau,rdeltav ! temporary scalars !!---------------------------------------------------------------------- IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN CALL ctl_stop('dyn_ldf_smag : requested workspace arrays unavailable') ; RETURN ENDIF IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ldf_dyn_smag : 3D lateral eddy viscosity coefficient' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF ! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operators ! ================= whatever its orientation is) IF( ln_dynldf_lap ) THEN ! define ahm1 and ahm2 at the right grid point position ! (USER: modify ahm1 and ahm2 following your desiderata) DO jk=1,jpk zue2(:,:)=un(:,:,jk)/e2u(:,:) zve1(:,:)=vn(:,:,jk)/e1v(:,:) zue1(:,:)=un(:,:,jk)/e1u(:,:) zve2(:,:)=vn(:,:,jk)/e2v(:,:) DO jj=2,jpj DO ji=2,jpi zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) ENDDO ENDDO DO jj=1,jpjm1 DO ji=1,jpi zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk) zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk) ENDDO ENDDO DO jj=2,jpjm1 DO ji=2,jpim1 rdeltat=2./(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2)) rdeltaf=2./(e1f(ji,jj)**(-2)+e2f(ji,jj)**(-2)) ahm1(ji,jj,jk)=(rn_cmsmag_1/3.14)**2*rdeltat* & sqrt( (zux(ji,jj)-zvy(ji,jj))**2+ & 0.0625*(zuy(ji,jj)+zuy(ji,jj-1)+zuy(ji-1,jj)+zuy(ji-1,jj-1)+ & zvx(ji,jj)+zvx(ji,jj-1)+zvx(ji-1,jj)+zvx(ji-1,jj-1))**2) ahm2(ji,jj,jk)=(rn_cmsmag_1/3.14)**2*rdeltaf* & sqrt( (zuy(ji,jj)+zvx(ji,jj))**2+ & 0.0625*(zux(ji,jj)+zux(ji,jj+1)+zux(ji+1,jj)+zux(ji+1,jj+1)- & zvy(ji,jj)-zvy(ji,jj+1)-zvy(ji+1,jj)-zvy(ji+1,jj+1))**2) ahm1(ji,jj,jk)=MAX(ahm1(ji,jj,jk),rn_ahm_0_lap) ahm2(ji,jj,jk)=MAX(ahm2(ji,jj,jk),rn_ahm_0_lap) ! stability criteria ahm1(ji,jj,jk)=MIN(ahm1(ji,jj,jk),rdeltat**2/(16*rdt)) ahm2(ji,jj,jk)=MIN(ahm2(ji,jj,jk),rdeltaf**2/(16*rdt)) ENDDO ENDDO ENDDO ! jpk ahm1(:,:,jpk) = ahm1(:,:,jpkm1) ahm2(:,:,jpk) = ahm2(:,:,jpkm1) IF(lwp.and.kt==nit000) WRITE(numout,'(36x," ahm ", 7x)') DO jk = 1, jpk IF(lwp.and.kt==nit000) WRITE(numout,'(30x,E10.2,8x,i3)') ahm1(jpi/2,jpj/2,jk), jk END DO CALL lbc_lnk( ahm1, 'T', 1. ) ! Lateral boundary conditions on ( ahtt ) CALL lbc_lnk( ahm2, 'F', 1. ) ! Lateral boundary conditions on ( ahtt ) ENDIF ! ln_dynldf ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator ! ================================ whatever its orientation is) ! (USER: modify ahm3 and ahm4 following your desiderata) ! Here: ahm is proportional to the cube of the maximum of the gridspacing ! in the to horizontal direction IF( ln_dynldf_bilap ) THEN DO jk=1,jpk zue2(:,:)=un(:,:,jk)/e2u(:,:) zve1(:,:)=vn(:,:,jk)/e1v(:,:) zue1(:,:)=un(:,:,jk)/e1u(:,:) zve2(:,:)=vn(:,:,jk)/e2v(:,:) DO jj=2,jpj DO ji=2,jpi zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) ENDDO ENDDO DO jj=1,jpjm1 DO ji=1,jpim1 zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk) zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk) ENDDO ENDDO DO jj=2,jpjm1 DO ji=2,jpim1 rdeltau=2./(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2)) rdeltav=2./(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2)) ahm3(ji,jj,jk)=MIN(rn_ahm_0_blp,(rn_cmsmag_2/3.14)**2/8*rdeltau**2* & sqrt(0.25*(zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj))**2+ & 0.25*(zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1))**2)) ahm4(ji,jj,jk)=MIN(rn_ahm_0_blp ,(rn_cmsmag_2/3.14)**2/8*rdeltav**2* & sqrt(0.25*(zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1))**2+ & 0.25*(zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj))**2)) ! stability criteria ahm3(ji,jj,jk)=MAX(ahm3(ji,jj,jk),-rdeltau**2/(16*rdt)) ahm4(ji,jj,jk)=MAX(ahm4(ji,jj,jk),-rdeltav**2/(16*rdt)) ENDDO ENDDO ENDDO ahm3(:,:,jpk) = ahm3(:,:,jpkm1) ahm4(:,:,jpk) = ahm4(:,:,jpkm1) DO jk = 1, jpk IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(jpi/2,jpj/2,jk), jk ENDIF END DO CALL lbc_lnk( ahm3, 'U', 1. ) ! Lateral boundary conditions CALL lbc_lnk( ahm4, 'V', 1. ) ENDIF IF( wrk_not_released(2, 1,2,3,4,5,6,7,8)) CALL ctl_stop('dyn_ldf_smag: failed to release workspace arrays') END SUBROUTINE ldf_dyn_smag #else !!---------------------------------------------------------------------- !! Default option Dummy module !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ldf_dyn_smag( kt ) ! Empty routine WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt END SUBROUTINE ldf_dyn_smag #endif END MODULE ldfdyn_smag