File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\sicdif_hyb_gc.f

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