File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\sicdif_sig.f
1 subroutine sicdife_sig(de,te,qe,xe,ye,ze,am,bm,tov,sv,dt,ue,ve,
2 x ls_node,
3 x snnp1ev,ndexev,locl,de_n)
4 use gfs_dyn_resol_def
5 use gfs_dyn_layout1
6 use gfs_dyn_matrix_sig_def
7 implicit none
8 real(kind=kind_evod) de_n(len_trie_ls,2,levs)
9
10 real(kind=kind_evod) de(len_trie_ls,2,levs),te(len_trie_ls,2,levs)
11 real(kind=kind_evod) xe(len_trie_ls,2,levs),ye(len_trie_ls,2,levs)
12 real(kind=kind_evod) ue(len_trie_ls,2,levs),ve(len_trie_ls,2,levs)
13 real(kind=kind_evod) qe(len_trie_ls,2), ze(len_trie_ls,2)
14 real(kind=kind_evod) am(levs,levs), bm(levs,levs)
15 real(kind=kind_evod) tov(levs), sv(levs)
16 real(kind=kind_evod) dt
17 integer ls_node(ls_dim,3)
18 real(kind=kind_evod) snnp1ev(len_trie_ls)
19 integer ndexev(len_trie_ls)
20 integer i,indev,indev1,indev2,j,k,l,locl,n
21 real(kind=kind_evod) qdtze(len_trie_ls,2),
22 . elne(len_trie_ls,2,levs)
23 real(kind=kind_evod) svdt, u1, u2
24 integer indlsev,jbasev
25 integer indlsod,jbasod
26 include 'function2'
27 real(kind=kind_evod) cons0,cons1,cons2
28
29 = 0.d0
30 = 1.d0
31 = 2.d0
32 = ls_node(locl,1)
33 jbasev = ls_node(locl,2)
34 indev1 = indlsev(L,L)
35 if (mod(L,2).eq.mod(jcap+1,2)) then
36 indev2 = indlsev(jcap+1,L)
37 else
38 indev2 = indlsev(jcap ,L)
39 endif
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58 do j=1,levs
59 do k=1,levs,2
60 do indev = indev1 , indev2
61 ye(indev,1,j) =
62 x ye(indev,1,j) - de_n(indev,1,k )*bm(j,k )
63 x - de_n(indev,1,k+1)*bm(j,k+1)
64 ye(indev,2,j) =
65 x ye(indev,2,j) - de_n(indev,2,k )*bm(j,k )
66 x - de_n(indev,2,k+1)*bm(j,k+1)
67 enddo
68 enddo
69 enddo
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88 do indev = indev1 , indev2
89 qdtze(indev,1) = qe(indev,1)
90 x + dt*ze(indev,1)
91 qdtze(indev,2) = qe(indev,2)
92 x + dt*ze(indev,2)
93 enddo
94 do k=1,levs
95 do indev = indev1 , indev2
96 elne(indev,1,k) = te(indev,1,k)
97 x + dt* ye(indev,1,k)
98 elne(indev,2,k) = te(indev,2,k)
99 x + dt* ye(indev,2,k)
100 enddo
101 enddo
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119 if ( kind_evod .eq. 8 ) then
120 call dgemm ('n', 't',
121 & indev2-indev1+1, levs, levs, cons1,
122 & elne(indev1,1,1), len_trie_ls*2,
123 & am(1,1), levs, cons0,
124 & ve(indev1,1,1), len_trie_ls*2)
125 call dgemm ('n', 't',
126 & indev2-indev1+1, levs, levs, cons1,
127 & elne(indev1,2,1), len_trie_ls*2,
128 & am(1,1), levs, cons0,
129 & ve(indev1,2,1), len_trie_ls*2)
130 else
131 call sgemm ('n', 't',
132 & indev2-indev1+1, levs, levs, cons1,
133 & elne(indev1,1,1), len_trie_ls*2,
134 & am(1,1), levs, cons0,
135 & ve(indev1,1,1), len_trie_ls*2)
136 call sgemm ('n', 't',
137 & indev2-indev1+1, levs, levs, cons1,
138 & elne(indev1,2,1), len_trie_ls*2,
139 & am(1,1), levs, cons0,
140 & ve(indev1,2,1), len_trie_ls*2)
141 endif
142
143 do 17 j=1,levs
144
145 do indev = indev1 , indev2
146 ve(indev,1,j) =
147 x ve(indev,1,j) + tor_sig(j)*qdtze(indev,1)
148 ve(indev,1,j) =
149 x ve(indev,1,j) * snnp1ev(indev)
150 ve(indev,1,j) =
151 x ve(indev,1,j) + xe(indev,1,j)
152 ue(indev,1,j) = de(indev,1,j)
153 x + ve(indev,1,j)*dt
154 ve(indev,2,j) =
155 x ve(indev,2,j) + tor_sig(j)*qdtze(indev,2)
156 ve(indev,2,j) =
157 x ve(indev,2,j) * snnp1ev(indev)
158 ve(indev,2,j) =
159 x ve(indev,2,j) + xe(indev,2,j)
160 ue(indev,2,j) = de(indev,2,j)
161 x + ve(indev,2,j)*dt
162 enddo
163 17 continue
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186 if ( kind_evod .eq. 8 ) then
187 do indev = indev1 , indev2
188 call dgemm ('n', 't',
189 & 1, levs, levs, cons1,
190 & ue(indev,1,1), len_trie_ls*2,
191 & D_m(1,1,ndexev(indev)+1), levs, cons0,
192 & ve(indev,1,1), len_trie_ls*2)
193 call dgemm ('n', 't',
194 & 1, levs, levs, cons1,
195 & ue(indev,2,1), len_trie_ls*2,
196 & D_m(1,1,ndexev(indev)+1), levs, cons0,
197 & ve(indev,2,1), len_trie_ls*2)
198 enddo
199 else
200 do indev = indev1 , indev2
201 call sgemm ('n', 't',
202 & 1, levs, levs, cons1,
203 & ue(indev,1,1), len_trie_ls*2,
204 & D_m(1,1,ndexev(indev)+1), levs, cons0,
205 & ve(indev,1,1), len_trie_ls*2)
206 call sgemm ('n', 't',
207 & 1, levs, levs, cons1,
208 & ue(indev,2,1), len_trie_ls*2,
209 & D_m(1,1,ndexev(indev)+1), levs, cons0,
210 & ve(indev,2,1), len_trie_ls*2)
211 enddo
212 endif
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227 if ( kind_evod .eq. 8 ) then
228 call dgemm ('n', 't',
229 & indev2-indev1+1, 1, levs, dt,
230 & ve(indev1,1,1), len_trie_ls*2,
231 & sv(1), 1, cons0,
232 & ue(indev1,1,1), len_trie_ls*2)
233 call dgemm ('n', 't',
234 & indev2-indev1+1, 1, levs, dt,
235 & ve(indev1,2,1), len_trie_ls*2,
236 & sv(1), 1, cons0,
237 & ue(indev1,2,1), len_trie_ls*2)
238 else
239 call sgemm ('n', 't',
240 & indev2-indev1+1, 1, levs, dt,
241 & ve(indev1,1,1), len_trie_ls*2,
242 & sv(1), 1, cons0,
243 & ue(indev1,1,1), len_trie_ls*2)
244 call sgemm ('n', 't',
245 & indev2-indev1+1, 1, levs, dt,
246 & ve(indev1,2,1), len_trie_ls*2,
247 & sv(1), 1, cons0,
248 & ue(indev1,2,1), len_trie_ls*2)
249 endif
250
251 do indev = indev1 , indev2
252 qdtze(indev,1) =
253 x qdtze(indev,1) + ue(indev,1,1)
254 ze(indev,1) = cons2*qdtze(indev,1)
255 - qe(indev,1)
256 qdtze(indev,2) =
257 x qdtze(indev,2) + ue(indev,2,1)
258 ze(indev,2) = cons2*qdtze(indev,2)
259 - qe(indev,2)
260 enddo
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278 if ( kind_evod .eq. 8 ) then
279 call dgemm ('n', 't',
280 & indev2-indev1+1, levs, levs, cons1,
281 & ve(indev1,1,1), len_trie_ls*2,
282 & bm(1,1), levs, cons0,
283 & ue(indev1,1,1), len_trie_ls*2)
284 call dgemm ('n', 't',
285 & indev2-indev1+1, levs, levs, cons1,
286 & ve(indev1,2,1), len_trie_ls*2,
287 & bm(1,1), levs, cons0,
288 & ue(indev1,2,1), len_trie_ls*2)
289 else
290 call sgemm ('n', 't',
291 & indev2-indev1+1, levs, levs, cons1,
292 & ve(indev1,1,1), len_trie_ls*2,
293 & bm(1,1), levs, cons0,
294 & ue(indev1,1,1), len_trie_ls*2)
295 call sgemm ('n', 't',
296 & indev2-indev1+1, levs, levs, cons1,
297 & ve(indev1,2,1), len_trie_ls*2,
298 & bm(1,1), levs, cons0,
299 & ue(indev1,2,1), len_trie_ls*2)
300 endif
301
302 do j=1,levs
303
304 do indev = indev1 , indev2
305 u1 = elne(indev,1,j) + dt * ue(indev,1,j)
306 ye(indev,1,j) = cons2*u1-te(indev,1,j)
307 (indev,1,j) = cons2* ve(indev,1,j)
308 -de(indev,1,j)
309 u2 = elne(indev,2,j) + dt * ue(indev,2,j)
310 ye(indev,2,j) = cons2*u2-te(indev,2,j)
311 (indev,2,j) = cons2* ve(indev,2,j)
312 -de(indev,2,j)
313 enddo
314 enddo
315 return
316 end
317 subroutine sicdifo_sig(do,to,qo,xo,yo,zo,am,bm,tov,sv,dt,uo,vo,
318 x ls_node,
319 x snnp1od,ndexod,locl,do_n)
320
321 use gfs_dyn_resol_def
322 use gfs_dyn_layout1
323 use gfs_dyn_matrix_sig_def
324 implicit none
325 real(kind=kind_evod) do_n(len_trio_ls,2,levs)
326
327 real(kind=kind_evod) do(len_trio_ls,2,levs),to(len_trio_ls,2,levs)
328 real(kind=kind_evod) xo(len_trio_ls,2,levs),yo(len_trio_ls,2,levs)
329 real(kind=kind_evod) uo(len_trio_ls,2,levs),vo(len_trio_ls,2,levs)
330 real(kind=kind_evod) qo(len_trio_ls,2), zo(len_trio_ls,2)
331 real(kind=kind_evod) am(levs,levs), bm(levs,levs)
332 real(kind=kind_evod) tov(levs), sv(levs)
333 real(kind=kind_evod) dt
334 integer ls_node(ls_dim,3)
335 real(kind=kind_evod) snnp1od(len_trio_ls)
336 integer ndexod(len_trio_ls)
337 integer i,indod,indod1,indod2,j,k,l,locl,n
338 real(kind=kind_evod) qdtzo(len_trio_ls,2),
339 . elno(len_trio_ls,2,levs)
340 real(kind=kind_evod) svdt, u1, u2
341 integer indlsev,jbasev
342 integer indlsod,jbasod
343 include 'function2'
344 real(kind=kind_evod) cons0,cons1,cons2
345 = 0.d0
346 = 1.d0
347 = 2.d0
348 = ls_node(locl,1)
349 jbasod = ls_node(locl,3)
350 indod1 = indlsod(L+1,L)
351 if (mod(L,2).eq.mod(jcap+1,2)) then
352 indod2 = indlsod(jcap ,L)
353 else
354 indod2 = indlsod(jcap+1,L)
355 endif
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374 do j=1,levs
375 do k=1,levs,2
376 do indod = indod1 , indod2
377 yo(indod,1,j) =
378 x yo(indod,1,j) - do_n(indod,1,k )*bm(j,k )
379 x - do_n(indod,1,k+1)*bm(j,k+1)
380 yo(indod,2,j) =
381 x yo(indod,2,j) - do_n(indod,2,k )*bm(j,k )
382 x - do_n(indod,2,k+1)*bm(j,k+1)
383 enddo
384 enddo
385 enddo
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404 do indod = indod1, indod2
405 qdtzo(indod,1) = qo(indod,1)
406 x + dt*zo(indod,1)
407 qdtzo(indod,2) = qo(indod,2)
408 x + dt*zo(indod,2)
409 enddo
410 do k=1,levs
411 do indod = indod1, indod2
412 elno(indod,1,k) = to(indod,1,k)
413 x + dt* yo(indod,1,k)
414
415 elno(indod,2,k) = to(indod,2,k)
416 x + dt* yo(indod,2,k)
417 enddo
418 enddo
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436 if ( kind_evod .eq. 8 ) then
437 call dgemm ('n', 't',
438 & indod2-indod1+1, levs, levs, cons1,
439 & elno(indod1,1,1), len_trio_ls*2,
440 & am(1,1), levs, cons0,
441 & vo(indod1,1,1), len_trio_ls*2)
442 call dgemm ('n', 't',
443 & indod2-indod1+1, levs, levs, cons1,
444 & elno(indod1,2,1), len_trio_ls*2,
445 & am(1,1), levs, cons0,
446 & vo(indod1,2,1), len_trio_ls*2)
447 else
448 call sgemm ('n', 't',
449 & indod2-indod1+1, levs, levs, cons1,
450 & elno(indod1,1,1), len_trio_ls*2,
451 & am(1,1), levs, cons0,
452 & vo(indod1,1,1), len_trio_ls*2)
453 call sgemm ('n', 't',
454 & indod2-indod1+1, levs, levs, cons1,
455 & elno(indod1,2,1), len_trio_ls*2,
456 & am(1,1), levs, cons0,
457 & vo(indod1,2,1), len_trio_ls*2)
458 endif
459
460 do 17 j=1,levs
461
462 do indod = indod1, indod2
463 vo(indod,1,j) =
464 x vo(indod,1,j) + tor_sig(j)*qdtzo(indod,1)
465
466 vo(indod,1,j) =
467 x vo(indod,1,j) * snnp1od(indod)
468
469 vo(indod,1,j) =
470 x vo(indod,1,j) + xo(indod,1,j)
471
472 uo(indod,1,j) = do(indod,1,j)
473 x + vo(indod,1,j)*dt
474 vo(indod,2,j) =
475 x vo(indod,2,j) + tor_sig(j)*qdtzo(indod,2)
476 vo(indod,2,j) =
477 x vo(indod,2,j) * snnp1od(indod)
478 vo(indod,2,j) =
479 x vo(indod,2,j) + xo(indod,2,j)
480 uo(indod,2,j) = do(indod,2,j)
481 x + vo(indod,2,j)*dt
482 enddo
483 17 continue
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506 if ( kind_evod .eq. 8 ) then
507 do indod = indod1 , indod2
508 call dgemm ('n', 't',
509 & 1, levs, levs, cons1,
510 & uo(indod,1,1), len_trio_ls*2,
511 & D_m(1,1,ndexod(indod)+1), levs, cons0,
512 & vo(indod,1,1), len_trio_ls*2)
513 call dgemm ('n', 't',
514 & 1, levs, levs, cons1,
515 & uo(indod,2,1), len_trio_ls*2,
516 & D_m(1,1,ndexod(indod)+1), levs, cons0,
517 & vo(indod,2,1), len_trio_ls*2)
518 enddo
519 else
520 do indod = indod1 , indod2
521 call sgemm ('n', 't',
522 & 1, levs, levs, cons1,
523 & uo(indod,1,1), len_trio_ls*2,
524 & D_m(1,1,ndexod(indod)+1), levs, cons0,
525 & vo(indod,1,1), len_trio_ls*2)
526 call sgemm ('n', 't',
527 & 1, levs, levs, cons1,
528 & uo(indod,2,1), len_trio_ls*2,
529 & D_m(1,1,ndexod(indod)+1), levs, cons0,
530 & vo(indod,2,1), len_trio_ls*2)
531 enddo
532 endif
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547 if ( kind_evod .eq. 8 ) then
548 call dgemm ('n', 't',
549 & indod2-indod1+1, 1, levs, dt,
550 & vo(indod1,1,1), len_trio_ls*2,
551 & sv(1), 1, cons0,
552 & uo(indod1,1,1), len_trio_ls*2)
553 call dgemm ('n', 't',
554 & indod2-indod1+1, 1, levs, dt,
555 & vo(indod1,2,1), len_trio_ls*2,
556 & sv(1), 1, cons0,
557 & uo(indod1,2,1), len_trio_ls*2)
558 else
559 call sgemm ('n', 't',
560 & indod2-indod1+1, 1, levs, dt,
561 & vo(indod1,1,1), len_trio_ls*2,
562 & sv(1), 1, cons0,
563 & uo(indod1,1,1), len_trio_ls*2)
564 call sgemm ('n', 't',
565 & indod2-indod1+1, 1, levs, dt,
566 & vo(indod1,2,1), len_trio_ls*2,
567 & sv(1), 1, cons0,
568 & uo(indod1,2,1), len_trio_ls*2)
569 endif
570
571 do indod = indod1, indod2
572 qdtzo(indod,1) =
573 x qdtzo(indod,1) + uo(indod,1,1)
574 zo(indod,1) = cons2*qdtzo(indod,1)
575 - qo(indod,1)
576 qdtzo(indod,2) =
577 x qdtzo(indod,2) + uo(indod,2,1)
578 zo(indod,2) = cons2*qdtzo(indod,2)
579 - qo(indod,2)
580 enddo
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598 if ( kind_evod .eq. 8 ) then
599 call dgemm ('n', 't',
600 & indod2-indod1+1, levs, levs, cons1,
601 & vo(indod1,1,1), len_trio_ls*2,
602 & bm(1,1), levs, cons0,
603 & uo(indod1,1,1), len_trio_ls*2)
604 call dgemm ('n', 't',
605 & indod2-indod1+1, levs, levs, cons1,
606 & vo(indod1,2,1), len_trio_ls*2,
607 & bm(1,1), levs, cons0,
608 & uo(indod1,2,1), len_trio_ls*2)
609 else
610 call sgemm ('n', 't',
611 & indod2-indod1+1, levs, levs, cons1,
612 & vo(indod1,1,1), len_trio_ls*2,
613 & bm(1,1), levs, cons0,
614 & uo(indod1,1,1), len_trio_ls*2)
615 call sgemm ('n', 't',
616 & indod2-indod1+1, levs, levs, cons1,
617 & vo(indod1,2,1), len_trio_ls*2,
618 & bm(1,1), levs, cons0,
619 & uo(indod1,2,1), len_trio_ls*2)
620 endif
621
622 do j=1,levs
623
624 do indod = indod1, indod2
625 u1 = elno(indod,1,j) + dt * uo(indod,1,j)
626 yo(indod,1,j) = cons2*u1-to(indod,1,j)
627 (indod,1,j) = cons2* vo(indod,1,j)
628 -do(indod,1,j)
629 u2 = elno(indod,2,j) + dt * uo(indod,2,j)
630 yo(indod,2,j) = cons2*u2-to(indod,2,j)
631 (indod,2,j) = cons2* vo(indod,2,j)
632 -do(indod,2,j)
633 enddo
634 enddo
635 return
636 end
637 subroutine matinv(a,m,n,d,p,r)
638 use gfs_dyn_machine
639 implicit none
640 integer i,j,k,l,m,n
641 real(kind=kind_evod) a(m,n,n),d(m),p(m),r(m)
642 real(kind=kind_evod) cons0,cons1
643 = 0.d0
644 = 1.d0
645 do 200 l=1,m
646 d(l)=cons1
647 continue
648 do 100 k=1,n
649 do 250 l=1,m
650 p(l)=a(l,k,k)
651 250 continue
652 do 300 l=1,m
653 r(l)=-cons1/p(l)
654 continue
655 do 350 l=1,m
656 a(l,k,k)=cons0
657 continue
658 do 20 i=1,n
659 do 400 l=1,m
660 a(l,i,k)=a(l,i,k)*r(l)
661 400 continue
662 20 continue
663 do 60 i=1,n
664 if(i.eq.k) go to 60
665 do 40 j=1,n
666 do 450 l=1,m
667 a(l,i,j)=a(l,i,k)*a(l,k,j)+a(l,i,j)
668 450 continue
669 40 continue
670 60 continue
671 do 600 l=1,m
672 r(l)=-r(l)
673 600 continue
674 do 80 j=1,n
675 do 650 l=1,m
676 a(l,k,j)=a(l,k,j)*r(l)
677 650 continue
678 80 continue
679 do 700 l=1,m
680 d(l)=d(l)*p(l)
681 700 continue
682 do 750 l=1,m
683 a(l,k,k)=r(l)
684 750 continue
685 100 continue
686 return
687 end
688