1 | MODULE dynldf_lap_blp |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE dynldf_lap_blp *** |
---|
4 | !! Ocean dynamics: lateral viscosity trend (laplacian and bilaplacian) |
---|
5 | !!====================================================================== |
---|
6 | !! History : OPA ! 1990-09 (G. Madec) Original code |
---|
7 | !! 4.0 ! 1991-11 (G. Madec) |
---|
8 | !! 6.0 ! 1996-01 (G. Madec) statement function for e3 and ahm |
---|
9 | !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module |
---|
10 | !! - ! 2004-08 (C. Talandier) New trends organization |
---|
11 | !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, |
---|
12 | !! ! add velocity dependent coefficient and optional read in file |
---|
13 | !!---------------------------------------------------------------------- |
---|
14 | |
---|
15 | !!---------------------------------------------------------------------- |
---|
16 | !! dyn_ldf_lap : update the momentum trend with the lateral viscosity using an iso-level laplacian operator |
---|
17 | !! dyn_ldf_blp : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator |
---|
18 | !!---------------------------------------------------------------------- |
---|
19 | USE oce ! ocean dynamics and tracers |
---|
20 | USE dom_oce ! ocean space and time domain |
---|
21 | USE ldfdyn ! lateral diffusion: eddy viscosity coef. |
---|
22 | USE ldfslp ! iso-neutral slopes |
---|
23 | USE zdf_oce ! ocean vertical physics |
---|
24 | ! |
---|
25 | USE in_out_manager ! I/O manager |
---|
26 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
27 | USE wrk_nemo ! Memory Allocation |
---|
28 | USE timing ! Timing |
---|
29 | |
---|
30 | IMPLICIT NONE |
---|
31 | PRIVATE |
---|
32 | |
---|
33 | PUBLIC dyn_ldf_lap ! called by dynldf.F90 |
---|
34 | PUBLIC dyn_ldf_blp ! called by dynldf.F90 |
---|
35 | |
---|
36 | !! * Substitutions |
---|
37 | # include "domzgr_substitute.h90" |
---|
38 | # include "vectopt_loop_substitute.h90" |
---|
39 | !!---------------------------------------------------------------------- |
---|
40 | !! NEMO/OPA 3.7 , NEMO Consortium (2014) |
---|
41 | !! $Id$ |
---|
42 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
43 | !!---------------------------------------------------------------------- |
---|
44 | CONTAINS |
---|
45 | |
---|
46 | SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) |
---|
47 | !!---------------------------------------------------------------------- |
---|
48 | !! *** ROUTINE dyn_ldf_lap *** |
---|
49 | !! |
---|
50 | !! ** Purpose : Compute the before horizontal momentum diffusive |
---|
51 | !! trend and add it to the general trend of momentum equation. |
---|
52 | !! |
---|
53 | !! ** Method : The Laplacian operator apply on horizontal velocity is |
---|
54 | !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) |
---|
55 | !! |
---|
56 | !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. |
---|
57 | !!---------------------------------------------------------------------- |
---|
58 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
59 | INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage |
---|
60 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity [m/s] |
---|
61 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! velocity trend [m/s2] |
---|
62 | ! |
---|
63 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
64 | REAL(wp) :: zsign ! local scalars |
---|
65 | REAL(wp) :: zua, zva ! local scalars |
---|
66 | REAL(wp), POINTER, DIMENSION(:,:) :: zcur, zdiv |
---|
67 | !!---------------------------------------------------------------------- |
---|
68 | ! |
---|
69 | IF( kt == nit000 .AND. lwp ) THEN |
---|
70 | WRITE(numout,*) |
---|
71 | WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass |
---|
72 | WRITE(numout,*) '~~~~~~~ ' |
---|
73 | ENDIF |
---|
74 | ! |
---|
75 | IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_lap') |
---|
76 | ! |
---|
77 | CALL wrk_alloc( jpi, jpj, zcur, zdiv ) |
---|
78 | ! |
---|
79 | IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign |
---|
80 | ELSE ; zsign = -1._wp ! (eddy viscosity coef. >0) |
---|
81 | ENDIF |
---|
82 | ! |
---|
83 | ! ! =============== |
---|
84 | DO jk = 1, jpkm1 ! Horizontal slab |
---|
85 | ! ! =============== |
---|
86 | DO jj = 2, jpj |
---|
87 | DO ji = fs_2, jpi ! vector opt. |
---|
88 | ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) |
---|
89 | !!gm open question here : fse3f at before or now ? probably now... |
---|
90 | !!gm note that ahmf has already been multiplied by fmask |
---|
91 | zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * fse3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & |
---|
92 | & * ( e2v(ji ,jj-1) * pvb(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) & |
---|
93 | & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) |
---|
94 | ! ! ahm * div (computed from 2 to jpi/jpj) |
---|
95 | zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / fse3t(ji,jj,jk) * tmask(ji,jj,jk) & |
---|
96 | & * ( e2u(ji,jj)*fse3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*fse3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & |
---|
97 | & + e1v(ji,jj)*fse3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*fse3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) |
---|
98 | END DO |
---|
99 | END DO |
---|
100 | ! |
---|
101 | DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) |
---|
102 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
103 | pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( & |
---|
104 | & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) & |
---|
105 | & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) |
---|
106 | ! |
---|
107 | pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( & |
---|
108 | & ( zcur(ji,jj ) - zcur(ji-1,jj) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) & |
---|
109 | & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) |
---|
110 | END DO |
---|
111 | END DO |
---|
112 | ! ! =============== |
---|
113 | END DO ! End of slab |
---|
114 | ! ! =============== |
---|
115 | CALL wrk_dealloc( jpi, jpj, zcur, zdiv ) |
---|
116 | ! |
---|
117 | IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_lap') |
---|
118 | ! |
---|
119 | END SUBROUTINE dyn_ldf_lap |
---|
120 | |
---|
121 | |
---|
122 | SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) |
---|
123 | !!---------------------------------------------------------------------- |
---|
124 | !! *** ROUTINE dyn_ldf_blp *** |
---|
125 | !! |
---|
126 | !! ** Purpose : Compute the before lateral momentum viscous trend |
---|
127 | !! and add it to the general trend of momentum equation. |
---|
128 | !! |
---|
129 | !! ** Method : The lateral viscous trends is provided by a bilaplacian |
---|
130 | !! operator applied to before field (forward in time). |
---|
131 | !! It is computed by two successive calls to dyn_ldf_lap routine |
---|
132 | !! |
---|
133 | !! ** Action : pta updated with the before rotated bilaplacian diffusion |
---|
134 | !!---------------------------------------------------------------------- |
---|
135 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
136 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity fields |
---|
137 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! momentum trend |
---|
138 | ! |
---|
139 | REAL(wp), POINTER, DIMENSION(:,:,:) :: zulap, zvlap ! laplacian at u- and v-point |
---|
140 | !!---------------------------------------------------------------------- |
---|
141 | ! |
---|
142 | IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_blp') |
---|
143 | ! |
---|
144 | CALL wrk_alloc( jpi, jpj, jpk, zulap, zvlap ) |
---|
145 | ! |
---|
146 | IF( kt == nit000 ) THEN |
---|
147 | IF(lwp) WRITE(numout,*) |
---|
148 | IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' |
---|
149 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' |
---|
150 | ENDIF |
---|
151 | ! |
---|
152 | zulap(:,:,:) = 0._wp |
---|
153 | zvlap(:,:,:) = 0._wp |
---|
154 | ! |
---|
155 | CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 ) ! rotated laplacian applied to ptb (output in zlap) |
---|
156 | ! |
---|
157 | CALL lbc_lnk( zulap(:,:,:) , 'U', -1. ) ! Lateral boundary conditions |
---|
158 | CALL lbc_lnk( zvlap(:,:,:) , 'V', -1. ) |
---|
159 | ! |
---|
160 | CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 ) ! rotated laplacian applied to zlap (output in pta) |
---|
161 | ! |
---|
162 | CALL wrk_dealloc( jpi, jpj, jpk, zulap, zvlap ) |
---|
163 | ! |
---|
164 | IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_blp') |
---|
165 | ! |
---|
166 | END SUBROUTINE dyn_ldf_blp |
---|
167 | |
---|
168 | !!====================================================================== |
---|
169 | END MODULE dynldf_lap_blp |
---|