Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

traldf_iso.F90 @ 13149

Last change on this file since 13149 was 13149, checked in by andmirek, 4 years ago
Ticket #2482: Dissable restart and use allocatable arrays
File size: 21.5 KB

Line
1	MODULE traldf_iso
2	!!======================================================================
3	!! * MODULE traldf_iso *
4	!! Ocean tracers: horizontal component of the lateral tracer mixing trend
5	!!======================================================================
6	!! History : OPA ! 1994-08 (G. Madec, M. Imbard)
7	!! 8.0 ! 1997-05 (G. Madec) split into traldf and trazdf
8	!! NEMO ! 2002-08 (G. Madec) Free form, F90
9	!! 1.0 ! 2005-11 (G. Madec) merge traldf and trazdf :-)
10	!! 3.3 ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC
11	!! 3.7 ! 2014-01 (G. Madec, S. Masson) restructuration/simplification of aht/aeiv specification
12	!! - ! 2014-02 (F. Lemarie, G. Madec) triad operator (Griffies) + Method of Stabilizing Correction
13	!!----------------------------------------------------------------------
14
15	!!----------------------------------------------------------------------
16	!! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator
17	!! and with the vertical part of the isopycnal or geopotential s-coord. operator
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and active tracers
20	USE dom_oce ! ocean space and time domain
21	USE trc_oce ! share passive tracers/Ocean variables
22	USE zdf_oce ! ocean vertical physics
23	USE ldftra ! lateral diffusion: tracer eddy coefficients
24	USE ldfslp ! iso-neutral slopes
25	USE diaptr ! poleward transport diagnostics
26	USE diaar5 ! AR5 diagnostics
27	!
28	USE in_out_manager ! I/O manager
29	USE iom ! I/O library
30	USE phycst ! physical constants
31	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
32
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC tra_ldf_iso ! routine called by step.F90
37
38	LOGICAL :: l_ptr ! flag to compute poleward transport
39	LOGICAL :: l_hst ! flag to compute heat transport
40	REAL(wp), DIMENSION(:, :) , ALLOCATABLE :: zdkt, zdk1t, z2d
41	REAL(wp), DIMENSION(:, :, :), ALLOCATABLE :: zdit, zdjt, zftu, zftv, ztfw
42
43	!! * Substitutions
44	# include "vectopt_loop_substitute.h90"
45	!!----------------------------------------------------------------------
46	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
47	!! $Id$
48	!! Software governed by the CeCILL license (see ./LICENSE)
49	!!----------------------------------------------------------------------
50	CONTAINS
51
52	SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv , &
53	& pgui, pgvi, &
54	& ptb , ptbb, pta , kjpt, kpass )
55	!!----------------------------------------------------------------------
56	!! * ROUTINE tra_ldf_iso *
57	!!
58	!! ** Purpose : Compute the before horizontal tracer (t & s) diffusive
59	!! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
60	!! add it to the general trend of tracer equation.
61	!!
62	!! ** Method : The horizontal component of the lateral diffusive trends
63	!! is provided by a 2nd order operator rotated along neural or geopo-
64	!! tential surfaces to which an eddy induced advection can be added
65	!! It is computed using before fields (forward in time) and isopyc-
66	!! nal or geopotential slopes computed in routine ldfslp.
67	!!
68	!! 1st part : masked horizontal derivative of T ( di[ t ] )
69	!! ======== with partial cell update if ln_zps=T
70	!! with top cell update if ln_isfcav
71	!!
72	!! 2nd part : horizontal fluxes of the lateral mixing operator
73	!! ========
74	!! zftu = pahu e2u*e3u/e1u di[ tb ]
75	!! - pahu e2u*uslp dk[ mi(mk(tb)) ]
76	!! zftv = pahv e1v*e3v/e2v dj[ tb ]
77	!! - pahv e2u*vslp dk[ mj(mk(tb)) ]
78	!! take the horizontal divergence of the fluxes:
79	!! difft = 1/(e1e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] }
80	!! Add this trend to the general trend (ta,sa):
81	!! ta = ta + difft
82	!!
83	!! 3rd part: vertical trends of the lateral mixing operator
84	!! ======== (excluding the vertical flux proportional to dk[t] )
85	!! vertical fluxes associated with the rotated lateral mixing:
86	!! zftw = - { mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
87	!! + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ] }
88	!! take the horizontal divergence of the fluxes:
89	!! difft = 1/(e1e2t*e3t) dk[ zftw ]
90	!! Add this trend to the general trend (ta,sa):
91	!! pta = pta + difft
92	!!
93	!! ** Action : Update pta arrays with the before rotated diffusion
94	!!----------------------------------------------------------------------
95	INTEGER , INTENT(in ) :: kt ! ocean time-step index
96	INTEGER , INTENT(in ) :: kit000 ! first time step index
97	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
98	INTEGER , INTENT(in ) :: kjpt ! number of tracers
99	INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage
100	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s]
101	REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels
102	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels
103	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! tracer (kpass=1) or laplacian of tracer (kpass=2)
104	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptbb ! tracer (only used in kpass=2)
105	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
106	!
107	INTEGER :: ji, jj, jk, jn ! dummy loop indices
108	INTEGER :: ikt
109	INTEGER :: ierr ! local integer
110	REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars
111	REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - -
112	REAL(wp) :: zcoef0, ze3w_2, zsign, z2dt, z1_2dt ! - -
113	! REAL(wp), DIMENSION(jpi,jpj) :: zdkt, zdk1t, z2d
114	! REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw
115	!!----------------------------------------------------------------------
116	!
117	IF( kpass == 1 .AND. kt == kit000 ) THEN
118	IF(lwp) WRITE(numout,*)
119	IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
120	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
121	!
122	akz (:,:,:) = 0._wp
123	ah_wslp2(:,:,:) = 0._wp
124	ALLOCATE(zdit(jpi,jpj,jpk), zdjt(jpi,jpj,jpk), zftu(jpi,jpj,jpk), zftv(jpi,jpj,jpk), ztfw(jpi,jpj,jpk))
125	ALLOCATE(zdkt(jpi,jpj), zdk1t(jpi,jpj), z2d(jpi,jpj))
126	ENDIF
127	!
128	l_hst = .FALSE.
129	l_ptr = .FALSE.
130	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
131	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
132	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
133	!
134	! ! set time step size (Euler/Leapfrog)
135	IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler)
136	ELSE ; z2dt = 2.* rdt ! (Leapfrog)
137	ENDIF
138	z1_2dt = 1._wp / z2dt
139	!
140	IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0)
141	ELSE ; zsign = -1._wp
142	ENDIF
143
144	!!----------------------------------------------------------------------
145	!! 0 - calculate ah_wslp2 and akz
146	!!----------------------------------------------------------------------
147	!
148	IF( kpass == 1 ) THEN !== first pass only ==!
149	!
150	DO jk = 2, jpkm1
151	DO jj = 2, jpjm1
152	DO ji = fs_2, fs_jpim1 ! vector opt.
153	!
154	zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &
155	& + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp )
156	zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &
157	& + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp )
158	!
159	zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) &
160	& + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku
161	zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) &
162	& + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv
163	!
164	ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
165	& + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
166	END DO
167	END DO
168	END DO
169	!
170	IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient
171	DO jk = 2, jpkm1
172	DO jj = 2, jpjm1
173	DO ji = fs_2, fs_jpim1
174	akz(ji,jj,jk) = 0.25_wp * ( &
175	& ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) &
176	& + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) &
177	& + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) &
178	& + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) )
179	END DO
180	END DO
181	END DO
182	!
183	IF( ln_traldf_blp ) THEN ! bilaplacian operator
184	DO jk = 2, jpkm1
185	DO jj = 1, jpjm1
186	DO ji = 1, fs_jpim1
187	akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) &
188	& * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) )
189	END DO
190	END DO
191	END DO
192	ELSEIF( ln_traldf_lap ) THEN ! laplacian operator
193	DO jk = 2, jpkm1
194	DO jj = 1, jpjm1
195	DO ji = 1, fs_jpim1
196	ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
197	zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 )
198	akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
199	END DO
200	END DO
201	END DO
202	ENDIF
203	!
204	ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
205	akz(:,:,:) = ah_wslp2(:,:,:)
206	ENDIF
207	ENDIF
208	!
209	! ! ===========
210	DO jn = 1, kjpt ! tracer loop
211	! ! ===========
212	!
213	!!----------------------------------------------------------------------
214	!! I - masked horizontal derivative
215	!!----------------------------------------------------------------------
216	!!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient....
217	zdit (1,:,:) = 0._wp ; zdit (jpi,:,:) = 0._wp
218	zdjt (1,:,:) = 0._wp ; zdjt (jpi,:,:) = 0._wp
219	!!end
220
221	! Horizontal tracer gradient
222	DO jk = 1, jpkm1
223	DO jj = 1, jpjm1
224	DO ji = 1, fs_jpim1 ! vector opt.
225	zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
226	zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
227	END DO
228	END DO
229	END DO
230	IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient
231	DO jj = 1, jpjm1 ! bottom correction (partial bottom cell)
232	DO ji = 1, fs_jpim1 ! vector opt.
233	zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
234	zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
235	END DO
236	END DO
237	IF( ln_isfcav ) THEN ! first wet level beneath a cavity
238	DO jj = 1, jpjm1
239	DO ji = 1, fs_jpim1 ! vector opt.
240	IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)
241	IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)
242	END DO
243	END DO
244	ENDIF
245	ENDIF
246	!
247	!!----------------------------------------------------------------------
248	!! II - horizontal trend (full)
249	!!----------------------------------------------------------------------
250	!
251	DO jk = 1, jpkm1 ! Horizontal slab
252	!
253	! !== Vertical tracer gradient
254	zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1
255	!
256	IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2)
257	ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
258	ENDIF
259	DO jj = 1 , jpjm1 !== Horizontal fluxes
260	DO ji = 1, fs_jpim1 ! vector opt.
261	zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
262	zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
263	!
264	zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) &
265	& + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. )
266	!
267	zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) &
268	& + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. )
269	!
270	zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
271	zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
272	!
273	zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) &
274	& + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) &
275	& + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk)
276	zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) &
277	& + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) &
278	& + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk)
279	END DO
280	END DO
281	!
282	DO jj = 2 , jpjm1 !== horizontal divergence and add to pta
283	DO ji = fs_2, fs_jpim1 ! vector opt.
284	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) &
285	& + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) &
286	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
287	END DO
288	END DO
289	END DO ! End of slab
290
291	!!----------------------------------------------------------------------
292	!! III - vertical trend (full)
293	!!----------------------------------------------------------------------
294	!
295	ztfw(1,:,:) = 0._wp ; ztfw(jpi,:,:) = 0._wp
296	!
297	! Vertical fluxes
298	! ---------------
299	! ! Surface and bottom vertical fluxes set to zero
300	ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp
301
302	DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1)
303	DO jj = 2, jpjm1
304	DO ji = fs_2, fs_jpim1 ! vector opt.
305	!
306	zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &
307	& + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp )
308	zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &
309	& + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp )
310	!
311	zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) &
312	& + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku
313	zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) &
314	& + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv
315	!
316	zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked
317	zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
318	!
319	ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) &
320	& + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) &
321	& + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) &
322	& + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) )
323	END DO
324	END DO
325	END DO
326	! !== add the vertical 33 flux ==!
327	IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz
328	DO jk = 2, jpkm1
329	DO jj = 1, jpjm1
330	DO ji = fs_2, fs_jpim1
331	ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) &
332	& * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) &
333	& * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
334	END DO
335	END DO
336	END DO
337	!
338	ELSE ! bilaplacian
339	SELECT CASE( kpass )
340	CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2
341	DO jk = 2, jpkm1
342	DO jj = 1, jpjm1
343	DO ji = fs_2, fs_jpim1
344	ztfw(ji,jj,jk) = ztfw(ji,jj,jk) &
345	& + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) &
346	& * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
347	END DO
348	END DO
349	END DO
350	CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp.
351	DO jk = 2, jpkm1
352	DO jj = 1, jpjm1
353	DO ji = fs_2, fs_jpim1
354	ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) &
355	& * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) &
356	& + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) )
357	END DO
358	END DO
359	END DO
360	END SELECT
361	ENDIF
362	!
363	DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==!
364	DO jj = 2, jpjm1
365	DO ji = fs_2, fs_jpim1 ! vector opt.
366	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) &
367	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
368	END DO
369	END DO
370	END DO
371	!
372	IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==!
373	( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass (bilaplacian) ==!
374	!
375	! ! "Poleward" diffusive heat or salt transports (T-S case only)
376	! note sign is reversed to give down-gradient diffusive transports )
377	IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:) )
378	! ! Diffusive heat transports
379	IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
380	!
381	ENDIF !== end pass selection ==!
382	!
383	! ! ===============
384	END DO ! end tracer loop
385	!
386	END SUBROUTINE tra_ldf_iso
387
388	!!==============================================================================
389	END MODULE traldf_iso

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE/TRA/traldf_iso.F90 @ 13149

Download in other formats: