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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 865 for trunk/NEMO/LIM_SRC_3/limthd_lac.F90 – NEMO

Ignore:
Timestamp:
2008-03-14T15:43:09+01:00 (16 years ago)
Author:
rblod
Message:

Correct a bug and clean comments in limthd_lac, see ticket #79

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/LIM_SRC_3/limthd_lac.F90

    r834 r865  
    227227! 3) Collection thickness of ice formed in leads and polynyas 
    228228!------------------------------------------------------------------------------!     
    229       ! hicol is the thickness of new ice. 
     229      ! hicol is the thickness of new ice formed in open water 
     230      ! hicol can be either prescribed (frazswi = 0) 
     231      ! or computed (frazswi = 1) 
    230232      ! Frazil ice forms in open water, is transported by wind 
    231233      ! accumulates at the edge of the consolidated ice edge 
    232234      ! where it forms aggregates of a specific thickness called 
    233235      ! collection thickness. 
     236 
     237      ! Note : the following algorithm currently breaks vectorization 
     238      !  
    234239 
    235240      zvrel(:,:) = 0.0 
     
    253258          zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag) 
    254259          zgamafr = 0.03 
    255  
    256           !+++++ 
    257           WRITE(numout,*) ' ztwogp : ', ztwogp 
    258           WRITE(numout,*) ' rau0   : ', rau0 
    259           WRITE(numout,*) ' grav   : ', grav 
    260           WRITE(numout,*) ' rhoic  : ', rhoic 
    261           !+++++ 
    262260 
    263261          DO jj = 1, jpj 
     
    275273            ! Square root of wind stress 
    276274            ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
    277 !           !+++++ 
    278             ! tension on B-grid 
    279 !           ztaux         = ( gtaux(ji  ,jj  ) * tmu(ji,jj    ) & 
    280 !                         +   gtaux(ji+1,jj  ) * tmu(ji+1,jj  ) & 
    281 !                         +   gtaux(ji,jj+1  ) * tmu(ji,jj+1  ) & 
    282 !                         +   gtaux(ji+1,jj+1) * tmu(ji+1,jj+1) ) / 4.0 
    283 !           ztauy         = ( gtauy(ji  ,jj  ) * tmu(ji,jj    ) & 
    284 !                         +   gtauy(ji+1,jj  ) * tmu(ji+1,jj  ) & 
    285 !                         +   gtauy(ji,jj+1  ) * tmu(ji,jj+1  ) & 
    286 !                         +   gtauy(ji+1,jj+1) * tmu(ji+1,jj+1) ) / 4.0 
    287 !           !+++++ 
    288 !           !+++++ 
    289 !           IF ( ( ji.EQ.jiindex ) .AND. ( jj.EQ.jjindex ) ) THEN 
    290 !              WRITE(numout,*) 
    291 !              WRITE(numout,*) ' ztaux : ', ztaux 
    292 !              WRITE(numout,*) ' ztauy : ', ztauy 
    293 !              WRITE(numout,*) ' |tau| : ', ztenagm 
    294 !           ENDIF 
    295 !           !+++++ 
    296275 
    297276            !--------------------- 
     
    320299            zvrel(ji,jj)  = SQRT(zvrel2) 
    321300 
    322             !+++++++++ 
    323 !           ! ice drift on B-grid 
    324 !           zindb = MAX(0.0, SIGN(1.0, at_i(ji,jj) )) 
    325 !           zvgx  = zindb * (  u_ice(ji,jj    )   * tmu(ji,jj  )     & 
    326 !                           +  u_ice(ji+1,jj  )   * tmu(ji+1,jj)     & 
    327 !                           +  u_ice(ji,jj+1  )   * tmu(ji,jj+1)     & 
    328 !                           +  u_ice(ji+1,jj+1)   * tmu(ji+1,jj+1) ) & 
    329 !                           / 4.0 
    330 !           zvgy  = zindb * (  v_ice(ji,jj    )   * tmu(ji,jj  ) & 
    331 !                           +  v_ice(ji+1,jj  )   * tmu(ji+1,jj) & 
    332 !                           +  v_ice(ji,jj+1  )   * tmu(ji,jj+1) & 
    333 !                           +  v_ice(ji+1,jj+1)   * tmu(ji+1,jj+1) ) & 
    334 !                           / 4.0 
    335  
    336 !           !+++++ 
    337 !           IF ( ( ji.EQ.jiindex ) .AND. ( jj.EQ.jjindex ) ) THEN 
    338 !              WRITE(numout,*) 
    339 !              WRITE(numout,*) ' zvfrx : ', zvfrx 
    340 !              WRITE(numout,*) ' zvfry : ', zvfry 
    341 !              WRITE(numout,*) ' zvgx  : ', zvgx  
    342 !              WRITE(numout,*) ' zvgy  : ', zvgx  
    343 !           ENDIF 
    344 !           !+++++ 
    345  
    346 !           !+++++ 
    347 !           IF ( ( ji.EQ.jiindex ) .AND. ( jj.EQ.jjindex ) ) THEN 
    348 !              WRITE(numout,*) 
    349 !              WRITE(numout,*) ' zvrel2: ', zvrel2 
    350 !              WRITE(numout,*) ' zvrel : ', zvrel(ji,jj) 
    351 !           ENDIF 
    352  
    353 !           hicol(ji,jj) = 10.0 ! starting value has to be high!!! 
    354  
    355301            !--------------------- 
    356302            ! Iterative procedure 
     
    372318               hicol(ji,jj)   = zhicol_new 
    373319 
    374                !++++++++++++++++++++ 
    375 !              IF ( ABS(-zf/zfp) .LT. 1.0d-9 ) THEN 
    376 !                 iterate_frazil = .false.  
    377 !              ENDIF 
    378  
    379 !              IF ( (ji.eq.jiindex).AND.(jj.eq.jjindex) ) THEN 
    380 !                      WRITE(numout,*) ' ------- iter = ', iter, & 
    381 !                      '----------------' 
    382 !                      WRITE(numout,*) ' iterate_frazil ', iterate_frazil 
    383 !                      WRITE(numout,*) ' zf, zfp : ', zf, zfp 
    384 !                      WRITE(numout,*) ' dhicol  : ', -zf/zfp 
    385 !                      WRITE(numout,*) ' hicol   : ', hicol(ji,jj) 
    386                !++++++++++++++++++++ 
    387 !              ENDIF 
    388  
    389320               iter = iter + 1 
    390321 
    391322            END DO ! do while 
    392 !         hicol(ji,jj) = MAX( hicol(ji,jj) , 0.03) 
    393 !         WRITE(numout,*) ' zvrel2 : ', zvrel2, ' zvrel : ', SQRT(zvrel2), ' hicol : ',hicol(ji,jj) 
    394323 
    395324          ENDIF ! end of selection of pixels where ice forms 
    396  
    397           !+++++++++++ 
    398           IF ( hicol(ji,jj) .GT. 2.00 .OR. hicol(ji,jj) .LT. 0.0 ) THEN 
    399                   WRITE(numout,*) ' ALERTE 125 : hicol too bad ' 
    400                   WRITE(numout,*) ' ji,jj : ', ji, jj 
    401                   WRITE(numout,*) ' lat   : ', gphit(ji,jj) 
    402                   WRITE(numout,*) ' lon   : ', glamt(ji,jj) 
    403                   WRITE(numout,*) ' hicol : ', hicol(ji,jj) 
    404           ENDIF 
    405           !+++++++++++ 
    406325 
    407326      END DO ! loop on ji ends 
    408327      END DO ! loop on jj ends 
    409  
    410       WRITE(numout,*) ' ' 
    411       WRITE(numout,*) ' Apres calcul : ' 
    412       WRITE(numout,*) ' hicol      :   ', hicol(jiindex,jjindex) 
    413       WRITE(numout,*) ' at_i       :   ', at_i(jiindex,jjindex) 
    414328 
    415329      ENDIF ! End of computation of frazil ice collection thickness 
     
    489403        END DO 
    490404        IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:) 
    491  
    492         !+++++++++++++ 
    493         DO ji = 1, nbpac 
    494            IF (zh_newice(ji) .LE. 0.0) THEN 
    495               zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    496               zjj                  = ( npac(ji) - 1 ) / jpi + 1 
    497               WRITE(numout,*) ' collection thickness  <= 0 ', zh_newice(ji), ji, zji, zjj 
    498            ENDIF 
    499         END DO 
    500         !+++++++++++++ 
    501405 
    502406        !---------------------- 
     
    563467           zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
    564468        END DO 
     469 
     470        !--------------------------------- 
     471        ! Salt flux due to new ice growth 
     472        !--------------------------------- 
     473        IF ( ( num_sal .EQ. 4 ) ) THEN  
     474           DO ji = 1, nbpac 
     475              zji            = MOD( npac(ji) - 1, jpi ) + 1 
     476              zjj            = ( npac(ji) - 1 ) / jpi + 1 
     477              fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
     478                               ( sss_io(zji,zjj) - bulk_sal      ) * rhoic *      & 
     479                               zv_newice(ji) / rdt_ice 
     480           END DO 
     481        ELSE 
     482           DO ji = 1, nbpac 
     483              zji            = MOD( npac(ji) - 1, jpi ) + 1 
     484              zjj            = ( npac(ji) - 1 ) / jpi + 1 
     485              fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
     486                               ( sss_io(zji,zjj) - zs_newice(ji) ) * rhoic *      & 
     487                               zv_newice(ji) / rdt_ice 
     488           END DO ! ji 
     489        ENDIF 
    565490 
    566491        !------------------------------------ 
     
    804729        ENDIF ! num_sal 
    805730   
    806         !--------------------------------- 
    807         ! Salt flux due to new ice growth 
    808         !--------------------------------- 
    809         IF ( ( num_sal .EQ. 4 ) ) THEN  
    810            DO ji = 1, nbpac 
    811               zji            = MOD( npac(ji) - 1, jpi ) + 1 
    812               zjj            = ( npac(ji) - 1 ) / jpi + 1 
    813               fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
    814                                ( sss_io(zji,zjj) - bulk_sal      ) * rhoic *      & 
    815                                zv_newice(ji) / rdt_ice 
    816            END DO 
    817         ELSE 
    818            DO ji = 1, nbpac 
    819               zji            = MOD( npac(ji) - 1, jpi ) + 1 
    820               zjj            = ( npac(ji) - 1 ) / jpi + 1 
    821               fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
    822                                ( sss_io(zji,zjj) - zs_newice(ji) ) * rhoic *      & 
    823                                zv_newice(ji) / rdt_ice 
    824            END DO ! ji 
    825         ENDIF 
    826731 
    827732!------------------------------------------------------------------------------! 
Note: See TracChangeset for help on using the changeset viewer.