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 4333 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 – NEMO

Ignore:
Timestamp:
2013-12-11T18:34:00+01:00 (10 years ago)
Author:
clem
Message:

remove remaining bugs in LIM3, so that it can run in both regional and global config

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4298 r4333  
    4343   PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
    4444 
    45    REAL(wp) ::   epsi11 = 1.e-11_wp   ! constant values 
     45   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
    4646   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values 
    4747   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values 
     
    5050   ! Variables shared among ridging subroutines 
    5151   !----------------------------------------------------------------------- 
    52    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    53    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    54    ! 
    55    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
    56    !                                                           !  closing associated w/ category n 
    57    ! 
    58    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    59    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    60    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
    61    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice  
    62    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
    63    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
     52   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
     53   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
     54   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
     55   !                                                     ! closing associated w/ category n 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
     57   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
     58   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
     59   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice  
     60   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
     61   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
    6462 
    6563   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
    6664   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer 
    67    REAL(wp), PARAMETER ::   kamax   = 1.0 
     65   REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99) 
    6866 
    6967   REAL(wp) ::   Cp                             !  
     
    261259               ! Reduce the closing rate if more than 100% of the open water  
    262260               ! would be removed.  Reduce the opening rate proportionately. 
    263                IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
     261               IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
    264262                  w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    265263                  IF ( w1 .GT. ato_i(ji,jj)) THEN 
     
    281279            DO jj = 1, jpj 
    282280               DO ji = 1, jpi 
    283                   IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
     281                  IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
    284282                     w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    285283                     IF ( w1  >  a_i(ji,jj,jl) ) THEN 
     
    312310         DO jj = 1, jpj 
    313311            DO ji = 1, jpi 
    314                IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 
     312               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN 
    315313                  closing_net(ji,jj) = 0._wp 
    316314                  opning     (ji,jj) = 0._wp 
     
    354352         DO ji = 1, jpi 
    355353 
    356             IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 
     354            IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true. 
    357355 
    358356            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
     
    373371      DO jj = 1, jpj 
    374372         DO ji = 1, jpi 
    375             IF( ABS( asum(ji,jj) - kamax)  >  epsi11 ) THEN   ! there is a bug 
     373            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
    376374               WRITE(numout,*) ' ' 
    377375               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
     
    396394 
    397395      !-----------------------------------------------------------------------------! 
    398       ! 6) Updating state variables and trend terms 
     396      ! 6) Updating state variables and trend terms (done in limupdate) 
    399397      !-----------------------------------------------------------------------------! 
    400  
    401398      CALL lim_var_glo2eqv 
    402399      CALL lim_itd_me_zapsmall 
     
    415412         END DO 
    416413      END DO 
    417  
    418       !----------------- 
    419       !  Trend terms 
    420       !----------------- 
    421  
    422       d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
    423       d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
    424       d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
    425       d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
    426       d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
    427       d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
    428       d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:) 
    429       d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    430       d_smv_i_trp(:,:,:)   = 0._wp 
    431       IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    432414 
    433415      IF(ln_ctl) THEN     ! Control print 
     
    487469      ! ------------------------------- 
    488470 
    489       !-------------------------! 
    490       ! Back to initial values 
    491       !-------------------------! 
    492  
    493       ! update of fields will be made later in lim update 
    494       u_ice(:,:)   = old_u_ice(:,:) 
    495       v_ice(:,:)   = old_v_ice(:,:) 
    496       a_i(:,:,:)   = old_a_i(:,:,:) 
    497       v_s(:,:,:)   = old_v_s(:,:,:) 
    498       v_i(:,:,:)   = old_v_i(:,:,:) 
    499       e_s(:,:,:,:) = old_e_s(:,:,:,:) 
    500       e_i(:,:,:,:) = old_e_i(:,:,:,:) 
    501       oa_i(:,:,:)  = old_oa_i(:,:,:) 
    502       IF(  num_sal == 2  )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    503  
    504       !----------------------------------------------------! 
    505       ! Advection of ice in a free cell, newly ridged ice 
    506       !----------------------------------------------------! 
    507  
    508       ! to allow for thermodynamics to melt new ice 
    509       ! we immediately advect ice in free cells 
    510  
    511       ! heat content has to be corrected before ice volume 
    512 !clem@order 
    513 !      DO jl = 1, jpl 
    514 !         DO jk = 1, nlay_i 
    515 !            DO jj = 1, jpj 
    516 !               DO ji = 1, jpi 
    517 !                  IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
    518 !                     ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
    519 !                     old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
    520 !                     d_e_i_trp(ji,jj,jk,jl) = 0._wp 
    521 !                  ENDIF 
    522 !               END DO 
    523 !            END DO 
    524 !         END DO 
    525 !      END DO 
    526 ! 
    527 !      DO jl = 1, jpl 
    528 !         DO jj = 1, jpj 
    529 !            DO ji = 1, jpi 
    530 !               IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. & 
    531 !                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
    532 !                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
    533 !                  d_v_i_trp (ji,jj,jl)   = 0._wp 
    534 !                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
    535 !                  d_a_i_trp (ji,jj,jl)   = 0._wp 
    536 !                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
    537 !                  d_v_s_trp (ji,jj,jl)   = 0._wp 
    538 !                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
    539 !                  d_e_s_trp (ji,jj,1,jl) = 0._wp 
    540 !                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
    541 !                  d_oa_i_trp(ji,jj,jl)   = 0._wp 
    542 !                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
    543 !                  d_smv_i_trp(ji,jj,jl)  = 0._wp 
    544 !               ENDIF 
    545 !            END DO 
    546 !         END DO 
    547 !      END DO 
    548 !clem@order 
    549471      ENDIF  ! ln_limdyn=.true. 
    550472      ! 
     
    601523               DO ji = 1, jpi 
    602524                  ! 
    603                   IF(  a_i(ji,jj,jl)    > epsi11  .AND.     & 
    604                        athorn(ji,jj,jl) > 0._wp   ) THEN 
     525                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 
    605526                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    606527                     !---------------------------- 
     
    620541                        * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )    
    621542!!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                     
    622                   ENDIF            ! aicen > epsi11 
     543                  ENDIF            ! aicen > epsi10 
    623544                  ! 
    624545               END DO ! ji 
     
    677598         DO jj = 2, jpj - 1 
    678599            DO ji = 2, jpi - 1 
    679                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 
     600               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 
    680601                  ! present 
    681602                  zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
     
    717638         DO jj = 1, jpj - 1 
    718639            DO ji = 1, jpi - 1 
    719                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN       ! ice is present 
     640               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present 
    720641                  numts_rm = 1 ! number of time steps for the running mean 
    721642                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
     
    790711      DO jj = 1, jpj 
    791712         DO ji = 1, jpi 
    792             IF( ato_i(ji,jj) > epsi11 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj) 
     713            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj) 
    793714            ELSE                               ;   Gsum(ji,jj,0) = 0._wp 
    794715            ENDIF 
     
    800721         DO jj = 1, jpj  
    801722            DO ji = 1, jpi 
    802                IF( a_i(ji,jj,jl) .GT. epsi11 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
     723               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
    803724               ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 
    804725               ENDIF 
     
    883804      IF ( raftswi == 1 ) THEN 
    884805 
    885          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN 
     806         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN 
    886807            DO jl = 1, jpl 
    887808               DO jj = 1, jpj 
    888809                  DO ji = 1, jpi 
    889                      IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 
    890                         epsi11 ) THEN 
     810                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN 
    891811                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
    892812                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
     
    934854            DO ji = 1, jpi 
    935855 
    936                IF (a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
     856               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
    937857                  hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    938858                  hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin) 
     
    988908      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    989909      INTEGER ::   icells            ! number of cells with aicen > puny 
    990       REAL(wp) ::   zeps, zindb, zsrdg2   ! local scalar 
     910      REAL(wp) ::   zindb, zsrdg2   ! local scalar 
    991911      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration 
    992912 
     
    1040960 
    1041961      IF( con_i ) THEN 
    1042          CALL lim_column_sum (jpl,   v_i, vice_init ) 
    1043          WRITE(numout,*) ' vice_init  : ', vice_init(jiindx,jjindx) 
     962         CALL lim_column_sum        (jpl,    v_i,       vice_init ) 
    1044963         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init ) 
    1045          WRITE(numout,*) ' eice_init  : ', eice_init(jiindx,jjindx) 
     964         DO ji = mi0(jiindx), mi1(jiindx) 
     965            DO jj = mj0(jjindx), mj1(jjindx) 
     966               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj) 
     967               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj) 
     968            END DO 
     969         END DO 
    1046970      ENDIF 
    1047  
    1048       zeps   = 1.e-20_wp 
    1049971 
    1050972      !------------------------------------------------------------------------------- 
     
    1058980            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        & 
    1059981               &                        + opning(ji,jj)                          * rdt_ice 
    1060             IF( ato_i(ji,jj) < -epsi11 ) THEN 
     982            IF( ato_i(ji,jj) < -epsi10 ) THEN 
    1061983               neg_ato_i = .TRUE. 
    1062984            ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error 
     
    1070992         DO jj = 1, jpj  
    1071993            DO ji = 1, jpi 
    1072                IF( ato_i(ji,jj) < -epsi11 ) THEN  
     994               IF( ato_i(ji,jj) < -epsi10 ) THEN  
    1073995                  WRITE(numout,*) ''   
    1074996                  WRITE(numout,*) 'Ridging error: ato_i < 0' 
    1075997                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 
    1076                ENDIF               ! ato_i < -epsi11 
     998               ENDIF               ! ato_i < -epsi10 
    1077999            END DO 
    10781000         END DO 
     
    11161038         DO jj = 1, jpj 
    11171039            DO ji = 1, jpi 
    1118                IF( aicen_init(ji,jj,jl1)  >  epsi11  .AND.  athorn(ji,jj,jl1)  >  0._wp       & 
    1119                   .AND. closing_gross(ji,jj) > 0._wp ) THEN 
     1040               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  & 
     1041                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN 
    11201042                  icells = icells + 1 
    11211043                  indxi(icells) = ji 
     
    11541076            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    11551077 
    1156             IF (afrac(ji,jj) > kamax + epsi11) THEN  !riging 
     1078            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging 
    11571079               large_afrac = .true. 
    11581080            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error 
    11591081               afrac(ji,jj) = kamax 
    11601082            ENDIF 
    1161             IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 
     1083            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 
    11621084               large_afrft = .true. 
    11631085            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error 
     
    12181140            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    12191141            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1220             !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
    12211142            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    12221143 
     
    12721193 
    12731194               ! corrected sea water salinity 
    1274                zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - zeps ) ) 
    1275                zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), zeps ) 
     1195               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 
     1196               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 
    12761197 
    12771198               ztmelts          = - tmut * zdummy + rtt 
     
    13081229               ji = indxi(ij) 
    13091230               jj = indxj(ij) 
    1310                IF( afrac(ji,jj) > kamax + epsi11 ) THEN  
     1231               IF( afrac(ji,jj) > kamax + epsi10 ) THEN  
    13111232                  WRITE(numout,*) '' 
    13121233                  WRITE(numout,*) ' ardg > a_i' 
     
    13201241               ji = indxi(ij) 
    13211242               jj = indxj(ij) 
    1322                IF( afrft(ji,jj) > kamax + epsi11 ) THEN  
     1243               IF( afrft(ji,jj) > kamax + epsi10 ) THEN  
    13231244                  WRITE(numout,*) '' 
    13241245                  WRITE(numout,*) ' arft > a_i' 
     
    14201341         fieldid = ' v_i : limitd_me ' 
    14211342         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)  
    1422          WRITE(numout,*) ' vice_init  : ', vice_init(jiindx,jjindx) 
    1423          WRITE(numout,*) ' vice_final : ', vice_final(jiindx,jjindx) 
    14241343 
    14251344         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final ) 
    14261345         fieldid = ' e_i : limitd_me ' 
    14271346         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)  
    1428          WRITE(numout,*) ' eice_init  : ', eice_init(jiindx,jjindx) 
    1429          WRITE(numout,*) ' eice_final : ', eice_final(jiindx,jjindx) 
     1347 
     1348         DO ji = mi0(jiindx), mi1(jiindx) 
     1349            DO jj = mj0(jjindx), mj1(jjindx) 
     1350               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj) 
     1351               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) 
     1352               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj) 
     1353               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj) 
     1354            END DO 
     1355         END DO 
    14301356      ENDIF 
    14311357      ! 
     
    15361462 
    15371463      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace 
    1538        
     1464      REAL(wp)                          ::   zmask_glo 
    15391465!!gm      REAL(wp) ::   xtmp      ! temporary variable 
    15401466      !!------------------------------------------------------------------- 
     
    15481474         ! Abort model in case of negative area. 
    15491475         !----------------------------------------------------------------- 
    1550          IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 .AND. ln_nicep ) THEN 
    1551             DO jj = 1, jpj 
    1552                DO ji = 1, jpi 
    1553                   IF ( a_i(ji,jj,jl) .LT. -epsi11 ) THEN 
    1554                      WRITE (numout,*) ' ALERTE 98 '  
    1555                      WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl 
    1556                      WRITE (numout,*) ' a_i    ', a_i(ji,jj,jl) 
    1557                   ENDIF 
    1558                END DO 
    1559             END DO 
    1560          ENDIF 
    1561  
    15621476         icells = 0 
    1563          zmask  = 0._wp 
     1477         zmask(:,:)  = 0._wp 
    15641478         DO jj = 1, jpj 
    15651479            DO ji = 1, jpi 
    1566                IF(  ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0._wp   ) .OR.   & 
    1567                   & ( a_i(ji,jj,jl) .GT. 0._wp   .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) .OR.   & 
    1568                   & ( v_i(ji,jj,jl)  ==  0._wp   .AND. a_i(ji,jj,jl) .GT. 0._wp   ) .OR.   & 
    1569                   & ( v_i(ji,jj,jl) .GT. 0._wp   .AND. v_i(ji,jj,jl) .LT. 1.e-12  )      )   zmask(ji,jj) = 1._wp 
    1570             END DO 
    1571          END DO 
    1572          IF( ln_nicep )   WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 
     1480               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   & 
     1481                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
     1482                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   & 
     1483                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
     1484            END DO 
     1485         END DO 
     1486         zmask_glo = glob_sum(zmask) 
     1487         !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean ' 
    15731488 
    15741489         !----------------------------------------------------------------- 
     
    15821497!!gm                  xtmp = xtmp * unit_fac 
    15831498                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    1584                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) 
     1499                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
    15851500               END DO 
    15861501            END DO 
     
    16041519!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ??????? 
    16051520 
    1606                t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 
     1521               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    16071522 
    16081523               !----------------------------------------------------------------- 
     
    16181533 
    16191534               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj) 
    1620                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1621                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1622                v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1623                t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1 - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
    1624                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1625                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1535               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1536               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1537               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1538               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
     1539               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1540               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    16261541               ! 
    16271542            END DO 
Note: See TracChangeset for help on using the changeset viewer.