1 | MODULE ldfc1d_c2d |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE ldfc1d_c2d *** |
---|
4 | !! Ocean physics: profile and horizontal shape of lateral eddy coefficients |
---|
5 | !!===================================================================== |
---|
6 | !! History : 3.7 ! 2013-12 (G. Madec) restructuration/simplification of aht/aeiv specification, |
---|
7 | !! ! add velocity dependent coefficient and optional read in file |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | |
---|
10 | !!---------------------------------------------------------------------- |
---|
11 | !! ldf_c1d : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m) |
---|
12 | !! ldf_c2d : ah = F(e1,e2) (laplacian or = F(e1^3,e2^3) (bilaplacian) |
---|
13 | !!---------------------------------------------------------------------- |
---|
14 | USE oce ! ocean dynamics and tracers |
---|
15 | USE dom_oce ! ocean space and time domain |
---|
16 | USE phycst ! physical constants |
---|
17 | ! |
---|
18 | USE in_out_manager ! I/O manager |
---|
19 | USE lib_mpp ! distribued memory computing library |
---|
20 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
21 | |
---|
22 | IMPLICIT NONE |
---|
23 | PRIVATE |
---|
24 | |
---|
25 | PUBLIC ldf_c1d ! called by ldftra and ldfdyn modules |
---|
26 | PUBLIC ldf_c2d ! called by ldftra and ldfdyn modules |
---|
27 | |
---|
28 | REAL(wp) :: r1_2 = 0.5_wp ! =1/2 |
---|
29 | REAL(wp) :: r1_4 = 0.25_wp ! =1/4 |
---|
30 | REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 |
---|
31 | |
---|
32 | !!---------------------------------------------------------------------- |
---|
33 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
---|
34 | !! $Id$ |
---|
35 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
36 | !!---------------------------------------------------------------------- |
---|
37 | CONTAINS |
---|
38 | |
---|
39 | SUBROUTINE ldf_c1d( cd_type, pahs1, pahs2, pah1, pah2 ) |
---|
40 | !!---------------------------------------------------------------------- |
---|
41 | !! *** ROUTINE ldf_c1d *** |
---|
42 | !! |
---|
43 | !! ** Purpose : 1D eddy diffusivity/viscosity coefficients |
---|
44 | !! |
---|
45 | !! ** Method : 1D eddy diffusivity coefficients F( depth ) |
---|
46 | !! Reduction by zratio from surface to bottom |
---|
47 | !! hyperbolic tangent profile with inflection point |
---|
48 | !! at zh=500m and a width of zw=200m |
---|
49 | !! |
---|
50 | !! cd_type = TRA pah1, pah2 defined at U- and V-points |
---|
51 | !! DYN pah1, pah2 defined at T- and F-points |
---|
52 | !!---------------------------------------------------------------------- |
---|
53 | CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers |
---|
54 | REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pahs1, pahs2 ! surface value of eddy coefficient [m2/s] |
---|
55 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pah1 , pah2 ! eddy coefficient [m2/s] |
---|
56 | ! |
---|
57 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
58 | REAL(wp) :: zh, zc, zdep1 ! local scalars |
---|
59 | REAL(wp) :: zw , zdep2 ! - - |
---|
60 | REAL(wp) :: zratio ! - - |
---|
61 | !!---------------------------------------------------------------------- |
---|
62 | ! |
---|
63 | IF(lwp) WRITE(numout,*) |
---|
64 | IF(lwp) WRITE(numout,*) ' ldf_c1d : set a given profile to eddy mixing coefficients' |
---|
65 | ! |
---|
66 | ! initialization of the profile |
---|
67 | zratio = 0.25_wp ! surface/bottom ratio |
---|
68 | zh = 500._wp ! depth of the inflection point [m] |
---|
69 | zw = 1._wp / 200._wp ! width^-1 - - - [1/m] |
---|
70 | ! ! associated coefficient [-] |
---|
71 | zc = ( 1._wp - zratio ) / ( 1._wp + TANH( zh * zw) ) |
---|
72 | ! |
---|
73 | ! |
---|
74 | SELECT CASE( cd_type ) ! point of calculation |
---|
75 | ! |
---|
76 | CASE( 'DYN' ) ! T- and F-points |
---|
77 | DO jk = jpkm1, 1, -1 ! pah1 at T-point |
---|
78 | pah1(:,:,jk) = pahs1(:,:) * ( zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) ) ) |
---|
79 | END DO |
---|
80 | DO jk = jpkm1, 1, -1 ! pah2 at F-point (zdep2 is an approximation in zps-coord.) |
---|
81 | DO jj = 1, jpjm1 |
---|
82 | DO ji = 1, jpim1 |
---|
83 | zdep2 = ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) & |
---|
84 | & + gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) * r1_4 |
---|
85 | pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) |
---|
86 | END DO |
---|
87 | END DO |
---|
88 | END DO |
---|
89 | CALL lbc_lnk( 'ldfc1d_c2d', pah2, 'F', 1. ) ! Lateral boundary conditions |
---|
90 | ! |
---|
91 | CASE( 'TRA' ) ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) |
---|
92 | DO jk = jpkm1, 1, -1 |
---|
93 | DO jj = 1, jpjm1 |
---|
94 | DO ji = 1, jpim1 |
---|
95 | zdep1 = ( gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk) ) * 0.5_wp |
---|
96 | zdep2 = ( gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk) ) * 0.5_wp |
---|
97 | pah1(ji,jj,jk) = pahs1(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) |
---|
98 | pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) |
---|
99 | END DO |
---|
100 | END DO |
---|
101 | END DO |
---|
102 | ! Lateral boundary conditions |
---|
103 | CALL lbc_lnk_multi( 'ldfc1d_c2d', pah1, 'U', 1. , pah2, 'V', 1. ) |
---|
104 | ! |
---|
105 | CASE DEFAULT ! error |
---|
106 | CALL ctl_stop( 'ldf_c1d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) |
---|
107 | END SELECT |
---|
108 | ! |
---|
109 | END SUBROUTINE ldf_c1d |
---|
110 | |
---|
111 | |
---|
112 | SUBROUTINE ldf_c2d( cd_type, pUfac, knn, pah1, pah2 ) |
---|
113 | !!---------------------------------------------------------------------- |
---|
114 | !! *** ROUTINE ldf_c2d *** |
---|
115 | !! |
---|
116 | !! ** Purpose : 2D eddy diffusivity/viscosity coefficients |
---|
117 | !! |
---|
118 | !! ** Method : 2D eddy diffusivity coefficients F( e1 , e2 ) |
---|
119 | !! laplacian operator : ah proportional to the scale factor [m2/s] |
---|
120 | !! bilaplacian operator : ah proportional to the (scale factor)^3 [m4/s] |
---|
121 | !! In both cases, pah0 is the maximum value reached by the coefficient |
---|
122 | !! at the Equator in case of e1=ra*rad= ~111km, not over the whole domain. |
---|
123 | !! |
---|
124 | !! cd_type = TRA pah1, pah2 defined at U- and V-points |
---|
125 | !! DYN pah1, pah2 defined at T- and F-points |
---|
126 | !!---------------------------------------------------------------------- |
---|
127 | CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers |
---|
128 | REAL(wp) , INTENT(in ) :: pUfac ! =1/2*Uc LAPlacian BiLaPlacian |
---|
129 | INTEGER , INTENT(in ) :: knn ! characteristic velocity [m/s] |
---|
130 | REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pah1, pah2 ! eddy coefficients [m2/s or m4/s] |
---|
131 | ! |
---|
132 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
133 | INTEGER :: inn ! local integer |
---|
134 | !!---------------------------------------------------------------------- |
---|
135 | ! |
---|
136 | IF(lwp) WRITE(numout,*) |
---|
137 | IF(lwp) WRITE(numout,*) ' ldf_c2d : aht = Ufac * max(e1,e2) with Ufac = ', pUfac, ' m/s' |
---|
138 | ! |
---|
139 | ! |
---|
140 | SELECT CASE( cd_type ) !== surface values ==! (chosen grid point function of DYN or TRA) |
---|
141 | ! |
---|
142 | CASE( 'DYN' ) ! T- and F-points |
---|
143 | DO jj = 1, jpj |
---|
144 | DO ji = 1, jpi |
---|
145 | pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn |
---|
146 | pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn |
---|
147 | END DO |
---|
148 | END DO |
---|
149 | CASE( 'TRA' ) ! U- and V-points |
---|
150 | DO jj = 1, jpj |
---|
151 | DO ji = 1, jpi |
---|
152 | pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn |
---|
153 | pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn |
---|
154 | END DO |
---|
155 | END DO |
---|
156 | CASE DEFAULT ! error |
---|
157 | CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) |
---|
158 | END SELECT |
---|
159 | ! !== deeper values = surface one ==! (except jpk) |
---|
160 | DO jk = 2, jpkm1 |
---|
161 | pah1(:,:,jk) = pah1(:,:,1) |
---|
162 | pah2(:,:,jk) = pah2(:,:,1) |
---|
163 | END DO |
---|
164 | ! |
---|
165 | END SUBROUTINE ldf_c2d |
---|
166 | |
---|
167 | !!====================================================================== |
---|
168 | END MODULE ldfc1d_c2d |
---|