Changeset 12523


Ignore:
Timestamp:
2020-03-09T11:59:47+01:00 (8 months ago)
Author:
clem
Message:

solve ticket #2398. This routine is now much much faster

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/release-4.0-HEAD/src/OCE/DIA/diaharm.F90

    r11536 r12523  
    2424    
    2525   INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo 
    26    INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24 
     26   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*366*24*2 ! 30min for a 1yr-long run 
    2727 
    2828   !                         !!** namelist variables ** 
     
    3535   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
    3636   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp 
    37    REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut   , vt   , ft 
    38    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta , out_u, out_v 
    39  
    40    INTEGER ::   ninco, nsparse 
     37   REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut, vt, ft 
     38   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta, out_u, out_v 
     39 
     40   INTEGER  ::   ninco, nsparse 
     41   REAL(wp) ::   z1_tmp3 
    4142   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse 
    4243   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1 
    4344   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse 
    44    REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7 
     45   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7, z1_pivot 
    4546   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier 
    46    REAL(wp), SAVE, DIMENSION(jpincomax)           ::   zpivot 
    4747 
    4848   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...) 
     
    7676         WRITE(numout,*) 
    7777         WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization' 
    78          WRITE(numout,*) '~~~~~~~ ' 
     78         WRITE(numout,*) '~~~~~~~~~~~~~ ' 
    7979      ENDIF 
    8080      ! 
     
    131131         ENDIF 
    132132 
    133          ALLOCATE(name    (nb_ana)) 
     133         ALLOCATE(name(nb_ana)) 
    134134         DO jh=1,nb_ana 
    135135            DO ji=1,jpmax_harmo 
     
    163163 
    164164 
    165    SUBROUTINE dia_harm ( kt ) 
     165   SUBROUTINE dia_harm( kt ) 
    166166      !!---------------------------------------------------------------------- 
    167167      !!                 ***  ROUTINE dia_harm  *** 
     
    172172      !! 
    173173      !!-------------------------------------------------------------------- 
    174       INTEGER, INTENT( IN ) :: kt 
    175       ! 
    176       INTEGER  :: ji, jj, jh, jc, nhc 
    177       REAL(wp) :: ztime, ztemp 
     174      INTEGER, INTENT( in ) ::  kt 
     175      ! 
     176      INTEGER  ::   ji, jj, jh, jc, nhc 
     177      REAL(wp) ::   ztime, ztemp 
    178178      !!-------------------------------------------------------------------- 
    179179      IF( ln_timing )   CALL timing_start('dia_harm') 
     
    187187            DO jc = 1, 2 
    188188               nhc = nhc+1 
    189                ztemp =   MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  & 
     189               ztemp =  (   MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  & 
    190190                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
    191191                  ! 
    192                DO jj = 1,jpj 
    193                   DO ji = 1,jpi 
    194                      ! Elevation 
    195                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)         
    196                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
    197                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
     192               DO jj = 2, jpjm1 
     193                  DO ji = 2, jpim1 
     194                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp * sshn(ji,jj) * ssmask (ji,jj) ! elevation       
     195                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp * un_b(ji,jj) * ssumask(ji,jj) ! u-vel 
     196                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp * vn_b(ji,jj) * ssvmask(ji,jj) ! v-vel 
    198197                  END DO 
    199198               END DO 
    200                ! 
    201             END DO 
    202          END DO 
    203          !        
     199            END DO 
     200         END DO 
    204201      END IF 
    205202      ! 
     
    220217      !! 
    221218      !!-------------------------------------------------------------------- 
    222       INTEGER :: ji, jj, jh, jc, jn, nhan, jl 
    223       INTEGER :: ksp, kun, keq 
    224       REAL(wp) :: ztime, ztime_ini, ztime_end 
    225       REAL(wp) :: X1, X2 
    226       REAL(wp), DIMENSION(jpi,jpj,jpmax_harmo,2) ::   ana_amp   ! workspace 
     219      INTEGER  ::   ji, jj, jh, jc, jn, nhan 
     220      INTEGER  ::   ksp, kun, keq 
     221      REAL(wp) ::   ztime, ztime_ini, ztime_end, z1_han 
    227222      !!-------------------------------------------------------------------- 
    228223      ! 
    229224      IF(lwp) WRITE(numout,*) 
    230       IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis' 
     225      IF(lwp) WRITE(numout,*) 'dia_harm_end: kt=nitend_han: Perform harmonic analysis' 
    231226      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     227       
     228      ALLOCATE( out_eta(jpi,jpj,2*nb_ana), out_u(jpi,jpj,2*nb_ana), out_v(jpi,jpj,2*nb_ana) ) 
    232229 
    233230      ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis 
    234231      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis 
    235232      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis 
    236  
     233      z1_han = 1._wp / REAL(nhan-1)  
     234       
    237235      ninco = 2*nb_ana 
    238236 
     
    240238      keq = 0         
    241239      DO jn = 1, nhan 
    242          ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1) 
     240         ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end ) * z1_han 
    243241         keq = keq + 1 
    244242         kun = 0 
     
    257255      nsparse = ksp 
    258256 
     257      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
     258      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
     259 
     260      CALL SUR_DETERMINE_INIT 
     261 
    259262      ! Elevation: 
    260       DO jj = 1, jpj 
    261          DO ji = 1, jpi 
     263      DO jj = 2, jpjm1 
     264         DO ji = 2, jpim1 
     265 
    262266            ! Fill input array 
    263             kun = 0 
    264             DO jh = 1, nb_ana 
    265                DO jc = 1, 2 
    266                   kun = kun + 1 
    267                   ztmp4(kun)=ana_temp(ji,jj,kun,1) 
    268                END DO 
    269             END DO 
    270  
    271             CALL SUR_DETERMINE(jj) 
    272  
     267            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1) 
     268            CALL SUR_DETERMINE 
     269             
    273270            ! Fill output array 
    274271            DO jh = 1, nb_ana 
    275                ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
    276                ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
    277             END DO 
    278          END DO 
    279       END DO 
    280  
    281       ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   &  
    282          &      out_u  (jpi,jpj,2*nb_ana),   & 
    283          &      out_v  (jpi,jpj,2*nb_ana)  ) 
    284  
    285       DO jj = 1, jpj 
    286          DO ji = 1, jpi 
    287             DO jh = 1, nb_ana  
    288                X1 = ana_amp(ji,jj,jh,1) 
    289                X2 =-ana_amp(ji,jj,jh,2) 
    290                out_eta(ji,jj,jh       ) = X1 * tmask_i(ji,jj) 
    291                out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj) 
     272               out_eta(ji,jj,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj) 
     273               out_eta(ji,jj,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 
    292274            END DO 
    293275         END DO 
     
    295277 
    296278      ! ubar: 
    297       DO jj = 1, jpj 
    298          DO ji = 1, jpi 
     279      DO jj = 2, jpjm1 
     280         DO ji = 2, jpim1 
     281 
    299282            ! Fill input array 
    300             kun=0 
    301             DO jh = 1,nb_ana 
    302                DO jc = 1,2 
    303                   kun = kun + 1 
    304                   ztmp4(kun)=ana_temp(ji,jj,kun,2) 
    305                END DO 
    306             END DO 
    307  
    308             CALL SUR_DETERMINE(jj+1) 
     283            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2) 
     284            CALL SUR_DETERMINE 
    309285 
    310286            ! Fill output array 
    311287            DO jh = 1, nb_ana 
    312                ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1) 
    313                ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2) 
    314             END DO 
    315  
    316          END DO 
    317       END DO 
    318  
    319       DO jj = 1, jpj 
    320          DO ji = 1, jpi 
    321             DO jh = 1, nb_ana  
    322                X1= ana_amp(ji,jj,jh,1) 
    323                X2=-ana_amp(ji,jj,jh,2) 
    324                out_u(ji,jj,       jh) = X1 * ssumask(ji,jj) 
    325                out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 
    326             ENDDO 
    327          ENDDO 
    328       ENDDO 
     288               out_u(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj) 
     289               out_u(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 
     290            END DO 
     291 
     292        END DO 
     293      END DO 
    329294 
    330295      ! vbar: 
    331       DO jj = 1, jpj 
    332          DO ji = 1, jpi 
     296      DO jj = 2, jpjm1 
     297         DO ji = 2, jpim1 
     298 
    333299            ! Fill input array 
    334             kun=0 
    335             DO jh = 1,nb_ana 
    336                DO jc = 1,2 
    337                   kun = kun + 1 
    338                   ztmp4(kun)=ana_temp(ji,jj,kun,3) 
    339                END DO 
    340             END DO 
    341  
    342             CALL SUR_DETERMINE(jj+1) 
     300            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3) 
     301            CALL SUR_DETERMINE 
    343302 
    344303            ! Fill output array 
    345304            DO jh = 1, nb_ana 
    346                ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
    347                ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
    348             END DO 
    349  
    350          END DO 
    351       END DO 
    352  
    353       DO jj = 1, jpj 
    354          DO ji = 1, jpi 
    355             DO jh = 1, nb_ana  
    356                X1=ana_amp(ji,jj,jh,1) 
    357                X2=-ana_amp(ji,jj,jh,2) 
    358                out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj) 
    359                out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 
    360             END DO 
    361          END DO 
    362       END DO 
     305               out_v(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 
     306               out_v(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 
     307            END DO 
     308 
     309         END DO 
     310      END DO 
     311      ! 
     312      ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj 
     313      !       but I think this is the most efficient 
     314      CALL lbc_lnk_multi( 'dia_harm_end', out_eta, 'T', 1., out_u, 'U', -1. , out_v, 'V', -1. ) 
    363315      ! 
    364316      CALL dia_wri_harm ! Write results in files 
     317      ! 
     318      DEALLOCATE( out_eta, out_u, out_v ) 
    365319      ! 
    366320   END SUBROUTINE dia_harm_end 
     
    373327      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file 
    374328      !!-------------------------------------------------------------------- 
    375       CHARACTER(LEN=lc) :: cltext 
    376       CHARACTER(LEN=lc) ::   & 
    377          cdfile_name_T   ,   & ! name of the file created (T-points) 
    378          cdfile_name_U   ,   & ! name of the file created (U-points) 
    379          cdfile_name_V         ! name of the file created (V-points) 
    380329      INTEGER  ::   jh 
    381330      !!---------------------------------------------------------------------- 
     
    383332      IF(lwp) WRITE(numout,*) '  ' 
    384333      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 
    385       IF(lwp) WRITE(numout,*) '  ' 
     334      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    386335 
    387336      ! A) Elevation 
    388337      !///////////// 
    389       ! 
    390338      DO jh = 1, nb_ana 
    391339      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) 
    392       CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) ) 
     340      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) ) 
    393341      END DO 
    394342 
    395343      ! B) ubar 
    396344      !///////// 
    397       ! 
    398345      DO jh = 1, nb_ana 
    399346      CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) 
    400       CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) ) 
     347      CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,jh+nb_ana) ) 
    401348      END DO 
    402349 
    403350      ! C) vbar 
    404351      !///////// 
    405       ! 
    406352      DO jh = 1, nb_ana 
    407353         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) ) 
     
    412358 
    413359 
    414    SUBROUTINE SUR_DETERMINE(init) 
     360   SUBROUTINE SUR_DETERMINE_INIT 
     361      !!--------------------------------------------------------------------------------- 
     362      !!                      *** ROUTINE SUR_DETERMINE_INIT *** 
     363      !!        
     364      !!--------------------------------------------------------------------------------- 
     365      INTEGER                        :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd 
     366      INTEGER                        :: ipivot 
     367      REAL(wp)                       :: zval1, zval2, zcol1, zcol2 
     368      INTEGER , DIMENSION(jpincomax) :: ipos2 
     369      !!--------------------------------------------------------------------------------- 
     370      !             
     371      ! 
     372      ztmp3(:,:) = 0._wp 
     373      ! 
     374      DO jh1_sd = 1, nsparse 
     375         DO jh2_sd = 1, nsparse 
     376            IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN 
     377               ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  & 
     378                  &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd) 
     379            ENDIF 
     380         END DO 
     381      END DO 
     382      ! 
     383      DO jj_sd = 1, ninco 
     384         ipos1(jj_sd) = jj_sd 
     385         ipos2(jj_sd) = jj_sd 
     386      END DO 
     387      ! 
     388      DO ji_sd = 1, ninco 
     389         ! 
     390         !find greatest non-zero pivot: 
     391         zval1 = ABS(ztmp3(ji_sd,ji_sd)) 
     392         ! 
     393         ipivot = ji_sd 
     394         DO jj_sd = ji_sd, ninco 
     395            zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
     396            IF( zval2 >= zval1 )THEN 
     397               ipivot = jj_sd 
     398               zval1  = zval2 
     399            ENDIF 
     400         END DO 
     401         ! 
     402         DO ji1_sd = 1, ninco 
     403            zcol1                = ztmp3(ji1_sd,ji_sd) 
     404            zcol2                = ztmp3(ji1_sd,ipivot) 
     405            ztmp3(ji1_sd,ji_sd)  = zcol2 
     406            ztmp3(ji1_sd,ipivot) = zcol1 
     407         END DO 
     408         ! 
     409         ipos2(ji_sd)    = ipos1(ipivot) 
     410         ipos2(ipivot)   = ipos1(ji_sd) 
     411         ipos1(ji_sd)    = ipos2(ji_sd) 
     412         ipos1(ipivot)   = ipos2(ipivot) 
     413         z1_pivot(ji_sd) = 1._wp / ztmp3(ji_sd,ji_sd) 
     414         DO jj_sd = 1, ninco 
     415            ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) * z1_pivot(ji_sd) 
     416         END DO 
     417         ! 
     418         DO ji2_sd = ji_sd+1, ninco 
     419            zpilier(ji2_sd,ji_sd) = ztmp3(ji2_sd,ji_sd) 
     420            DO jj_sd=1,ninco 
     421               ztmp3(ji2_sd,jj_sd) = ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 
     422            END DO 
     423         END DO 
     424         ! 
     425      END DO 
     426      ! 
     427      z1_tmp3 = 1._wp / ztmp3(ninco,ninco) 
     428      ! 
     429   END SUBROUTINE SUR_DETERMINE_INIT 
     430 
     431    
     432   SUBROUTINE SUR_DETERMINE 
    415433      !!--------------------------------------------------------------------------------- 
    416434      !!                      *** ROUTINE SUR_DETERMINE *** 
    417435      !!     
    418       !!     
    419       !!        
    420       !!--------------------------------------------------------------------------------- 
    421       INTEGER, INTENT(in) ::   init  
    422       ! 
    423       INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd 
    424       REAL(wp)                        :: zval1, zval2, zx1 
    425       REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2 
    426       INTEGER , DIMENSION(jpincomax) :: ipos2, ipivot 
    427       !--------------------------------------------------------------------------------- 
     436      !!--------------------------------------------------------------------------------- 
     437      INTEGER                        :: ji_sd, jj_sd, ji1_sd, ji2_sd 
     438      REAL(wp)                       :: zx1 
     439      REAL(wp), DIMENSION(jpincomax) :: ztmpx 
     440      !!--------------------------------------------------------------------------------- 
    428441      !             
    429       IF( init == 1 ) THEN 
    430          IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
    431          IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
    432          ! 
    433          ztmp3(:,:) = 0._wp 
    434          ! 
    435          DO jh1_sd = 1, nsparse 
    436             DO jh2_sd = 1, nsparse 
    437                nisparse(jh2_sd) = nisparse(jh2_sd) 
    438                njsparse(jh2_sd) = njsparse(jh2_sd) 
    439                IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN 
    440                   ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  & 
    441                      &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd) 
    442                ENDIF 
    443             END DO 
    444          END DO 
    445          ! 
    446          DO jj_sd = 1 ,ninco 
    447             ipos1(jj_sd) = jj_sd 
    448             ipos2(jj_sd) = jj_sd 
    449          END DO 
    450          ! 
    451          DO ji_sd = 1 , ninco 
    452             ! 
    453             !find greatest non-zero pivot: 
    454             zval1 = ABS(ztmp3(ji_sd,ji_sd)) 
    455             ! 
    456             ipivot(ji_sd) = ji_sd 
    457             DO jj_sd = ji_sd, ninco 
    458                zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
    459                IF( zval2 >= zval1 )THEN 
    460                   ipivot(ji_sd) = jj_sd 
    461                   zval1         = zval2 
    462                ENDIF 
    463             END DO 
    464             ! 
    465             DO ji1_sd = 1, ninco 
    466                zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd) 
    467                zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd)) 
    468                ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd) 
    469                ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 
    470             END DO 
    471             ! 
    472             ipos2(ji_sd)         = ipos1(ipivot(ji_sd)) 
    473             ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 
    474             ipos1(ji_sd)         = ipos2(ji_sd) 
    475             ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 
    476             zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd) 
    477             DO jj_sd = 1, ninco 
    478                ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 
    479             END DO 
    480             ! 
    481             DO ji2_sd = ji_sd+1, ninco 
    482                zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 
    483                DO jj_sd=1,ninco 
    484                   ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 
    485                END DO 
    486             END DO 
    487             ! 
    488          END DO 
    489          ! 
    490       ENDIF ! End init==1 
    491  
    492442      DO ji_sd = 1, ninco 
    493          ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 
     443         ztmp4(ji_sd) = ztmp4(ji_sd) * z1_pivot(ji_sd) 
    494444         DO ji2_sd = ji_sd+1, ninco 
    495445            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 
     
    498448 
    499449      !system solving:  
    500       ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 
    501       ji_sd = ninco 
     450      ztmpx(ninco) = ztmp4(ninco) * z1_tmp3 
    502451      DO ji_sd = ninco-1, 1, -1 
    503452         zx1 = 0._wp 
     
    505454            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 
    506455         END DO 
    507          ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 
    508       END DO 
    509  
    510       DO jj_sd =1, ninco 
    511          ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 
     456         ztmpx(ji_sd) = ztmp4(ji_sd) - zx1 
     457      END DO 
     458 
     459      DO jj_sd = 1, ninco 
     460         ztmp7(ipos1(jj_sd)) = ztmpx(jj_sd) 
    512461      END DO 
    513462      ! 
    514463   END SUBROUTINE SUR_DETERMINE 
    515464 
     465    
    516466   !!====================================================================== 
    517467END MODULE diaharm 
Note: See TracChangeset for help on using the changeset viewer.