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 58 for trunk/NEMO/LIM_SRC – NEMO

Changeset 58 for trunk/NEMO/LIM_SRC


Ignore:
Timestamp:
2004-04-22T12:10:45+02:00 (20 years ago)
Author:
opalod
Message:

CT : BUGFIX032 : # Add a MAX function when computing temporary zusden array to avoid a floating-point zero division

# Remove the hard coded vertical lenght dz in limmsh.F90 and use fse3t(:,:,1) in limthd.F90
# Change variable name jeq to njeq and use a new algorithm to compute njeq based on Coriolis value

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/LIM_SRC/limmsh.F90

    r12 r58  
    11MODULE limmsh 
    2 #if defined key_ice_lim 
    32   !!====================================================================== 
    43   !!                     ***  MODULE  limmsh  *** 
    54   !! LIM ice model :   definition of the ice mesh parameters 
    65   !!====================================================================== 
    7  
     6#if defined key_ice_lim 
     7   !!---------------------------------------------------------------------- 
     8   !!   'key_ice_lim'                                     LIM sea-ice model 
    89   !!---------------------------------------------------------------------- 
    910   !!   lim_msh   : definition of the ice mesh 
     
    1415   USE dom_ice 
    1516   USE lbclnk 
     17   USE in_out_manager 
    1618 
    1719   IMPLICIT NONE 
     
    5153         !                   ! (resp. y direction) (defined at the center) 
    5254      REAL(wp) ::         & 
    53          zh1p  , zh2p  ,  &  ! Idem zh1, zh2 for the bottom left corner of the grid 
    54          zd2d1p, zd1d2p,  &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid 
    55          zusden, zusden2, &  ! temporary scalars 
    56          zaire4              !    "         " 
     55         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid 
     56         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid 
     57         zusden, zusden2     ! temporary scalars 
    5758      !!--------------------------------------------------------------------- 
    58       !!  LIM 2.0, UCL-LODYC-IPSL (2002) 
    59       !!--------------------------------------------------------------------- 
     59 
     60      IF(lwp) THEN 
     61         WRITE(numout,*) 
     62         WRITE(numout,*) 'lim_msh : LIM sea-ice model, mesh initialization' 
     63         WRITE(numout,*) '~~~~~~~' 
     64      ENDIF 
    6065       
    6166      !----------------------------------------------------------                           
     
    6368      !------------------------------------------------------------------  
    6469       
    65       jeq   = INT( jpj / 2 )   !i bug mpp potentiel 
    66       jeqm1 = jeq - 1  
     70      njeq   = INT( jpj / 2 )   !i bug mpp potentiel 
     71      njeqm1 = njeq - 1  
    6772 
    6873      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor 
    6974  
     75!i    DO jj = 1, jpj 
     76!i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line 
     77!!ii     write(numout,*) jj, zind(jj) 
     78!i    END DO 
     79 
     80      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere 
     81         l_jeq = .TRUE. 
     82         njeq  = 1 
     83         DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 ) 
     84            njeq = njeq + 1 
     85         END DO 
     86         IF(lwp ) WRITE(numout,*) '          the equator is inside the domain at about njeq = ', njeq 
     87      ELSEIF( fcor(1,1) < 0.e0 ) THEN 
     88         l_jeq = .FALSE. 
     89         njeq = jpj + 10 
     90         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the southern hemisphere: njeq = ', njeq 
     91      ELSE 
     92         l_jeq = .FALSE. 
     93         njeq = -10 
     94         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the northern hemisphere: njeq = ', njeq 
     95      ENDIF 
     96 
     97      njeqm1 = njeq - 1 
     98 
    7099 
    71100      !   For each grid, definition of geometric tables  
     
    73102       
    74103      !------------------- 
    75       ! Conventions :    | 
     104      ! Conventions :    ! 
    76105      !------------------- 
    77106      !  indices 1 \ 2 <-> localisation in the 2 direction x \ y 
     
    80109      !  3 = corner SW x(i-1/2),y(j-1/2) 
    81110      !------------------- 
     111!!ibug ??? 
     112      akappa(:,:,:,:) = 0.e0 
     113      wght(:,:,:,:) = 0.e0 
     114      alambd(:,:,:,:,:,:) = 0.e0 
     115      tmu(:,:) = 0.e0 
     116!!i 
    82117       
    83118       
     
    86121      !                                                       ! akappa 
    87122      DO jj = 2, jpj 
    88          DO ji = 1, jpi 
    89             zd1d2(ji,jj) = e1v(ji,jj) - e1v(ji,jj-1) 
    90          END DO 
     123         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1) 
    91124      END DO 
    92125      CALL lbc_lnk( zd1d2, 'T', -1. ) 
    93126 
    94       DO jj = 1, jpj 
    95          DO ji = 2, jpi 
    96             zd2d1(ji,jj) = e2u(ji,jj) - e2u(ji-1,jj) 
    97          END DO 
     127      DO ji = 2, jpi 
     128         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:) 
    98129      END DO 
    99130      CALL lbc_lnk( zd2d1, 'T', -1. ) 
    100131 
    101       DO jj = 1, jpj 
    102          DO ji = 1, jpi 
    103             zaire4            = 4.0 * e1t(ji,jj) * e2t(ji,jj) 
    104             akappa(ji,jj,1,1) =          1.0 / ( 2.0 * e1t(ji,jj) ) 
    105             akappa(ji,jj,1,2) = zd1d2(ji,jj) / zaire4 
    106             akappa(ji,jj,2,1) = zd2d1(ji,jj) / zaire4 
    107             akappa(ji,jj,2,2) =          1.0 / ( 2.0 * e2t(ji,jj) ) 
    108          END DO 
    109       END DO 
     132      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) ) 
     133      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) 
     134      akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) 
     135      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) ) 
    110136       
    111137      !                                                      ! weights (wght) 
     
    140166               &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 
    141167 
    142 ! better writen but change the last digit and thus solver in less than 100 timestep 
     168! better written but change the last digit and thus solver in less than 100 timestep 
    143169!           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   & 
    144170!              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)  
     
    147173!              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1) 
    148174 
    149             zusden = 1.0 / ( zh1p * zh2p * 4.e0 ) 
     175!!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 ) 
     176            zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) 
    150177            zusden2 = zusden * 2.0  
    151178 
     
    195222      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      ! 
    196223             
    197       ! Definition of scale dephts : bathymetry                
    198       !------------------------------------------- 
    199 !i bug  dz forced to 10 meters 
    200       dz = 10                  !!bug   potential bug if first level not equal to 10 m 
    201 !!!      dz = gdept(1) 
    202  
    203        
     224 
    204225      ! Initialization of ice masks 
    205226      !---------------------------- 
     
    221242      ! unmasked and masked area of T-grid cell 
    222243      area(:,:) = e1t(:,:) * e2t(:,:) 
    223       aire(:,:) = area(:,:) * tms(:,:) 
    224244       
    225245   END SUBROUTINE lim_msh 
     246 
    226247#else 
    227    !!============================================================================== 
    228    !!                       ***  MODULE limmsh   *** 
    229    !!                            No sea ice 
    230    !!============================================================================== 
     248   !!---------------------------------------------------------------------- 
     249   !!   Default option            Dummy Module         NO LIM sea-ice model 
     250   !!---------------------------------------------------------------------- 
    231251CONTAINS 
    232    SUBROUTINE lim_msh           ! Empty routine 
     252   SUBROUTINE lim_msh           ! Dummy routine 
    233253   END SUBROUTINE lim_msh 
    234  
    235254#endif 
     255 
    236256   !!====================================================================== 
    237257END MODULE limmsh 
Note: See TracChangeset for help on using the changeset viewer.