1 | MODULE ldftra_smag |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE ldftrasmag *** |
---|
4 | !! Ocean physics: variable eddy induced velocity coefficients |
---|
5 | !!====================================================================== |
---|
6 | !! last modified : Maria Luneva, October 2012 |
---|
7 | !!---------------------------------------------------------------------- |
---|
8 | #if defined key_traldf_smag |
---|
9 | !!---------------------------------------------------------------------- |
---|
10 | !! 'key_traldf_smag' and smagorinsky diffusivity |
---|
11 | !!---------------------------------------------------------------------- |
---|
12 | !! ldf_tra_smag : compute the smagorinski eddy coefficients |
---|
13 | !!---------------------------------------------------------------------- |
---|
14 | USE oce ! ocean dynamics and tracers |
---|
15 | USE dom_oce ! ocean space and time domain |
---|
16 | USE sbc_oce ! surface boundary condition: ocean |
---|
17 | USE sbcrnf ! river runoffs |
---|
18 | USE ldftra ! ocean tracer lateral physics |
---|
19 | USE phycst ! physical constants |
---|
20 | USE ldfslp ! iso-neutral slopes |
---|
21 | ! |
---|
22 | USE in_out_manager ! I/O manager |
---|
23 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
24 | USE prtctl ! Print control |
---|
25 | USE iom ! |
---|
26 | USE ioipsl ! |
---|
27 | USE wrk_nemo ! |
---|
28 | |
---|
29 | IMPLICIT NONE |
---|
30 | PRIVATE |
---|
31 | |
---|
32 | PUBLIC ldf_tra_smag ! routine called by step.F90 |
---|
33 | |
---|
34 | !! * Substitutions |
---|
35 | # include "domzgr_substitute.h90" |
---|
36 | # include "vectopt_loop_substitute.h90" |
---|
37 | !!---------------------------------------------------------------------- |
---|
38 | !! NEMO/OPA 3.6 , LOCEAN-IPSL (2014) |
---|
39 | !! $Id: ldf_tra_smag.F90 1482 2010-06-13 15:28:06Z $ |
---|
40 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
41 | !!---------------------------------------------------------------------- |
---|
42 | CONTAINS |
---|
43 | |
---|
44 | SUBROUTINE ldf_tra_smag( kt ) |
---|
45 | !!---------------------------------------------------------------------- |
---|
46 | !! *** ROUTINE ldf_tra_smag *** |
---|
47 | !! |
---|
48 | !! ** Purpose : initializations of the horizontal ocean physics |
---|
49 | !! |
---|
50 | !! ** Method : 3D eddy viscosity coef. |
---|
51 | !! M.Griffies, R.Hallberg AMS, 2000 |
---|
52 | !! for laplacian: |
---|
53 | !! Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=1 for tracers, C=3-4 for viscosity |
---|
54 | !! D^2= rm_smsh*(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates |
---|
55 | !! IF rm_smsh = 0 , only shear is used, recommended for tidal flows |
---|
56 | !! in general case du/dx ==> e2 d(u/e2)/dx; du/dy ==> e1 d(u/e1)/dy; |
---|
57 | !! dv/dx ==> e2 d(v/e2)/dx; dv/dy ==> e1 d(v/e1)/dy |
---|
58 | !! for bilaplacian: now this option is deleted as unstable or non-conservative |
---|
59 | !! - delta{Bsmag (delta(T)} = -Bsmag* delta{delta(T)} - delta(Bsmag)*delta( T ) |
---|
60 | !! second term is of arbitrary sign on the edge of fronts and can induce instability |
---|
61 | !! Bsmag=Asmag*dx*dy/8 |
---|
62 | !! |
---|
63 | !! laplacian operator : ahm1, ahm2 defined at T- and F-points |
---|
64 | !! ahm3, ahm4 never used |
---|
65 | !! bilaplacian operator : ahm1, ahm2 never used |
---|
66 | !! : ahm3, ahm4 defined at U- and V-points |
---|
67 | !! ??? explanation of the default is missing |
---|
68 | !!---------------------------------------------------------------------- |
---|
69 | INTEGER, INTENT( in ) :: kt ! ocean time-step inedx |
---|
70 | ! |
---|
71 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
72 | REAL(wp) :: zdeltau, zhsmag ! local scalars |
---|
73 | REAL(wp) :: zdeltav, zsmsh , zcoef ! - - |
---|
74 | REAL(wp), POINTER , DIMENSION (:,:) :: zux, zvx , zuy , zvy |
---|
75 | REAL(wp), POINTER , DIMENSION (:,:) :: zue1, zue2 , zve1 , zve2 |
---|
76 | |
---|
77 | CALL wrk_alloc (jpi,jpj,zux, zvx , zuy , zvy ) |
---|
78 | CALL wrk_alloc (jpi,jpj,zue1, zue2 , zve1 , zve2 ) |
---|
79 | !!---------------------------------------------------------------------- |
---|
80 | ! |
---|
81 | IF( kt == nit000 ) THEN |
---|
82 | IF(lwp) WRITE(numout,*) |
---|
83 | IF(lwp) WRITE(numout,*) ' ldf_tra_smag : 3D eddy smagorinsky diffusivity ' |
---|
84 | IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~ -- ' |
---|
85 | ENDIF |
---|
86 | |
---|
87 | zhsmag = rn_chsmag |
---|
88 | zsmsh = rn_smsh |
---|
89 | zux(:,:) = 0._wp ; zuy(:,:) = 0._wp ; zvx(:,:) = 0._wp ; zvy(:,:) = 0._wp |
---|
90 | |
---|
91 | ! ------------------- |
---|
92 | ahtt(:,:,:) = rn_aht_0 |
---|
93 | IF( ln_traldf_bilap ) THEN |
---|
94 | IF( lwp .AND. kt == nit000) WRITE(numout,* ) 'ldf_tra_smag :no bilaplacian Smagorinsky diffusivity' |
---|
95 | IF( lwp .AND. kt == nit000) WRITE(numout,* ) 'ldf_tra_smag :bilaplacian diffusivity set to constant' |
---|
96 | ENDIF |
---|
97 | |
---|
98 | |
---|
99 | |
---|
100 | ! harmonic operator (U-, V-, W-points) |
---|
101 | ! ----------------- |
---|
102 | ahtu(:,:,:) = rn_aht_0 ! set ahtu , ahtv at u- and v-points, |
---|
103 | ahtv(:,:,:) = rn_aht_0 ! and ahtw at w-point |
---|
104 | |
---|
105 | IF( ln_traldf_lap ) THEN |
---|
106 | ! |
---|
107 | DO jk = 1 , jpkm1 |
---|
108 | zue2(:,:) = un(:,:,jk) / e2u(:,:) !!gm for stability reason use of before instead of now here !!!! |
---|
109 | zve1(:,:) = vn(:,:,jk) / e1v(:,:) |
---|
110 | zue1(:,:) = un(:,:,jk) / e1u(:,:) |
---|
111 | zve2(:,:) = vn(:,:,jk) / e2v(:,:) |
---|
112 | ! |
---|
113 | DO jj = 2, jpj !!gm multiplication by tmask useless as un, vn maked field ! |
---|
114 | DO ji= 2, jpi |
---|
115 | zux(ji,jj) = ( zue2(ji,jj) - zue2(ji-1,jj ) ) / e1e2t(ji,jj) * tmask(ji,jj,jk) * zsmsh |
---|
116 | zvy(ji,jj) = ( zve1(ji,jj) - zve1(ji ,jj-1) ) / e1e2t(ji,jj) * tmask(ji,jj,jk) * zsmsh |
---|
117 | END DO |
---|
118 | END DO |
---|
119 | ! |
---|
120 | DO jj = 1, jpjm1 |
---|
121 | DO ji = 1, jpim1 |
---|
122 | zuy(ji,jj) = ( zue1(ji ,jj+1) - zue1(ji,jj) ) / e2f(ji,jj) *e1f(ji,jj) * fmask(ji,jj,jk) |
---|
123 | zvx(ji,jj) = ( zve2(ji+1,jj ) - zve2(ji,jj) ) / e1f(ji,jj) *e2f(ji,jj) * fmask(ji,jj,jk) |
---|
124 | END DO |
---|
125 | END DO |
---|
126 | ! |
---|
127 | DO jj = 2, jpjm1 |
---|
128 | DO ji = 2, jpim1 |
---|
129 | zdeltau = 2._wp / ( e1u(ji,jj)**(-2) + e2u(ji,jj)**(-2) ) |
---|
130 | zdeltav = 2._wp / ( e1v(ji,jj)**(-2) + e2v(ji,jj)**(-2) ) |
---|
131 | ! |
---|
132 | ahtu(ji,jj,jk) = MAX( rn_aht_0 , (zhsmag/rpi)**2*zdeltau* & |
---|
133 | & SQRT( 0.25_wp*( zux(ji,jj)+zux(ji+1,jj )-zvy(ji,jj)-zvy(ji+1,jj ) )**2 & |
---|
134 | & + 0.25_wp*( zuy(ji,jj)+zuy(ji ,jj-1)+zvx(ji,jj)+zvx(ji ,jj-1) )**2 ) ) |
---|
135 | ! |
---|
136 | ahtv(ji,jj,jk) = MAX( rn_aht_0 , (zhsmag/rpi)**2*zdeltav* & |
---|
137 | & SQRT( 0.25_wp*( zux(ji,jj)+zux(ji ,jj+1)-zvy(ji ,jj)-zvy(ji,jj+1) )**2 & |
---|
138 | & + 0.25_wp*( zuy(ji,jj)+zuy(ji-1,jj )+zvx(ji-1,jj)+zvx(ji,jj ) )**2 ) ) |
---|
139 | ! |
---|
140 | ! stability criteria: aht<delta**2/(4*dt) dt=2*rdt , positiveness require aht<delta**2/(8*dt) |
---|
141 | ahtu(ji,jj,jk) = MIN( ahtu(ji,jj,jk) , zdeltau / (16*rdt) , rn_aht_m ) |
---|
142 | ahtv(ji,jj,jk) = MIN( ahtv(ji,jj,jk) , zdeltav / (16*rdt) , rn_aht_m ) |
---|
143 | END DO |
---|
144 | END DO |
---|
145 | END DO |
---|
146 | ENDIF |
---|
147 | ahtu(:,:,jpk) = ahtu(:,:,jpkm1) |
---|
148 | ahtv(:,:,jpk) = ahtv(:,:,jpkm1) |
---|
149 | CALL lbc_lnk( ahtu, 'U', 1. ) ! Lateral boundary conditions |
---|
150 | CALL lbc_lnk( ahtv, 'V', 1. ) |
---|
151 | ! |
---|
152 | CALL wrk_dealloc ( jpi,jpj,zux, zvx , zuy , zvy ) |
---|
153 | CALL wrk_dealloc ( jpi,jpj,zue1, zue2 , zve1 , zve2 ) |
---|
154 | ! |
---|
155 | END SUBROUTINE ldf_tra_smag |
---|
156 | |
---|
157 | #else |
---|
158 | !!---------------------------------------------------------------------- |
---|
159 | !! Default option Dummy module |
---|
160 | !!---------------------------------------------------------------------- |
---|
161 | CONTAINS |
---|
162 | SUBROUTINE ldf_tra_smag( kt ) ! Empty routine |
---|
163 | WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt |
---|
164 | END SUBROUTINE ldf_tra_smag |
---|
165 | #endif |
---|
166 | |
---|
167 | !!====================================================================== |
---|
168 | END MODULE ldftra_smag |
---|