Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

paresp.F90 @ 1885

Last change on this file since 1885 was 1885, checked in by rblod, 14 years ago
add TAM sources
Property svn:executable set to ``*
File size: 9.6 KB

Line
1	MODULE paresp
2	!!======================================================================
3	!! * MODULE paresp *
4	!! NEMOVAR : Weights for an energy-type scalar product
5	!!======================================================================
6
7	!!----------------------------------------------------------------------
8	!! par_esp : Compute and store weights for an energy-type scalar product
9	!!----------------------------------------------------------------------
10	!! * Modules used
11	USE par_kind, ONLY : & ! Precision variables
12	& wp
13	USE par_oce, ONLY : & ! Ocean space and time domain variables
14	& jpk
15	USE eosbn2, ONLY : & ! Equation of state
16	& ralpha, &
17	& rbeta
18	USE phycst, ONLY : & ! Physical constants
19	& grav, &
20	& rcp
21	USE zdf_oce, ONLY : & ! Vertical mixing parameters
22	& avt0, &
23	& avm0
24	USE dom_oce, ONLY : & ! Ocean space domain parameters
25	& gdept_0, &
26	& e3w_0
27	USE in_out_manager ! I/O stuff
28
29	IMPLICIT NONE
30
31	!! * Routine accessibility
32	PRIVATE
33
34	PUBLIC &
35	& par_esp !: Compute and store energy weights
36
37	!! * Share module variables
38	REAL(wp), DIMENSION(jpk), PUBLIC :: &
39	& wesp_t, & !: Normalized energy weights for temperature
40	& wesp_s !: Normalized energy weights for salinity
41	REAL(wp), PUBLIC :: &
42	& wesp_u, & !: Normalized energy weights for velocity
43	& wesp_ssh, & !: Normalized energy weights for SSH
44	& wesp_tau, & !: Normalized energy weights for windstress
45	& wesp_q, & !: Normalized energy weights for heat flux
46	& wesp_emp !: Normalized energy weights for EmP
47
48	CONTAINS
49
50	SUBROUTINE par_esp
51	!!-----------------------------------------------------------------------
52	!!
53	!! * ROUTINE par_esp *
54	!!
55	!! ** Purpose : Compute and store weights for an energy-type scalar
56	!! product.
57	!!
58	!! ** Method :
59	!! The weighting coefficients are computed from an analytical
60	!! linear density profile which is similar to the one used in
61	!! default option in OPA direct.
62	!!
63	!! The weighting coefficients for the forcing fields are estimated
64	!! from the surface boundary conditions.
65	!!
66	!! The weighting coefficients for the velocity variables are
67	!! normalized to one.
68	!!
69	!! wesp_t(k) = ( g / N(k) )^2 x alpha^2
70	!! wesp_s(k) = ( g / N(k) )^2 x beta^2
71	!! wesp_u = 1
72	!! wesp_ssh = g
73	!! wesp_tau = ( ( e3w(1) / ( avu x rho0 ) )^2
74	!! wesp_q = ( ( g / N(1) x alpha
75	!! x e3w(1) / ( avT x rho0 x cp ) )^2
76	!! wesp_emp = ( ( g / N(1) x beta
77	!! x e3w(1) x S(1) / avT )^2
78	!!
79	!! where
80	!!
81	!! N(k) is the Brunt-Vaisala frequency (depth-dependent)
82	!! alpha is the thermal expansion coefficient
83	!! beta is the salinity expansion coefficient
84	!! g is gravity
85	!! avu is vertical eddy viscosity coefficent for dynamics
86	!! avT is vertical eddy diffusivity coefficent for tracers
87	!! e3w(1) is the vertical scale factor at top w-point
88	!! rho0 is a reference density
89	!! S(1) is a reference salinity at the surface
90	!!
91	!! The complete scalar product < ; > for the state vector
92	!!
93	!! dx = ( du , dv , dT , dS , deta , dq , demp , dtaux , dtauy )
94	!!
95	!! involves the volume/area elements:
96	!!
97	!! < dx ; dx > = sum_i,j ( deta^2 * e1T * e2T * wesp_ssh
98	!! + dq^2 * e1T * e2T * e3w(1) * wesp_q
99	!! + demp^2 * e1T * e2T * e3w(1) * wesp_emp
100	!! + dtaux^2 * e1u * e2u * e3uw(1) * wesp_tau
101	!! + dtauy^2 * e1v * e2v * e3vw(1) * wesp_tau
102	!! + sum_k ( du^2 * e1u * e2u * e3u * wesp_u
103	!! + dv^2 * e1v * e2v * e3v * wesp_u
104	!! + dT^2 * e1T * e2T * e3T * wesp_T
105	!! + dS^2 * e1T * e2T * e3T * wesp_S ) )
106	!!
107	!! ** Action :
108	!!
109	!! History :
110	!! ! 97-06 (A. Weaver, J. Vialard) OPAVAR version
111	!! ! 07-11 (A. Weaver) NEMOVAR version based on OPAVAR paresp.F
112	!!-----------------------------------------------------------------------
113	!! * Modules used
114
115	!! * Arguments
116
117	!! * Local declarations
118	REAL(KIND=wp) :: &
119	& zgrau0, &
120	& zk, &
121	& zt0, &
122	& zs0, &
123	& zrau0
124	REAL(KIND=wp), DIMENSION(jpk) :: &
125	& ztan, &
126	& zsan, &
127	& zrhan, &
128	& zbn2an
129	INTEGER :: &
130	& jk
131
132	!--------------------------------------------------------------------
133	! Local constant initialization
134	!--------------------------------------------------------------------
135
136	zt0 = 19.0_wp ! Reference temperature
137	zs0 = 35.0_wp ! Reference salinity
138	zrau0 = 1025.0_wp ! Reference density
139
140	zgrau0 = -1.0_wp * grav / zrau0
141
142	!--------------------------------------------------------------------
143	! Initialize rho with an analytical profile
144	!--------------------------------------------------------------------
145
146	! Define analytical temperature profile
147
148	DO jk = 1, jpk
149	ztan(jk) = 7.5_wp * ( 1.0_wp - TANH( ( gdept_0(jk) - 80.0_wp ) &
150	& / 30.0_wp ) ) &
151	& + 10.0_wp * ( 5300.0_wp - gdept_0(jk) ) / 5300.0_wp
152	ztan(jk) = ABS( ztan(jk) )
153	END DO
154	#if defined key_pomme_r025
155	ztan(44) = 0.20
156	ztan(45) = 0.15
157	ztan(46) = 0.10
158	#endif
159
160	! Define analytical salinity profile
161
162	DO jk = 1, jpk
163	zsan(jk) = 35.50_wp
164	END DO
165
166	! Define analytical density profile
167
168	DO jk = 1, jpk
169	zrhan(jk) = zrau0 * ( 1.0_wp - ralpha * ( ztan(jk) - zt0 ) &
170	& + rbeta * ( zsan(jk) - zs0 ) )
171	END DO
172
173	!--------------------------------------------------------------------
174	! Calculate N^2 at T and S points
175	!--------------------------------------------------------------------
176
177	DO jk = 2, jpkm1
178	zbn2an(jk) = ( zgrau0 / 2.0_wp ) &
179	& * ( ( ( zrhan(jk-1) - zrhan(jk ) ) / e3w_0(jk ) ) &
180	& + ( ( zrhan(jk ) - zrhan(jk+1) ) / e3w_0(jk+1) ) )
181	END DO
182
183	zbn2an(1) = zbn2an(2)
184	zbn2an(jpk) = zbn2an(jpkm1)
185
186	!--------------------------------------------------------------------
187	! Calculate energy-type weights
188	!--------------------------------------------------------------------
189
190	DO jk = 1, jpk
191	zk = ( grav * grav ) / zbn2an(jk)
192	wesp_t(jk) = zk * ralpha * ralpha
193	wesp_s(jk) = zk * rbeta * rbeta
194	END DO
195
196	wesp_u = 1.0_wp
197	wesp_ssh = grav
198	wesp_tau = ( e3w_0(1) / ( avm0 * zrau0 ) )**2
199	wesp_q = wesp_t(1) * ( e3w_0(1) / ( avt0 * zrau0 * rcp ) )**2
200	wesp_emp = wesp_s(1) * ( e3w_0(1) * zsan(1) / avt0 )**2
201
202	!--------------------------------------------------------------------
203	! Print
204	!--------------------------------------------------------------------
205
206	IF (lwp) THEN
207	WRITE(numout,*)
208	WRITE(numout,*) ' par_esp: Analytical T, S, rho, N^2 profiles'
209	WRITE(numout,*) ' -------'
210	WRITE(numout,*)
211	WRITE(numout,995)
212	WRITE(numout,996)
213	DO jk = 1, jpk
214	IF(lwp)WRITE(numout,1007) jk, ztan(jk), zsan(jk), &
215	& zrhan(jk), zbn2an(jk)
216	END DO
217	WRITE(numout,*)
218	WRITE(numout,*) ' par_esp: Energy scalar product weights', &
219	& ' for T and S'
220	WRITE(numout,*) ' -------'
221	WRITE(numout,*)
222	WRITE(numout,998)
223	WRITE(numout,999)
224	DO jk = 1, jpk
225	WRITE(numout,1006) jk, wesp_t(jk), wesp_s(jk)
226	END DO
227	WRITE(numout,*)
228	WRITE(numout,*) ' par_esp: Energy scalar product weights', &
229	& ' for u, SSH, tau, Q and EmP'
230	WRITE(numout,*) ' -------'
231	WRITE(numout,*)
232	WRITE(numout,*) ' wesp_u = ', wesp_u
233	WRITE(numout,*) ' wesp_ssh = ', wesp_ssh
234	WRITE(numout,*) ' wesp_tau = ', wesp_tau
235	WRITE(numout,*) ' wesp_q = ', wesp_q
236	WRITE(numout,*) ' wesp_emp = ', wesp_emp
237
238	995 FORMAT(' level T S ', &
239	& ' rho N^2 ')
240	996 FORMAT(' ----- -------- --------', &
241	& ' ---------- ----------')
242	998 FORMAT(' level wesp_T wesp_S ')
243	999 FORMAT(' ----- ---------- ------------')
244	1005 FORMAT(4X,I2,3X,E12.5,1X,E12.5)
245	1006 FORMAT(4X,I2,3X,E12.5,1X,E12.5)
246	1007 FORMAT(4X,I2,6X,F10.5,5X,F10.5,5X,F10.5,5X,E10.5)
247	CALL flush( numout )
248
249	ENDIF
250
251	DO jk = 1, jpk
252	IF ( zbn2an(jk) < 0.0_wp ) THEN
253	IF(lwp)WRITE(numout,990) zbn2an(jk), jk
254	990 FORMAT(' paresp : Error unstable density profile;', &
255	& ' zbn2an = ',E10.5,' at level ',I3)
256	CALL ctl_stop( 'paresp; unstable density profile' )
257	ENDIF
258	ENDDO
259
260	END SUBROUTINE par_esp
261
262	END MODULE paresp

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/paresp.F90 @ 1885

Download in other formats: