Changeset 4683


Ignore:
Timestamp:
2014-06-21T13:23:16+02:00 (7 years ago)
Author:
gm
Message:

#1233 : correct a bug in diaharm.F90 in the v3.6alpha (trunk)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r4624 r4683  
    1818   USE daymod 
    1919   USE tide_mod 
     20   ! 
    2021   USE in_out_manager  ! I/O units 
    2122   USE iom             ! I/0 library 
     
    3435   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24 
    3536 
    36    !                            !!!namelist variables 
     37   !                         !!** namelist variables ** 
    3738   INTEGER ::   nit000_han    ! First time step used for harmonic analysis 
    3839   INTEGER ::   nitend_han    ! Last time step used for harmonic analysis 
    3940   INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis 
    40    INTEGER ::   nb_ana           ! Number of harmonics to analyse 
     41   INTEGER ::   nb_ana        ! Number of harmonics to analyse 
    4142 
    4243   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
     
    119120            ENDIF 
    120121         END DO 
    121       ENDDO 
     122      END DO 
    122123      ! 
    123124      IF(lwp) THEN 
     
    158159      ! ---------------------------- 
    159160      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    160       ana_temp(:,:,:,:) = 0.e0 
     161      ana_temp(:,:,:,:) = 0._wp 
    161162 
    162163   END SUBROUTINE dia_harm_init 
     
    179180      IF( nn_timing == 1 )   CALL timing_start('dia_harm') 
    180181 
    181       IF ( kt == nit000 ) CALL dia_harm_init 
    182  
    183       IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & 
    184            (MOD(kt,nstep_han).EQ.0) ) THEN 
    185  
    186         ztime = (kt-nit000+1)*rdt  
     182      IF( kt == nit000 ) CALL dia_harm_init 
     183 
     184      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 
     185 
     186         ztime = (kt-nit000+1) * rdt  
    187187        
    188         nhc = 0 
    189         DO jh = 1,nb_ana 
    190           DO jc = 1,2 
    191             nhc = nhc+1 
    192             ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  & 
    193                     +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
    194  
    195             DO jj = 1,jpj 
    196               DO ji = 1,jpi 
    197                 ! Elevation 
    198                 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)         
     188         nhc = 0 
     189         DO jh = 1, nb_ana 
     190            DO jc = 1, 2 
     191               nhc = nhc+1 
     192               ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  & 
     193                    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
     194 
     195               DO jj = 1, jpj 
     196                  DO ji = 1, jpi 
     197                     ! Elevation 
     198                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1) 
    199199#if defined key_dynspg_ts 
    200                 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 
    201                 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 
    202 #endif 
    203               END DO 
    204             END DO 
    205  
    206           END DO 
    207         END DO 
    208         
     200                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 
     201                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 
     202#endif 
     203                  END DO 
     204               END DO 
     205               ! 
     206            END DO 
     207         END DO 
     208         !        
    209209      END IF 
    210210 
     
    249249         keq = keq + 1 
    250250         kun = 0 
    251          DO jh = 1,nb_ana 
    252             DO jc = 1,2 
     251         DO jh = 1, nb_ana 
     252            DO jc = 1, 2 
    253253               kun = kun + 1 
    254254               ksp = ksp + 1 
     
    296296               out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1) 
    297297               out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1) 
    298             ENDDO 
    299          ENDDO 
    300       ENDDO 
     298            END DO 
     299         END DO 
     300      END DO 
    301301 
    302302      ! ubar: 
     
    309309                  kun = kun + 1 
    310310                  ztmp4(kun)=ana_temp(ji,jj,kun,2) 
    311                ENDDO 
    312             ENDDO 
     311               END DO 
     312            END DO 
    313313 
    314314            CALL SUR_DETERMINE(jj+1) 
     
    316316            ! Fill output array 
    317317            DO jh = 1, nb_ana 
    318                ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
    319                ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
     318               ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1) 
     319               ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2) 
    320320            END DO 
    321321 
     
    326326         DO ji = 1, jpi 
    327327            DO jh = 1, nb_ana  
    328                X1=ana_amp(ji,jj,jh,1) 
    329                X2=-ana_amp(ji,jj,jh,2) 
    330                out_u(ji,jj,jh) = X1 * umask(ji,jj,1) 
    331                out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1) 
    332             ENDDO 
    333          ENDDO 
    334       ENDDO 
     328               X1 = ana_amp(ji,jj,jh,1) 
     329               X2 =-ana_amp(ji,jj,jh,2) 
     330               out_u(ji,jj,jh       ) = X1 * umask(ji,jj,1) 
     331               out_u(ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1) 
     332            END DO 
     333         END DO 
     334      END DO 
    335335 
    336336      ! vbar: 
     
    343343                  kun = kun + 1 
    344344                  ztmp4(kun)=ana_temp(ji,jj,kun,3) 
    345                ENDDO 
    346             ENDDO 
     345               END DO 
     346            END DO 
    347347 
    348348            CALL SUR_DETERMINE(jj+1) 
     
    364364               out_v(ji,jj,jh)=X1 * vmask(ji,jj,1) 
    365365               out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1) 
    366             ENDDO 
    367          ENDDO 
    368       ENDDO 
     366            END DO 
     367         END DO 
     368      END DO 
    369369 
    370370      CALL dia_wri_harm ! Write results in files 
     
    437437#else 
    438438      DO jh = 1, nb_ana 
    439          CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh       ) ) 
    440          CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,jh+nb_ana) ) 
    441       END DO 
    442 #endif 
    443  
     439         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) ) 
     440         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 
     441      END DO 
     442#endif 
     443      ! 
    444444   END SUBROUTINE dia_wri_harm 
    445445 
    446446 
    447447   SUBROUTINE SUR_DETERMINE(init) 
    448    !!--------------------------------------------------------------------------------- 
    449    !!                      *** ROUTINE SUR_DETERMINE *** 
    450    !!     
    451    !!     
    452    !!        
    453    !!--------------------------------------------------------------------------------- 
    454    INTEGER, INTENT(in) ::   init  
    455    ! 
    456    INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
    457    REAL(wp)                        :: zval1, zval2, zx1 
    458    REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2 
    459    INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot 
    460    !--------------------------------------------------------------------------------- 
    461    CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
    462    CALL wrk_alloc( jpincomax , ipos2 , ipivot        ) 
     448      !!--------------------------------------------------------------------------------- 
     449      !!                      *** ROUTINE SUR_DETERMINE *** 
     450      !!     
     451      !!     
     452      !!        
     453      !!--------------------------------------------------------------------------------- 
     454      INTEGER, INTENT(in) ::   init  
     455      ! 
     456      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
     457      REAL(wp)                        :: zval1, zval2, zx1 
     458      REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2 
     459      INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot 
     460      !--------------------------------------------------------------------------------- 
     461      CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
     462      CALL wrk_alloc( jpincomax , ipos2 , ipivot        ) 
    463463             
    464    IF( init == 1 ) THEN 
    465       IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
    466       IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
    467       ! 
    468       ztmp3(:,:) = 0._wp 
    469       ! 
    470       DO jk1_sd = 1, nsparse 
    471          DO jk2_sd = 1, nsparse 
    472             nisparse(jk2_sd) = nisparse(jk2_sd) 
    473             njsparse(jk2_sd) = njsparse(jk2_sd) 
    474             IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
    475                ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
    476                                                         + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
    477             ENDIF 
    478          END DO 
    479       END DO 
    480  
    481       DO jj_sd = 1 ,ninco 
    482           ipos1(jj_sd) = jj_sd 
    483           ipos2(jj_sd) = jj_sd 
    484       ENDDO 
    485  
    486       DO ji_sd = 1 , ninco 
    487  
    488          !find greatest non-zero pivot: 
    489          zval1 = ABS(ztmp3(ji_sd,ji_sd)) 
    490  
    491          ipivot(ji_sd) = ji_sd 
    492          DO jj_sd = ji_sd, ninco 
    493             zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
    494             IF( zval2.GE.zval1 )THEN 
    495                ipivot(ji_sd) = jj_sd 
    496                zval1         = zval2 
    497             ENDIF 
    498          ENDDO 
    499  
    500          DO ji1_sd = 1, ninco 
    501             zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd) 
    502             zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd)) 
    503             ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd) 
    504             ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 
    505          ENDDO 
    506  
    507          ipos2(ji_sd)         = ipos1(ipivot(ji_sd)) 
    508          ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 
    509          ipos1(ji_sd)         = ipos2(ji_sd) 
    510          ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 
    511          zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd) 
    512          DO jj_sd = 1, ninco 
    513             ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 
    514          ENDDO 
    515  
     464      IF( init == 1 ) THEN 
     465         IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
     466         IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
     467         ! 
     468         ztmp3(:,:) = 0._wp 
     469         ! 
     470         DO jk1_sd = 1, nsparse 
     471            DO jk2_sd = 1, nsparse 
     472               nisparse(jk2_sd) = nisparse(jk2_sd) 
     473               njsparse(jk2_sd) = njsparse(jk2_sd) 
     474               IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
     475                  ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
     476                     &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
     477               ENDIF 
     478            END DO 
     479         END DO 
     480         ! 
     481         DO jj_sd = 1 ,ninco 
     482            ipos1(jj_sd) = jj_sd 
     483            ipos2(jj_sd) = jj_sd 
     484         END DO 
     485         ! 
     486         DO ji_sd = 1 , ninco 
     487            ! 
     488            !find greatest non-zero pivot: 
     489            zval1 = ABS(ztmp3(ji_sd,ji_sd)) 
     490            ! 
     491            ipivot(ji_sd) = ji_sd 
     492            DO jj_sd = ji_sd, ninco 
     493               zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
     494               IF( zval2.GE.zval1 )THEN 
     495                  ipivot(ji_sd) = jj_sd 
     496                  zval1         = zval2 
     497               ENDIF 
     498            END DO 
     499            ! 
     500            DO ji1_sd = 1, ninco 
     501               zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd) 
     502               zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd)) 
     503               ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd) 
     504               ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 
     505            END DO 
     506            ! 
     507            ipos2(ji_sd)         = ipos1(ipivot(ji_sd)) 
     508            ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 
     509            ipos1(ji_sd)         = ipos2(ji_sd) 
     510            ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 
     511            zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd) 
     512            DO jj_sd = 1, ninco 
     513               ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 
     514            END DO 
     515            ! 
     516            DO ji2_sd = ji_sd+1, ninco 
     517               zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 
     518               DO jj_sd=1,ninco 
     519                  ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 
     520               END DO 
     521            END DO 
     522            ! 
     523         END DO 
     524         ! 
     525      ENDIF ! End init==1 
     526 
     527      DO ji_sd = 1, ninco 
     528         ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 
    516529         DO ji2_sd = ji_sd+1, ninco 
    517             zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 
    518             DO jj_sd=1,ninco 
    519                ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 
    520             ENDDO 
    521          ENDDO 
    522  
    523       ENDDO 
    524  
    525    ENDIF ! End init==1 
    526  
    527    DO ji_sd = 1, ninco 
    528       ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 
    529       DO ji2_sd = ji_sd+1, ninco 
    530          ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 
    531       ENDDO 
    532    ENDDO 
    533  
    534    !system solving:  
    535    ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 
    536    ji_sd = ninco 
    537    DO ji_sd = ninco-1, 1, -1 
    538       zx1=0. 
    539       DO jj_sd = ji_sd+1, ninco 
    540          zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 
    541       ENDDO 
    542       ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 
    543    ENDDO 
    544  
    545    DO jj_sd =1, ninco 
    546       ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 
    547    ENDDO 
    548  
    549    CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
    550    CALL wrk_dealloc( jpincomax , ipos2 , ipivot        ) 
    551  
    552   END SUBROUTINE SUR_DETERMINE 
     530            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 
     531         END DO 
     532      END DO 
     533 
     534      !system solving:  
     535      ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 
     536      ji_sd = ninco 
     537      DO ji_sd = ninco-1, 1, -1 
     538         zx1 = 0._wp 
     539         DO jj_sd = ji_sd+1, ninco 
     540            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 
     541         END DO 
     542         ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 
     543      END DO 
     544 
     545      DO jj_sd =1, ninco 
     546         ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 
     547      END DO 
     548 
     549      CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
     550      CALL wrk_dealloc( jpincomax , ipos2 , ipivot        ) 
     551      ! 
     552   END SUBROUTINE SUR_DETERMINE 
    553553 
    554554#else 
Note: See TracChangeset for help on using the changeset viewer.