Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sfs_nba.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 43.2 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2
3	!==============================================================================
4	!
5	! © 2009. Lawrence Livermore National Security, LLC. All rights reserved.
6	! This work was produced at the Lawrence Livermore National Laboratory (LLNL) under
7	! contract no. DE-AC52-07NA27344 (Contract 44) between the U.S. Department of Energy (DOE)
8	! and Lawrence Livermore National Security, LLC (LLNS) for the operation of LLNL. Copyright
9	! is reserved to Lawrence Livermore National Security, LLC for purposes of controlled
10	! dissemination, commercialization through formal licensing, or other disposition under
11	! terms of Contract 44; DOE policies, regulations and orders; and U.S. statutes. The rights
12	! of the Federal Government are reserved under Contract 44.
13	!
14	! DISCLAIMER
15	! This work was prepared as an account of work sponsored by an agency of the United States
16	! Government. Neither the United States Government nor Lawrence Livermore National
17	! Security, LLC nor any of their employees, makes any warranty, express or implied, or
18	! assumes any liability or responsibility for the accuracy, completeness, or usefulness of
19	! any information, apparatus, product, or process disclosed, or represents that its use
20	! would not infringe privately-owned rights. Reference herein to any specific commercial
21	! products, process, or service by trade name, trademark, manufacturer or otherwise does
22	! not necessarily constitute or imply its endorsement, recommendation, or favoring by the
23	! United States Government or Lawrence Livermore National Security, LLC. The views and
24	! opinions of authors expressed herein do not necessarily state or reflect those of the
25	! United States Government or Lawrence Livermore National Security, LLC, and shall not be
26	! used for advertising or product endorsement purposes.
27	!
28	! LICENSING REQUIREMENTS
29	! Any use, reproduction, modification, or distribution of this software or documentation
30	! for commercial purposes requires a license from Lawrence Livermore National Security,
31	! LLC. Contact: Lawrence Livermore National Laboratory, Industrial Partnerships Office,
32	! P.O. Box 808, L-795, Livermore, CA 94551
33	!
34	!=============================================================================
35	!
36	! Modification History:
37	!
38	! Implemented 12/2009 by Jeff Mirocha, jmirocha@llnl.gov
39	!
40	!=============================================================================
41
42	MODULE module_sfs_nba
43
44	USE module_configure, ONLY : grid_config_rec_type
45
46	IMPLICIT NONE
47
48	REAL :: c1, c2, c3, ce, cb, cs ! global model parameters
49
50	CONTAINS
51
52	!=============================================================================
53
54	SUBROUTINE calc_mij_constants( )
55
56	!-----------------------------------------------------------------------------
57	!
58	! PURPOSE: Compute constants for Mij calculations
59	!
60	!-----------------------------------------------------------------------------
61
62	IMPLICIT NONE
63
64	REAL :: sk, pi ! local model parameters
65
66	!-----------------------------------------------------------------------------
67
68	sk = 0.5
69	pi = 3.1415927
70	cb = 0.36
71
72	cs = ( ( 8.0( 1.0+cb ) )/( 27.0pi2 ) )0.5
73	c1 = ( ( 960.0*0.5 )cb )/( 7.0( 1.0+cb )sk )
74	c2 = c1
75	ce = ( ( 8.0pi/27.0 )( 1.0/3.0 ) )cs**( 4.0/3.0 )
76	c3 = ( ( 27.0/( 8.0pi ) )( 1.0/3.0 ) )cs**( 2.0/3.0 )
77
78	RETURN
79
80	END SUBROUTINE calc_mij_constants
81
82	!=============================================================================
83
84	SUBROUTINE calc_smnsmn( smnsmn, &
85	s11, s22, s33, &
86	s12, s13, s23, &
87	config_flags, &
88	ids, ide, jds, jde, kds, kde, &
89	ims, ime, jms, jme, kms, kme, &
90	ips, ipe, jps, jpe, kps, kpe, &
91	its, ite, jts, jte, kts, kte )
92
93	!-----------------------------------------------------------------------------
94	!
95	! PURPOSE: Compute SmnSmn = S11^2 + S22^2 + S33^2 + 2(S12^2 + S13^2 +S23^2)
96	!
97	!-----------------------------------------------------------------------------
98
99	IMPLICIT NONE
100
101	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
102	:: smnsmn ! Smn*Smn (s-2)
103
104	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
105	:: s11 & ! 2*deformation element 11 (s-1)
106	, s22 & ! 2*deformation element 22 (s-1)
107	, s33 & ! 2*deformation element 33 (s-1)
108	, s12 & ! 2*deformation element 12 (s-1)
109	, s13 & ! 2*deformation element 13 (s-1)
110	, s23 ! 2*deformation element 23 (s-1)
111
112	TYPE (grid_config_rec_type), INTENT( IN ) &
113	:: config_flags
114
115	INTEGER, INTENT( IN ) &
116	:: ids, ide, jds, jde, kds, kde, &
117	ims, ime, jms, jme, kms, kme, &
118	ips, ipe, jps, jpe, kps, kpe, &
119	its, ite, jts, jte, kts, kte
120
121	! LOCAL VARIABLES ------------------------------------------------------------
122
123	REAL :: tmp
124
125	INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf
126
127	!-----------------------------------------------------------------------------
128	!
129	! Set loop limits,
130	! taken from /dyn_em/module_diffusion_em.F SUBROUTINE smag_km
131	!
132	!-----------------------------------------------------------------------------
133
134	ktf = min(kte,kde-1)
135
136	i_start = its
137	i_end = MIN(ite,ide-1)
138	j_start = jts
139	j_end = MIN(jte,jde-1)
140
141	IF ( config_flags%open_xs .or. config_flags%specified .or. &
142	config_flags%nested) i_start = MAX(ids+1,its)
143	IF ( config_flags%open_xe .or. config_flags%specified .or. &
144	config_flags%nested) i_end = MIN(ide-2,ite)
145	IF ( config_flags%open_ys .or. config_flags%specified .or. &
146	config_flags%nested) j_start = MAX(jds+1,jts)
147	IF ( config_flags%open_ye .or. config_flags%specified .or. &
148	config_flags%nested) j_end = MIN(jde-2,jte)
149	IF ( config_flags%periodic_x ) i_start = its
150	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
151
152	!-----------------------------------------------------------------------------
153
154	! Below the 0.25 factor divides the incoming WRF deformations,
155	! which are multiplied by a factor of 2, by 2
156
157	DO j=j_start,j_end
158	DO k=kts,ktf
159	DO i=i_start,i_end
160
161	smnsmn(i,k,j) = 0.25( s11(i,k,j)s11(i,k,j) + &
162	s22(i,k,j)*s22(i,k,j) + &
163	s33(i,k,j)*s33(i,k,j) )
164
165	END DO
166	END DO
167	END DO
168
169	! Below the 0.125 factor accounts for the four-point averaging (0.25)
170	! and divides the incoming WRF deformation elements by 2 (0.5)
171
172	DO j=j_start,j_end
173	DO k=kts,ktf
174	DO i=i_start,i_end
175
176	tmp = 0.125*( s12(i ,k,j) + s12(i ,k,j+1) + &
177	s12(i+1,k,j) + s12(i+1,k,j+1) )
178	smnsmn(i,k,j) = smnsmn(i,k,j) + 2.0tmptmp
179
180	END DO
181	END DO
182	END DO
183
184	DO j=j_start,j_end
185	DO k=kts,ktf
186	DO i=i_start,i_end
187
188	tmp = 0.125*( s13(i ,k+1,j) + s13(i ,k,j) + &
189	s13(i+1,k+1,j) + s13(i+1,k,j) )
190	smnsmn(i,k,j) = smnsmn(i,k,j) + 2.0tmptmp
191
192	END DO
193	END DO
194	END DO
195
196	DO j=j_start,j_end
197	DO k=kts,ktf
198	DO i=i_start,i_end
199
200	tmp = 0.125*( s23(i,k+1,j ) + s23(i,k,j ) + &
201	s23(i,k+1,j+1) + s23(i,k,j+1) )
202	smnsmn(i,k,j) = smnsmn(i,k,j) + 2.0tmptmp
203
204	END DO
205	END DO
206	END DO
207
208	RETURN
209
210	END SUBROUTINE calc_smnsmn
211
212	!=============================================================================
213
214	SUBROUTINE calc_mii( m11, m22, m33, &
215	s11, s22, s33, &
216	s12, s13, s23, &
217	r12, r13, r23, smnsmn, &
218	tke, rdzw, dx, dy, &
219	config_flags, &
220	ids, ide, jds, jde, kds, kde, &
221	ims, ime, jms, jme, kms, kme, &
222	ips, ipe, jps, jpe, kps, kpe, &
223	its, ite, jts, jte, kts, kte )
224
225	!-----------------------------------------------------------------------------
226	!
227	! PURPOSE: Compute Mij for i = j
228	!
229	!-----------------------------------------------------------------------------
230
231	IMPLICIT NONE
232
233	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
234	:: m11 & ! NBA stress element 11 (m2 s-2)
235	, m22 & ! NBA stress element 22 (m2 s-2)
236	, m33 ! NBA stress element 33 (m2 s-2)
237
238	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
239	:: s11 & ! 2*deformation element 11 (s-1)
240	, s22 & ! 2*deformation element 22 (s-1)
241	, s33 & ! 2*deformation element 33 (s-1)
242	, s12 & ! 2*deformation element 12 (s-1)
243	, s13 & ! 2*deformation element 13 (s-1)
244	, s23 & ! 2*deformation element 23 (s-1)
245	, r12 & ! 2*rotation element 12 (s-1)
246	, r13 & ! 2*rotation element 13 (s-1)
247	, r23 & ! 2*rotation element 23 (s-1)
248	, smnsmn & ! Smn*Smn (s-2)
249	, tke & ! tke (m2 s-2)
250	, rdzw ! 1/dz at w-levels (m-1)
251
252	REAL, INTENT( IN ) &
253	:: dx & ! grid spacing in x (m)
254	, dy ! grid spacing in y (m)
255
256	TYPE (grid_config_rec_type), INTENT( IN ) &
257	:: config_flags
258
259	INTEGER, INTENT( IN ) &
260	:: ids, ide, jds, jde, kds, kde, &
261	ims, ime, jms, jme, kms, kme, &
262	ips, ipe, jps, jpe, kps, kpe, &
263	its, ite, jts, jte, kts, kte
264
265	! LOCAL VARIABLES ------------------------------------------------------------
266
267	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
268	:: ss11 &
269	, ss22 &
270	, ss33 &
271	, ss12 &
272	, ss13 &
273	, ss23 &
274	, rr12 &
275	, rr13 &
276	, rr23
277
278	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to c
279	:: ss12c &
280	, rr12c &
281	, ss13c &
282	, rr13c &
283	, ss23c &
284	, rr23c
285
286	REAL :: delta, a, b
287
288	INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, is_ext, js_ext
289
290	!-----------------------------------------------------------------------------
291	!
292	! Set loop limits,
293	! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_11_22_33
294	!
295	!-----------------------------------------------------------------------------
296
297	ktf = MIN( kte, kde-1 )
298
299	i_start = its
300	i_end = ite
301	j_start = jts
302	j_end = jte
303
304	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
305	config_flags%nested) i_start = MAX( ids+1, its )
306	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
307	config_flags%nested) i_end = MIN( ide-1, ite )
308	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
309	config_flags%nested) j_start = MAX( jds+1, jts )
310	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
311	config_flags%nested) j_end = MIN( jde-1, jte )
312	IF ( config_flags%periodic_x ) i_start = its
313	IF ( config_flags%periodic_x ) i_end = ite
314
315	is_ext = 1
316	js_ext = 1
317
318	i_start = i_start - is_ext
319	j_start = j_start - js_ext
320
321	!-----------------------------------------------------------------------------
322	!
323	! Divide WRF deformations, which are multiplied by 2, by 2
324	!
325	!-----------------------------------------------------------------------------
326
327	DO j=j_start,j_end+1
328	DO k=kts,ktf
329	DO i=i_start,i_end+1
330
331	ss11(i,k,j)=s11(i,k,j)/2.0
332	ss22(i,k,j)=s22(i,k,j)/2.0
333	ss33(i,k,j)=s33(i,k,j)/2.0
334	ss12(i,k,j)=s12(i,k,j)/2.0
335	ss13(i,k,j)=s13(i,k,j)/2.0
336	ss23(i,k,j)=s23(i,k,j)/2.0
337	rr12(i,k,j)=r12(i,k,j)/2.0
338	rr13(i,k,j)=r13(i,k,j)/2.0
339	rr23(i,k,j)=r23(i,k,j)/2.0
340
341	END DO
342	END DO
343	END DO
344
345	DO j=j_start,j_end+1
346	DO i=i_start,i_end+1
347
348	ss13(i,kde,j) = 0.0
349	ss23(i,kde,j) = 0.0
350	rr13(i,kde,j) = 0.0
351	rr23(i,kde,j) = 0.0
352
353	END DO
354	END DO
355
356	!-----------------------------------------------------------------------------
357	!
358	! Project variables to c
359	!
360	!-----------------------------------------------------------------------------
361
362	DO j = j_start, j_end
363	DO k = kts, ktf
364	DO i = i_start, i_end
365
366	ss12c(i,k,j) = 0.25*( ss12(i ,k ,j ) + ss12(i ,k ,j+1) + &
367	ss12(i+1,k ,j ) + ss12(i+1,k ,j+1) )
368
369	rr12c(i,k,j) = 0.25*( rr12(i ,k ,j ) + rr12(i ,k ,j+1) + &
370	rr12(i+1,k ,j ) + rr12(i+1,k ,j+1) )
371
372	ss13c(i,k,j) = 0.25*( ss13(i ,k+1,j ) + ss13(i ,k ,j ) + &
373	ss13(i+1,k+1,j ) + ss13(i+1,k ,j ) )
374
375	rr13c(i,k,j) = 0.25*( rr13(i ,k+1,j ) + rr13(i ,k ,j ) + &
376	rr13(i+1,k+1,j ) + rr13(i+1,k ,j ) )
377
378	ss23c(i,k,j) = 0.25*( ss23(i ,k+1,j ) + ss23(i ,k ,j ) + &
379	ss23(i ,k+1,j+1) + ss23(i ,k ,j+1) )
380
381	rr23c(i,k,j) = 0.25*( rr23(i ,k+1,j ) + rr23(i ,k ,j ) + &
382	rr23(i ,k+1,j+1) + rr23(i ,k ,j+1) )
383
384	ENDDO
385	ENDDO
386	ENDDO
387
388	!-----------------------------------------------------------------------------
389	!
390	! Calculate M11, M22 and M33
391	!
392	!-----------------------------------------------------------------------------
393
394	IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
395
396	DO j=j_start,j_end
397	DO k=kts,ktf
398	DO i=i_start,i_end
399
400	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
401	a = -1.0( csdelta )**2
402
403	m11(i,k,j) = a( 2.0sqrt( 2.0smnsmn(i,k,j) )ss11(i,k,j) &
404	+ c1( ss11(i,k,j) ss11(i,k,j) &
405	+ ss12c(i,k,j)*ss12c(i,k,j) &
406	+ ss13c(i,k,j)*ss13c(i,k,j) &
407	- smnsmn(i,k,j)/3.0 &
408	) &
409	+ c2( -2.0( ss12c(i,k,j)*rr12c(i,k,j) &
410	+ ss13c(i,k,j)*rr13c(i,k,j) &
411	) &
412	) &
413	)
414
415	m22(i,k,j) = a( 2.0sqrt( 2.0smnsmn(i,k,j) )ss22(i,k,j) &
416	+ c1( ss22(i,k,j) ss22(i,k,j) &
417	+ ss12c(i,k,j)*ss12c(i,k,j) &
418	+ ss23c(i,k,j)*ss23c(i,k,j) &
419	- smnsmn(i,k,j)/3.0 &
420	) &
421	+ c2( 2.0( ss12c(i,k,j)*rr12c(i,k,j) &
422	- ss23c(i,k,j)*rr23c(i,k,j) &
423	) &
424	) &
425	)
426
427	m33(i,k,j) = a( 2.0sqrt( 2.0smnsmn(i,k,j) )ss33(i,k,j) &
428	+ c1( ss33(i,k,j) ss33(i,k,j) &
429	+ ss13c(i,k,j)*ss13c(i,k,j) &
430	+ ss23c(i,k,j)*ss23c(i,k,j) &
431	- smnsmn(i,k,j)/3.0 &
432	) &
433	+ c2( 2.0( ss13c(i,k,j)*rr13c(i,k,j) &
434	+ ss23c(i,k,j)*rr23c(i,k,j) &
435	) &
436	) &
437	)
438
439	ENDDO
440	ENDDO
441	ENDDO
442
443	ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
444
445	DO j=j_start,j_end
446	DO k=kts,ktf
447	DO i=i_start,i_end
448
449	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
450	a = -1.0cedelta
451	b = c3*delta
452
453	m11(i,k,j) = a( 2.0sqrt( tke(i,k,j) )*ss11(i,k,j) &
454	+ b*( &
455	c1( ss11(i,k,j) ss11(i,k,j) &
456	+ ss12c(i,k,j)*ss12c(i,k,j) &
457	+ ss13c(i,k,j)*ss13c(i,k,j) &
458	- smnsmn(i,k,j)/3.0 &
459	) &
460	+ c2( -2.0( ss12c(i,k,j)*rr12c(i,k,j) &
461	+ ss13c(i,k,j)*rr13c(i,k,j) &
462	) &
463	) &
464	) &
465	)
466
467	m22(i,k,j) = a( 2.0sqrt( tke(i,k,j) )*ss22(i,k,j) &
468	+ b*( &
469	c1( ss22(i,k,j) ss22(i,k,j) &
470	+ ss12c(i,k,j)*ss12c(i,k,j) &
471	+ ss23c(i,k,j)*ss23c(i,k,j) &
472	- smnsmn(i,k,j)/3.0 &
473	) &
474	+ c2( 2.0( ss12c(i,k,j)*rr12c(i,k,j) &
475	- ss23c(i,k,j)*rr23c(i,k,j) &
476	) &
477	) &
478	) &
479	)
480
481	m33(i,k,j) = a( 2.0sqrt( tke(i,k,j) )*ss33(i,k,j) &
482	+ b*( &
483	c1( ss33(i,k,j) ss33(i,k,j) &
484	+ ss13c(i,k,j)*ss13c(i,k,j) &
485	+ ss23c(i,k,j)*ss23c(i,k,j) &
486	- smnsmn(i,k,j)/3.0 &
487	) &
488	+ c2( 2.0( ss13c(i,k,j)*rr13c(i,k,j) &
489	+ ss23c(i,k,j)*rr23c(i,k,j) &
490	) &
491	) &
492	) &
493	)
494
495
496	ENDDO
497	ENDDO
498	ENDDO
499
500	ENDIF
501
502	RETURN
503
504	END SUBROUTINE calc_mii
505
506	!=============================================================================
507
508	SUBROUTINE calc_m12( m12, &
509	s11, s22, &
510	s12, s13, s23, &
511	r12, r13, r23, smnsmn, &
512	tke, rdzw, dx, dy, &
513	config_flags, &
514	ids, ide, jds, jde, kds, kde, &
515	ims, ime, jms, jme, kms, kme, &
516	ips, ipe, jps, jpe, kps, kpe, &
517	its, ite, jts, jte, kts, kte )
518
519	!-----------------------------------------------------------------------------
520	!
521	! PURPOSE: Compute M12
522	!
523	!-----------------------------------------------------------------------------
524
525	IMPLICIT NONE
526
527	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
528	:: m12 ! NBA stress element 12 (m2 s-2)
529
530	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
531	:: s11 & ! 2*deformation element 11 (s-1)
532	, s22 & ! 2*deformation element 22 (s-1)
533	, s12 & ! 2*deformation element 12 (s-1)
534	, s13 & ! 2*deformation element 13 (s-1)
535	, s23 & ! 2*deformation element 23 (s-1)
536	, r12 & ! 2*rotation element 12 (s-1)
537	, r13 & ! 2*rotation element 13 (s-1)
538	, r23 & ! 2*rotation element 23 (s-1)
539	, smnsmn & ! Smn*Smn (s-2)
540	, tke & ! tke (m2 s-2)
541	, rdzw ! 1/dz at w-levels (m-1)
542
543	REAL, INTENT( IN ) &
544	:: dx & ! grid spacing in x (m)
545	, dy ! grid spacing in y (m)
546
547	TYPE (grid_config_rec_type), INTENT( IN ) &
548	:: config_flags
549
550	INTEGER, INTENT( IN ) &
551	:: ids, ide, jds, jde, kds, kde, &
552	ims, ime, jms, jme, kms, kme, &
553	ips, ipe, jps, jpe, kps, kpe, &
554	its, ite, jts, jte, kts, kte
555
556	! LOCAL VARIABLES ------------------------------------------------------------
557
558	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
559	:: ss11 &
560	, ss22 &
561	, ss12 &
562	, ss13 &
563	, ss23 &
564	, rr12 &
565	, rr13 &
566	, rr23
567
568
569	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to d
570	:: tked &
571	, ss11d &
572	, ss22d &
573	, ss13d &
574	, ss23d &
575	, rr13d &
576	, rr23d &
577	, smnsmnd
578
579	REAL :: delta, a, b
580
581	INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, je_ext, ie_ext
582
583	!-----------------------------------------------------------------------------
584	!
585	! Set loop limits,
586	! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_12_21
587	!
588	!-----------------------------------------------------------------------------
589
590	ktf = MIN( kte, kde-1 )
591
592	! Needs one more point in the x and y directions.
593
594	i_start = its
595	i_end = ite
596	j_start = jts
597	j_end = jte
598
599	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
600	config_flags%nested ) i_start = MAX( ids+1, its )
601	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
602	config_flags%nested ) i_end = MIN( ide-1, ite )
603	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
604	config_flags%nested ) j_start = MAX( jds+1, jts )
605	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
606	config_flags%nested ) j_end = MIN( jde-1, jte )
607	IF ( config_flags%periodic_x ) i_start = its
608	IF ( config_flags%periodic_x ) i_end = ite
609
610	je_ext = 1
611	ie_ext = 1
612
613	i_end = i_end + ie_ext
614	j_end = j_end + je_ext
615
616	!-----------------------------------------------------------------------------
617	!
618	! Divide WRF deformations, which are multiplied by 2, by 2
619	!
620	!-----------------------------------------------------------------------------
621
622	DO j=j_start-1,j_end
623	DO k=kts,ktf
624	DO i=i_start-1,i_end
625
626	ss11(i,k,j)=s11(i,k,j)/2.0
627	ss22(i,k,j)=s22(i,k,j)/2.0
628	ss12(i,k,j)=s12(i,k,j)/2.0
629	ss13(i,k,j)=s13(i,k,j)/2.0
630	ss23(i,k,j)=s23(i,k,j)/2.0
631	rr12(i,k,j)=r12(i,k,j)/2.0
632	rr13(i,k,j)=r13(i,k,j)/2.0
633	rr23(i,k,j)=r23(i,k,j)/2.0
634
635	END DO
636	END DO
637	END DO
638
639	DO j=j_start-1,j_end
640	DO i=i_start-1,i_end
641
642	ss13(i,kde,j) = 0.0
643	ss23(i,kde,j) = 0.0
644	rr13(i,kde,j) = 0.0
645	rr23(i,kde,j) = 0.0
646
647	END DO
648	END DO
649
650	!-----------------------------------------------------------------------------
651	!
652	! Project variables to d
653	!
654	!-----------------------------------------------------------------------------
655
656	DO j = j_start, j_end
657	DO k = kts, ktf
658	DO i = i_start, i_end
659
660	tked(i,k,j) = 0.25*( tke(i-1,k ,j ) + tke(i ,k ,j ) + &
661	tke(i-1,k ,j-1) + tke(i ,k ,j-1) )
662
663	smnsmnd(i,k,j) = 0.25*( smnsmn(i-1,k ,j ) + smnsmn(i ,k ,j ) + &
664	smnsmn(i-1,k ,j-1) + smnsmn(i ,k ,j-1) )
665
666	ss11d(i,k,j) = 0.25*( ss11(i-1,k ,j ) + ss11(i ,k ,j ) + &
667	ss11(i-1,k ,j-1) + ss11(i ,k ,j-1) )
668
669	ss22d(i,k,j) = 0.25*( ss22(i-1,k ,j ) + ss22(i ,k ,j ) + &
670	ss22(i-1,k ,j-1) + ss22(i ,k ,j-1) )
671
672	ss13d(i,k,j) = 0.25*( ss13(i ,k+1,j ) + ss13(i ,k+1,j-1) + &
673	ss13(i ,k ,j ) + ss13(i ,k ,j-1) )
674
675	rr13d(i,k,j) = 0.25*( rr13(i ,k+1,j ) + rr13(i ,k+1,j-1) + &
676	rr13(i ,k ,j ) + rr13(i ,k ,j-1) )
677
678	ss23d(i,k,j) = 0.25*( ss23(i ,k+1,j ) + ss23(i-1,k+1,j ) + &
679	ss23(i ,k ,j ) + ss23(i-1,k ,j ) )
680
681	rr23d(i,k,j) = 0.25*( rr23(i ,k+1,j ) + rr23(i-1,k+1,j ) + &
682	rr23(i ,k ,j ) + rr23(i-1,k ,j ) )
683
684	END DO
685	END DO
686	END DO
687
688	!-----------------------------------------------------------------------------
689	!
690	! Calculate M12
691	!
692	!-----------------------------------------------------------------------------
693
694	IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
695
696	DO j=j_start,j_end
697	DO k=kts,ktf
698	DO i=i_start,i_end
699
700	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
701	a = -1.0( csdelta )**2
702
703	m12(i,k,j) = a( 2.0sqrt( 2.0smnsmnd(i,k,j) )ss12(i,k,j) &
704	+ c1( ss11d(i,k,j)ss12(i,k,j) &
705	+ ss22d(i,k,j)*ss12(i,k,j) &
706	+ ss13d(i,k,j)*ss23d(i,k,j) &
707	) &
708	+ c2( ss11d(i,k,j)rr12(i,k,j) &
709	- ss13d(i,k,j)*rr23d(i,k,j) &
710	- ss22d(i,k,j)*rr12(i,k,j) &
711	- ss23d(i,k,j)*rr13d(i,k,j) &
712	) &
713	)
714
715	ENDDO
716	ENDDO
717	ENDDO
718
719	ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
720
721	DO j=j_start,j_end
722	DO k=kts,ktf
723	DO i=i_start,i_end
724
725	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
726	a = -1.0cedelta
727	b = c3*delta
728
729	m12(i,k,j) = a( 2.0sqrt( tked(i,k,j) )*s12(i,k,j) &
730	+ b*( &
731	c1( ss11d(i,k,j)ss12(i,k,j) &
732	+ ss22d(i,k,j)*ss12(i,k,j) &
733	+ ss13d(i,k,j)*ss23d(i,k,j) &
734	) &
735	+ c2( ss11d(i,k,j)rr12(i,k,j) &
736	- ss13d(i,k,j)*rr23d(i,k,j) &
737	- ss22d(i,k,j)*rr12(i,k,j) &
738	- ss23d(i,k,j)*rr13d(i,k,j) &
739	) &
740	) &
741	)
742	ENDDO
743	ENDDO
744	ENDDO
745
746	ENDIF
747
748	RETURN
749
750	END SUBROUTINE calc_m12
751
752	!=============================================================================
753
754	SUBROUTINE calc_m13( m13, &
755	s11, s33, &
756	s12, s13, s23, &
757	r12, r13, r23, smnsmn, &
758	tke, rdzw, dx, dy, &
759	fnm, fnp, &
760	config_flags, &
761	ids, ide, jds, jde, kds, kde, &
762	ims, ime, jms, jme, kms, kme, &
763	ips, ipe, jps, jpe, kps, kpe, &
764	its, ite, jts, jte, kts, kte )
765
766	!-----------------------------------------------------------------------------
767	!
768	! PURPOSE: Compute M13
769	!
770	!-----------------------------------------------------------------------------
771
772	IMPLICIT NONE
773
774	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
775	:: m13 ! (m2 s-2)
776
777	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
778	:: s11 & ! 2*deformation element 11 (s-1)
779	, s33 & ! 2*deformation element 33 (s-1)
780	, s12 & ! 2*deformation element 12 (s-1)
781	, s13 & ! 2*deformation element 13 (s-1)
782	, s23 & ! 2*deformation element 23 (s-1)
783	, r12 & ! 2*rotation element 12 (s-1)
784	, r13 & ! 2*rotation element 13 (s-1)
785	, r23 & ! 2*rotation element 23 (s-1)
786	, smnsmn & ! Smn*Smn (s-2)
787	, tke & ! tke (m2 s-2)
788	, rdzw ! 1/dz at w-levels (m-1)
789
790	REAL, INTENT( IN ) &
791	:: dx & ! grid spacing in x (m)
792	, dy ! grid spacing in y (m)
793
794	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
795	:: fnm & ! vertical interpolation coefficients
796	, fnp !
797
798	TYPE (grid_config_rec_type), INTENT( IN ) &
799	:: config_flags
800
801	INTEGER, INTENT( IN ) &
802	:: ids, ide, jds, jde, kds, kde, &
803	ims, ime, jms, jme, kms, kme, &
804	ips, ipe, jps, jpe, kps, kpe, &
805	its, ite, jts, jte, kts, kte
806
807	! LOCAL VARIABLES ------------------------------------------------------------
808
809	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
810	:: ss11 &
811	, ss33 &
812	, ss12 &
813	, ss13 &
814	, ss23 &
815	, rr12 &
816	, rr13 &
817	, rr23
818
819	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to e
820	:: tkee &
821	, ss11e &
822	, ss33e &
823	, ss12e &
824	, ss23e &
825	, rr12e &
826	, rr23e &
827	, smnsmne
828
829	REAL :: delta, a, b
830
831	INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, ie_ext
832
833	!-----------------------------------------------------------------------------
834	!
835	! Set loop limits,
836	! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_13_31
837	!
838	!-----------------------------------------------------------------------------
839
840	ktf = MIN( kte, kde-1 )
841
842	! Find ide-1 and jde-1 for averaging to p point.
843
844	i_start = its
845	i_end = ite
846	j_start = jts
847	j_end = MIN( jte, jde-1 )
848
849	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
850	config_flags%nested) i_start = MAX( ids+1, its )
851	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
852	config_flags%nested) i_end = MIN( ide-1, ite )
853	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
854	config_flags%nested) j_start = MAX( jds+1, jts )
855	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
856	config_flags%nested) j_end = MIN( jde-2, jte )
857	IF ( config_flags%periodic_x ) i_start = its
858	IF ( config_flags%periodic_x ) i_end = ite
859
860	ie_ext = 1
861	i_end = i_end + ie_ext
862
863	!-----------------------------------------------------------------------------
864	!
865	! Divide WRF deformations, which are multiplied by 2, by 2
866	!
867	!-----------------------------------------------------------------------------
868
869	DO j=j_start,j_end+1
870	DO k=kts,ktf
871	DO i=i_start-1,i_end
872
873	ss11(i,k,j)=s11(i,k,j)/2.0
874	ss33(i,k,j)=s33(i,k,j)/2.0
875	ss12(i,k,j)=s12(i,k,j)/2.0
876	ss13(i,k,j)=s13(i,k,j)/2.0
877	ss23(i,k,j)=s23(i,k,j)/2.0
878	rr12(i,k,j)=r12(i,k,j)/2.0
879	rr13(i,k,j)=r13(i,k,j)/2.0
880	rr23(i,k,j)=r23(i,k,j)/2.0
881
882	END DO
883	END DO
884	END DO
885
886	!-----------------------------------------------------------------------------
887	!
888	! Project variables to e
889	!
890	!-----------------------------------------------------------------------------
891
892	DO j = j_start, j_end
893	DO k = kts+1, ktf
894	DO i = i_start, i_end
895
896	tkee(i,k,j) = 0.5( fnm(k)( tke(i,k ,j) + tke(i-1,k ,j) ) + &
897	fnp(k)*( tke(i,k-1,j) + tke(i-1,k-1,j) ) )
898
899	smnsmne(i,k,j) = 0.5( fnm(k)( smnsmn(i,k ,j) + smnsmn(i-1,k ,j) ) + &
900	fnp(k)*( smnsmn(i,k-1,j) + smnsmn(i-1,k-1,j) ) )
901
902	ss11e(i,k,j) = 0.5( fnm(k)( ss11(i ,k ,j ) + ss11(i-1,k ,j ) ) + &
903	fnp(k)*( ss11(i ,k-1,j ) + ss11(i-1,k-1,j ) ) )
904
905	ss33e(i,k,j) = 0.5( fnm(k)( ss33(i ,k ,j ) + ss33(i-1,k ,j ) ) + &
906	fnp(k)*( ss33(i ,k-1,j ) + ss33(i-1,k-1,j ) ) )
907
908	ss12e(i,k,j) = 0.5( fnm(k)( ss12(i ,k ,j ) + ss12(i ,k ,j+1) ) + &
909	fnp(k)*( ss12(i ,k-1,j ) + ss12(i ,k-1,j+1) ) )
910
911	rr12e(i,k,j) = 0.5( fnm(k)( rr12(i ,k ,j ) + rr12(i ,k ,j+1) ) + &
912	fnp(k)*( rr12(i ,k-1,j ) + rr12(i ,k-1,j+1) ) )
913
914	ss23e(i,k,j) = 0.25*( ss23(i ,k ,j) + ss23(i ,k ,j+1) + &
915	ss23(i-1,k ,j) + ss23(i-1,k ,j+1) )
916
917	rr23e(i,k,j) = 0.25*( rr23(i ,k ,j) + rr23(i ,k ,j+1) + &
918	rr23(i-1,k ,j) + rr23(i-1,k ,j+1) )
919
920	END DO
921	END DO
922	END DO
923
924	!-----------------------------------------------------------------------------
925	!
926	! Calculate M_13
927	!
928	!-----------------------------------------------------------------------------
929
930
931	IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
932
933	DO j=j_start,j_end
934	DO k=kts+1,ktf
935	DO i=i_start,i_end
936
937	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
938	a = -1.0( csdelta )**2
939
940	m13(i,k,j) = a( 2.0sqrt( 2.0smnsmne(i,k,j) )ss13(i,k,j) &
941	+ c1( ss11e(i,k,j)ss13(i,k,j) &
942	+ ss12e(i,k,j)*ss23e(i,k,j) &
943	+ ss13(i,k,j)*ss33e(i,k,j) &
944	) &
945	+ c2( ss11e(i,k,j)rr13(i,k,j) &
946	+ ss12e(i,k,j)*rr23e(i,k,j) &
947	- ss23e(i,k,j)*rr12e(i,k,j) &
948	- ss33e(i,k,j)*rr13(i,k,j) &
949	) &
950	)
951
952	ENDDO
953	ENDDO
954	ENDDO
955
956	ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
957
958	DO j=j_start,j_end
959	DO k=kts+1,ktf
960	DO i=i_start,i_end
961
962	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
963	a = -1.0cedelta
964	b = c3*delta
965
966	m13(i,k,j) = a( 2.0sqrt( tkee(i,k,j) )*ss13(i,k,j) &
967	+ b*( &
968	c1( ss11e(i,k,j)ss13(i,k,j) &
969	+ ss12e(i,k,j)*ss23e(i,k,j) &
970	+ ss13(i,k,j)*ss33e(i,k,j) &
971	) &
972	+ c2( ss11e(i,k,j)rr13(i,k,j) &
973	+ ss12e(i,k,j)*rr23e(i,k,j) &
974	- ss23e(i,k,j)*rr12e(i,k,j) &
975	- ss33e(i,k,j)*rr13(i,k,j) &
976	) &
977	) &
978	)
979
980	ENDDO
981	ENDDO
982	ENDDO
983
984	ENDIF
985
986	RETURN
987
988	END SUBROUTINE calc_m13
989
990	!=============================================================================
991
992	SUBROUTINE calc_m23( m23, &
993	s22, s33, &
994	s12, s13, s23, &
995	r12, r13, r23, smnsmn, &
996	tke, rdzw, dx, dy, &
997	fnm, fnp, &
998	config_flags, &
999	ids, ide, jds, jde, kds, kde, &
1000	ims, ime, jms, jme, kms, kme, &
1001	ips, ipe, jps, jpe, kps, kpe, &
1002	its, ite, jts, jte, kts, kte )
1003
1004	!-----------------------------------------------------------------------------
1005	!
1006	! PURPOSE: Compute M23
1007	!
1008	!-----------------------------------------------------------------------------
1009
1010	IMPLICIT NONE
1011
1012	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
1013	:: m23 ! (m2 s-2)
1014
1015	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
1016	:: s22 & ! 2*deformation element 22 (s-1)
1017	, s33 & ! 2*deformation element 33 (s-1)
1018	, s12 & ! 2*deformation element 12 (s-1)
1019	, s13 & ! 2*deformation element 13 (s-1)
1020	, s23 & ! 2*deformation element 23 (s-1)
1021	, r12 & ! 2*rotation element 12 (s-1)
1022	, r13 & ! 2*rotation element 13 (s-1)
1023	, r23 & ! 2*rotation element 23 (s-1)
1024	, smnsmn & ! Smn*Smn (s-2)
1025	, tke & ! tke (m2 s-2)
1026	, rdzw ! 1/dz at w-levels (m-1)
1027
1028	REAL, INTENT( IN ) &
1029	:: dx & ! grid spacing in x (m)
1030	, dy ! grid spacing in y (m)
1031
1032	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
1033	:: fnm & ! vertical interpolation coefficients
1034	, fnp !
1035
1036	TYPE (grid_config_rec_type), INTENT( IN ) &
1037	:: config_flags
1038
1039	INTEGER, INTENT( IN ) &
1040	:: ids, ide, jds, jde, kds, kde, &
1041	ims, ime, jms, jme, kms, kme, &
1042	ips, ipe, jps, jpe, kps, kpe, &
1043	its, ite, jts, jte, kts, kte
1044
1045	! LOCAL VARIABLES ------------------------------------------------------------
1046
1047	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
1048	:: ss22 &
1049	, ss33 &
1050	, ss12 &
1051	, ss13 &
1052	, ss23 &
1053	, rr12 &
1054	, rr13 &
1055	, rr23
1056
1057	REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to f
1058	:: tkef &
1059	, ss22f &
1060	, ss33f &
1061	, ss12f &
1062	, ss13f &
1063	, rr12f &
1064	, rr13f &
1065	, smnsmnf
1066
1067	REAL :: delta, a, b
1068
1069	INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, je_ext
1070
1071	!-----------------------------------------------------------------------------
1072	!
1073	! Set loop limits,
1074	! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_23_32
1075	!
1076	!-----------------------------------------------------------------------------
1077
1078	ktf = MIN( kte, kde-1 )
1079
1080	! Find ide-1 and jde-1 for averaging to p point.
1081
1082	i_start = its
1083	i_end = MIN( ite, ide-1 )
1084	j_start = jts
1085	j_end = jte
1086
1087	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
1088	config_flags%nested) i_start = MAX( ids+1, its )
1089	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
1090	config_flags%nested) i_end = MIN( ide-2, ite )
1091	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
1092	config_flags%nested) j_start = MAX( jds+1, jts )
1093	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
1094	config_flags%nested) j_end = MIN( jde-1, jte )
1095	IF ( config_flags%periodic_x ) i_start = its
1096	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1097
1098	je_ext = 1
1099	j_end = j_end + je_ext
1100
1101	!-----------------------------------------------------------------------------
1102	!
1103	! Divide WRF deformations, which are multiplied by 2, by 2
1104	!
1105	!-----------------------------------------------------------------------------
1106
1107	DO j=j_start-1,j_end
1108	DO k=kts,ktf
1109	DO i=i_start,i_end+1
1110
1111	ss22(i,k,j)=s22(i,k,j)/2.0
1112	ss33(i,k,j)=s33(i,k,j)/2.0
1113	ss12(i,k,j)=s12(i,k,j)/2.0
1114	ss13(i,k,j)=s13(i,k,j)/2.0
1115	ss23(i,k,j)=s23(i,k,j)/2.0
1116	rr12(i,k,j)=r12(i,k,j)/2.0
1117	rr13(i,k,j)=r13(i,k,j)/2.0
1118	rr23(i,k,j)=r23(i,k,j)/2.0
1119
1120	END DO
1121	END DO
1122	END DO
1123
1124	!-----------------------------------------------------------------------------
1125	!
1126	! Project variables to f
1127	!
1128	!-----------------------------------------------------------------------------
1129
1130	DO j = j_start, j_end
1131	DO k = kts+1, ktf
1132	DO i = i_start, i_end
1133
1134	tkef(i,k,j) = 0.5( fnm(k)( tke(i ,k ,j ) + tke(i ,k ,j-1) ) + &
1135	fnp(k)*( tke(i ,k-1,j ) + tke(i ,k-1,j-1) ) )
1136
1137	smnsmnf(i,k,j) = 0.5( fnm(k)( smnsmn(i ,k ,j ) + smnsmn(i ,k ,j-1) ) + &
1138	fnp(k)*( smnsmn(i ,k-1,j ) + smnsmn(i ,k-1,j-1) ) )
1139
1140	ss22f(i,k,j) = 0.5( fnm(k)( ss22(i ,k ,j ) + ss22(i ,k ,j-1) ) + &
1141	fnp(k)*( ss22(i ,k-1,j ) + ss22(i ,k-1,j-1) ) )
1142
1143	ss33f(i,k,j) = 0.5( fnm(k)( ss33(i ,k ,j ) + ss33(i ,k ,j-1) ) + &
1144	fnp(k)*( ss33(i ,k-1,j ) + ss33(i ,k-1,j-1) ) )
1145
1146	ss12f(i,k,j) = 0.5( fnm(k)( ss12(i ,k ,j ) + ss12(i+1,k ,j ) ) + &
1147	fnp(k)*( ss12(i ,k-1,j ) + ss12(i+1,k-1,j ) ) )
1148
1149	rr12f(i,k,j) = 0.5( fnm(k)( rr12(i ,k ,j ) + rr12(i+1,k ,j ) ) + &
1150	fnp(k)*( rr12(i ,k-1,j ) + rr12(i+1,k-1,j ) ) )
1151
1152	ss13f(i,k,j) = 0.25*( ss13(i ,k ,j ) + ss13(i ,k ,j-1) + &
1153	ss13(i+1,k ,j-1) + ss13(i+1,k ,j ) )
1154
1155	rr13f(i,k,j) = 0.25*( rr13(i ,k ,j ) + rr13(i ,k ,j-1) + &
1156	rr13(i+1,k ,j-1) + rr13(i+1,k ,j ) )
1157
1158	END DO
1159	END DO
1160	END DO
1161
1162	!-----------------------------------------------------------------------------
1163	!
1164	! Calculate M23
1165	!
1166	!-----------------------------------------------------------------------------
1167
1168	IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
1169
1170	DO j=j_start,j_end
1171	DO k=kts+1,ktf
1172	DO i=i_start,i_end
1173
1174	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
1175	a = -1.0( csdelta )**2
1176
1177	m23(i,k,j) = a( 2.0sqrt( 2.0smnsmnf(i,k,j) )ss23(i,k,j) &
1178	+ c1( ss12f(i,k,j)ss13f(i,k,j) &
1179	+ ss22f(i,k,j)*ss23(i,k,j) &
1180	+ ss23(i,k,j) *ss33f(i,k,j) &
1181	) &
1182	+ c2( ss12f(i,k,j)rr13f(i,k,j) &
1183	+ ss22f(i,k,j)*rr23(i,k,j) &
1184	+ ss13f(i,k,j)*rr12f(i,k,j) &
1185	- ss33f(i,k,j)*rr23(i,k,j) &
1186	) &
1187	)
1188
1189	ENDDO
1190	ENDDO
1191	ENDDO
1192
1193	ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
1194
1195	DO j=j_start,j_end
1196	DO k=kts+1,ktf
1197	DO i=i_start,i_end
1198
1199	delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
1200	a = -1.0cedelta
1201	b = c3*delta
1202
1203	m23(i,k,j) = a( 2.0sqrt( tkef(i,k,j) )*ss23(i,k,j) &
1204	+ b*( &
1205	c1( ss12f(i,k,j)ss13f(i,k,j) &
1206	+ ss22f(i,k,j)*ss23(i,k,j) &
1207	+ ss23(i,k,j) *ss33f(i,k,j) &
1208	) &
1209	+ c2( ss12f(i,k,j)rr13f(i,k,j) &
1210	+ ss22f(i,k,j)*rr23(i,k,j) &
1211	+ ss13f(i,k,j)*rr12f(i,k,j) &
1212	- ss33f(i,k,j)*rr23(i,k,j) &
1213	) &
1214	) &
1215	)
1216
1217	ENDDO
1218	ENDDO
1219	ENDDO
1220
1221	ENDIF
1222
1223	RETURN
1224
1225	END SUBROUTINE calc_m23
1226
1227	!=============================================================================
1228
1229	END MODULE module_sfs_nba

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

Download in other formats:

Original Format