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 – NEMO

Changeset 4333


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

Location:
branches/2013/dev_MERGE_2013/NEMOGCM
Files:
1 deleted
30 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/cfg.txt

    r4310 r4333  
    11GYRE_BFM OPA_SRC TOP_SRC 
    22GYRE_PISCES OPA_SRC TOP_SRC 
    3 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
    43ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    54GYRE OPA_SRC 
    … …  
    109ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    1110ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
     11ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4205 r4333  
    292292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s 
    293293 
    294    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   patho_case ! number of the pathological case (if any, of course) 
    295  
    296294   !!-------------------------------------------------------------------------- 
    297295   !! * Ice global state variables 
    … …  
    465463         &      fsup2D     (jpi,jpj) , focea2D    (jpi,jpj) , q_s   (jpi,jpj) , STAT=ierr(ii) ) 
    466464 
    467       ii = ii + 1 
    468       ALLOCATE( patho_case(jpi, jpj, jpl)  , STAT=ierr(ii) ) 
    469  
    470465      ! * Ice global state variables 
    471466      ii = ii + 1 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r4319 r4333  
    9090      ! 
    9191      !                                ! Initial sea-ice state 
    92       IF( .NOT.ln_rstart ) THEN              ! start from rest 
     92      IF( .NOT. ln_rstart ) THEN              ! start from rest 
    9393         numit = 0 
    9494         numit = nit000 - 1 
    … …  
    105105      hi_max(jpl) = 99._wp  ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 
    106106      ! 
    107       CALL lim_sbc_init                ! ice surface boundary condition    
    108       ! 
    109       fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction 
    110       tn_ice(:,:,:) = t_su(:,:,:) 
     107      CALL lim_sbc_init                 ! ice surface boundary condition    
     108      ! 
     109      fr_i(:,:) = at_i(:,:)             ! initialisation of sea-ice fraction 
     110      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
    111111      ! 
    112112      nstart = numit  + nn_fsbc       
    … …  
    143143      WRITE ( numoni, namicerun ) 
    144144      ! 
    145       IF( lk_mpp .AND. ln_nicep ) THEN 
    146          ln_nicep = .FALSE. 
    147          CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
    148       ENDIF 
     145      !IF( lk_mpp .AND. ln_nicep ) THEN 
     146      !   ln_nicep = .FALSE. 
     147      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
     148      !ENDIF 
    149149      ! 
    150150      IF(lwp) THEN                        ! control print 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90

    r4298 r4333  
    4545      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
    4646      REAL(wp) :: epsi06 = 1.0e-6 
    47       REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars 
     47      REAL(wp) :: zc1, zc2, zc3, zx1, zdh   ! local scalars 
    4848      REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space 
    4949  
    … …  
    8585      ! ---------------------------------------- 
    8686      DO ji = 1, ijpij 
    87         ! snow thickness in each category  
    88          zht_s(ji,1:jpl) = zhts(ji) 
    8987          
     88         IF( zhti(ji) > 0._wp ) THEN 
     89 
    9090         ! initialisation of tests 
    9191         ztest_1 = 0 
    … …  
    108108               zht_i(ji,1) = zhti(ji) 
    109109               za_i (ji,1) = zai (ji) 
    110                 
    111110               ! *** case ice is thicker: fill categories >1 
    112111            ELSE 
    … …  
    178177         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 
    179178         !ENDIF 
    180              
     179          
     180         ENDIF ! if zhti > 0 
     181 
    181182      END DO ! i loop 
     183 
     184      ! ------------------------------------------------ 
     185      ! Adding Snow in each category where za_i is not 0 
     186      ! ------------------------------------------------  
     187      DO jl = 1, jpl 
     188         DO ji = 1, ijpij 
     189            IF( za_i(ji,jl) > 0._wp ) THEN 
     190               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     191               ! In case snow load is in excess that would lead to transformation from snow to ice 
     192               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     193               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )  
     194               ! recompute ht_i, ht_s avoiding out of bounds values 
     195               zht_i(ji,jl) = MIN( zhi_max(jl), zht_i(ji,jl) + zdh ) 
     196               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 
     197            ENDIF 
     198         ENDDO 
     199      ENDDO 
    182200 
    183201      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D') 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r4161 r4333  
    2929 
    3030   LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call) 
    31    REAL(wp) ::   epsi04 = 1e-04              ! constant 
     31   REAL(wp) ::   epsi04 = 1.e-04              ! constant 
    3232   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient 
    3333 
    … …  
    5454      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied 
    5555      ! 
    56       INTEGER  ::  ji, jj            ! dummy loop indices 
    57       INTEGER  ::  its, iter, ierr   ! local integers 
    58       REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars 
     56      INTEGER  ::  ji, jj                   ! dummy loop indices 
     57      INTEGER  ::  its, iter, ierr          ! local integers 
     58      REAL(wp) ::   zalfa, zrlxint, zconv   ! local scalars 
    5959      REAL(wp), POINTER, DIMENSION(:,:) ::   zrlx, zflu, zflv, zdiv0, zdiv, ztab0 
    6060      CHARACTER(lc) ::   charout   ! local character 
    … …  
    7979      zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
    8080      its   = 100                         ! Maximum number of iteration 
    81       zeps  =  2._wp * epsi04 
    8281      ! 
    8382      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization 
    … …  
    9493      iter  = 0 
    9594      ! 
    96       DO WHILE( zconv > zeps .AND. iter <= its )   ! Sub-time step loop 
     95      DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its )   ! Sub-time step loop 
    9796         ! 
    9897         iter = iter + 1                                 ! incrementation of the sub-time step number 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4205 r4333  
    7676      !!                6) Lateral boundary conditions 
    7777      !! 
     78      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even 
     79      !!              where there is no ice (clem: I do not know why but it is mandatory)  
     80      !! 
    7881      !! History : 
    7982      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
    … …  
    8487      !! * Local variables 
    8588      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
    86       REAL(wp)   :: epsi06, epsi20, ztmelts 
     89      REAL(wp)   :: epsi20, ztmelts, zdh 
    8790      INTEGER    :: i_hemis, i_fill, jl0   
    8891      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
    … …  
    98101      CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
    99102 
    100       epsi06   = 1.0e-6 
    101103      epsi20   = 1.0e-20 
    102104      IF(lwp) WRITE(numout,*) 
    … …  
    308310      !--------------------------------------------------------------------- 
    309311 
    310       ! Ice concentration, thickness and volume, snow depth, ice 
    311       ! salinity, ice age, surface temperature 
     312      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 
    312313      DO jl = 1, jpl ! loop over categories 
    313314         DO jj = 1, jpj 
    … …  
    315316               a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
    316317               ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness 
    317                ht_s(ji,jj,jl)  = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj))     ! snow depth 
    318                sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity 
     318               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) )  ! snow depth 
     319               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity 
    319320               o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 
    320                t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp 
    321   
    322                ! ice volume, snow volume, salt content, age content 
     321               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! surf temp 
     322 
     323               ! This case below should not be used if (ht_s/ht_i) is ok in namelist 
     324               ! In case snow load is in excess that would lead to transformation from snow to ice 
     325               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     326               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )  
     327               ! recompute ht_i, ht_s avoiding out of bounds values 
     328               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
     329               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
     330 
     331               ! ice volume, salt content, age content 
    323332               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
    324333               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 

  • 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 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r4293 r4333  
    4545   PUBLIC   lim_itd_shiftice 
    4646 
    47    REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
    48    REAL(wp) ::   epsi13 = 1.e-13_wp   ! 
    4947   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
     48   REAL(wp) ::   epsi06 = 1.e-6_wp   ! 
    5049 
    5150   !!---------------------------------------------------------------------- 
    … …  
    110109      CALL lim_var_glo2eqv    ! only for info 
    111110 
    112       !---------------------------------------------------------------------------------------- 
    113       !  4) Computation of trend terms and get back to old values       
    114       !---------------------------------------------------------------------------------------- 
    115  
    116       !- Trend terms 
    117       d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:)  
    118       d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:) 
    119       d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)   
    120       d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
    121       d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    122       !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:) 
    123       d_smv_i_thd(:,:,:) = 0._wp 
    124       IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
    125  
    126       ! diag only (clem) 
    127       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    128  
    129       IF(ln_ctl) THEN   ! Control print 
     111     IF(ln_ctl) THEN   ! Control print 
    130112         CALL prt_ctl_info(' ') 
    131113         CALL prt_ctl_info(' - Cell values : ') 
    … …  
    183165      ! ------------------------------- 
    184166      ! 
    185       !- Recover Old values 
    186       a_i(:,:,:)   = old_a_i (:,:,:) 
    187       v_s(:,:,:)   = old_v_s (:,:,:) 
    188       v_i(:,:,:)   = old_v_i (:,:,:) 
    189       e_s(:,:,:,:) = old_e_s (:,:,:,:) 
    190       e_i(:,:,:,:) = old_e_i (:,:,:,:) 
    191       !?? oa_i(:,:,:)  = old_oa_i(:,:,:) 
    192       IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    193       ! 
    194       IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
     167     IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
    195168   END SUBROUTINE lim_itd_th 
    196169   ! 
    … …  
    281254         DO jj = 1, jpj 
    282255            DO ji = 1, jpi 
    283                zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10))     !0 if no ice and 1 if yes 
    284                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 
    285                zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
    286                zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 
    287                IF( a_i(ji,jj,jl) > 1e-6 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
     256               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
     257               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 
     258               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
     259               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb 
     260               IF( a_i(ji,jj,jl) > epsi06 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
    288261            END DO 
    289262         END DO 
    … …  
    329302            ij = nind_j(ji) 
    330303            ! 
    331             IF ( ( zht_i_o(ii,ij,jl)  .GT.epsi10 ) .AND. &  
    332                ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN 
     304            IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &  
     305               ( zht_i_o(ii,ij,jl+1) .GT. epsi10 ) ) THEN 
    333306               !interpolate between adjacent category growth rates 
    334307               zslope = ( zdhice(ii,ij,jl+1)     - zdhice(ii,ij,jl) ) / & 
    … …  
    449422                     * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    450423                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 
    451                   v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ? 
     424                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 
    452425               ENDIF     ! zetamax > 0 
    453426               ! ji, a_i > epsi10 
    … …  
    541514            a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim  
    542515            ht_i(ii,ij,1) = hiclim 
    543             v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless 
     516            v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless 
    544517         ENDIF 
    545518      END DO !ji 
    … …  
    604577      REAL(wp) ::   zdhr         ! 1 / (hR - hL) 
    605578      REAL(wp) ::   zwk1, zwk2   ! temporary variables 
    606       REAL(wp) ::   zacrith      ! critical minimum concentration in an ice category 
    607       !!------------------------------------------------------------------ 
    608       ! 
    609       zacrith       = 1.0e-6 
     579      !!------------------------------------------------------------------ 
     580      ! 
    610581      ! 
    611582      DO jj = 1, jpj 
    612583         DO ji = 1, jpi 
    613584            ! 
    614             IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > zacrith   & 
     585            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   & 
    615586               &                        .AND. hice(ji,jj)        > 0._wp     ) THEN 
    616587 
    … …  
    1015986                  ! end TECLIM change  
    1016987                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary 
    1017                   zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) ) / ht_i(ji,jj,jl)   
     988                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl)   
    1018989                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 
    1019990               ENDIF 
    … …  
    10431014         zshiftflag = 0 
    10441015 
     1016!clem-change 
     1017!         DO jj = 1, jpj 
     1018!            DO ji = 1, jpi 
     1019!               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   & 
     1020!                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
     1021!                  ! 
     1022!                  zshiftflag = 1 
     1023!                  zdonor(ji,jj,jl) = jl + 1 
     1024!                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
     1025!                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     1026!               ENDIF 
     1027!            END DO                 ! ji 
     1028!         END DO                 ! jj 
     1029! 
     1030!         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     1031!          
     1032!         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     1033!            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     1034!            ! Reset shift parameters 
     1035!            zdonor(:,:,jl) = 0 
     1036!            zdaice(:,:,jl) = 0._wp 
     1037!            zdvice(:,:,jl) = 0._wp 
     1038!         ENDIF 
     1039!clem-change 
     1040 
     1041         ! clem-change begin: why not doing that? 
    10451042         DO jj = 1, jpj 
    10461043            DO ji = 1, jpi 
    10471044               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   & 
    10481045                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    1049                   ! 
    1050                   zshiftflag = 1 
    1051                   zdonor(ji,jj,jl) = jl + 1 
    1052                   zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
    1053                   zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     1046                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
     1047                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    10541048               ENDIF 
    10551049            END DO                 ! ji 
    10561050         END DO                 ! jj 
    1057  
    1058          IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
    1059           
    1060          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    1061             CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    1062             ! Reset shift parameters 
    1063             zdonor(:,:,jl) = 0 
    1064             zdaice(:,:,jl) = 0._wp 
    1065             zdvice(:,:,jl) = 0._wp 
    1066          ENDIF 
     1051         ! clem-change end 
    10671052 
    10681053      END DO                    ! jl 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r4269 r4333  
    5050   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2) 
    5151 
     52   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    5253   REAL(wp) ::   rzero   = 0._wp   ! constant values 
    5354   REAL(wp) ::   rone    = 1._wp   ! constant values 
    … …  
    426427               !-Calculate stress tensor components zs1 and zs2  
    427428               !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    428                !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   & 
    429                !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &        
    430                !   &          / ( 1._wp + alphaevp * dtotel ) 
    431  
    432                !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   & 
    433                !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   & 
    434                !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
     429               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   & 
     430                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &        
     431                  &          / ( 1._wp + alphaevp * dtotel ) 
     432 
     433               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   & 
     434                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   & 
     435                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
    435436 
    436437               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
    437                zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  & 
    438                   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
    439                zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  & 
    440                   &         / ( 1._wp + dtotel ) 
     438               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  & 
     439               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
     440               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  & 
     441               !   &         / ( 1._wp + dtotel ) 
    441442 
    442443            END DO 
    … …  
    472473 
    473474               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 
    474                !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  & 
    475                !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
    476                !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )  
     475               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  & 
     476                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     477                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )  
    477478 
    478479               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
    479                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
    480                   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
    481                   &          / ( 1.0 + dtotel )  
     480               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
     481               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     482               !   &          / ( 1.0 + dtotel )  
    482483 
    483484            END DO ! ji 
    … …  
    485486 
    486487         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    487  
    488 !#if defined key_bdy 
    489 !         ! clem: change zs1, zs2, zs12 at the boundary for each iteration 
    490 !         CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 ) 
    491 !         CALL lbc_lnk( zs1 (:,:), 'T', 1. ) 
    492 !         CALL lbc_lnk( zs2 (:,:), 'T', 1. ) 
    493 !         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    494 !#endif          
    495488 
    496489         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
    … …  
    544537            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    545538#endif 
     539#if defined key_bdy 
     540         ! clem: change u_ice and v_ice at the boundary for each iteration 
     541         CALL bdy_ice_lim_dyn( 'U' ) 
     542#endif          
    546543 
    547544!CDIR NOVERRCHK 
    … …  
    572569            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
    573570#endif 
     571#if defined key_bdy 
     572         ! clem: change u_ice and v_ice at the boundary for each iteration 
     573         CALL bdy_ice_lim_dyn( 'V' ) 
     574#endif          
    574575 
    575576         ELSE  
    … …  
    599600            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    600601#if defined key_agrif && defined key_lim2 
    601             CALL agrif_rhg_lim2( jter, nevp , 'V' ) 
    602 #endif 
     602            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
     603#endif 
     604#if defined key_bdy 
     605         ! clem: change u_ice and v_ice at the boundary for each iteration 
     606         CALL bdy_ice_lim_dyn( 'V' ) 
     607#endif          
    603608 
    604609!CDIR NOVERRCHK 
    … …  
    632637            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    633638#endif 
     639#if defined key_bdy 
     640         ! clem: change u_ice and v_ice at the boundary for each iteration 
     641         CALL bdy_ice_lim_dyn( 'U' ) 
     642#endif          
    634643 
    635644         ENDIF 
    636645          
    637 !#if defined key_bdy 
    638 !         ! clem: change u_ice and v_ice at the boundary for each iteration 
    639 !         CALL bdy_ice_lim_dyn( 1 ) 
    640 !#endif          
    641  
    642646         IF(ln_ctl) THEN 
    643647            !---  Convergence test. 
    … …  
    666670!CDIR NOVERRCHK 
    667671         DO ji = fs_2, fs_jpim1 
    668             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    669             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     672            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     673            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    670674            zdummy = vt_i(ji,jj) 
    671675            IF ( zdummy .LE. hminrhg ) THEN 
    … …  
    684688#if defined key_bdy 
    685689      ! clem: change u_ice and v_ice at the boundary 
    686       CALL bdy_ice_lim_dyn( 1 ) 
     690      CALL bdy_ice_lim_dyn( 'U' ) 
     691      CALL bdy_ice_lim_dyn( 'V' ) 
    687692#endif          
    688693 
    689694      DO jj = k_j1+1, k_jpj-1  
    690695         DO ji = fs_2, fs_jpim1 
    691             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    692             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     696            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     697            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    693698            zdummy = vt_i(ji,jj) 
    694699            IF ( zdummy .LE. hminrhg ) THEN 
    … …  
    714719            !- zdd(:,:), zdt(:,:): divergence and tension at centre  
    715720            !- zds(:,:): shear on northeast corner of grid cells 
    716             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    717             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     721            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     722            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    718723            zdummy = vt_i(ji,jj) 
    719724            IF ( zdummy .LE. hminrhg ) THEN 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4302 r4333  
    5151   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim 
    5252 
    53    REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values 
    5453   REAL(wp)  ::   rzero  = 0._wp     
    5554   REAL(wp)  ::   rone   = 1._wp 
    … …  
    147146 
    148147            !   computation the solar flux at ocean surface 
    149             IF (lk_cpl) THEN ! be carfeful: not being tested yet 
     148            IF (lk_cpl) THEN ! be carfeful: not been tested yet 
    150149               ! original line 
    151150               !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) 
    … …  
    190189            qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux 
    191190            !                           ! fdtcn : turbulent oceanic heat flux 
    192  
    193             !!gm   this IF prevents the vertorisation of the whole loop 
    194           !  IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
    195           !     WRITE(numout,*) ' lim_sbc : heat fluxes ' 
    196           !     WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx) 
    197           !     WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx) 
    198           !     WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx) 
    199           !     WRITE(numout,*) 
    200           !     WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx) 
    201           !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    202           !     WRITE(numout,*) ' ifral     : ', ifral 
    203           !     WRITE(numout,*) ' ial       : ', ial   
    204           !     WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx) 
    205           !     WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx) 
    206           !     !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 
    207           !     !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 
    208           !     WRITE(numout,*) ' ifrdv     : ', ifrdv 
    209           !     WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx) 
    210           !     WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx) 
    211           !     !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 
    212           !     !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 
    213           !     WRITE(numout,*) ' ' 
    214           !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    215           !     WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx) 
    216           !     WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx) 
    217           !     WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx) 
    218           !     WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
    219           !  ENDIF 
    220             !!gm   end 
    221191         END DO 
    222192      END DO 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r4295 r4333  
    5050   PUBLIC   lim_thd_init   ! called by iceini module 
    5151 
    52    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    53    REAL(wp) ::   epsi16 = 1e-16_wp   ! 
    54    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    55    REAL(wp) ::   epsi06 = 1e-06_wp   ! 
    56    REAL(wp) ::   epsi04 = 1e-04_wp   ! 
     52   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    5753   REAL(wp) ::   zzero  = 0._wp      ! 
    5854   REAL(wp) ::   zone   = 1._wp      ! 
    … …  
    126122               DO ji = 1, jpi 
    127123                  !Energy of melting q(S,T) [J.m-3] 
    128                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i ) 
     124                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 
    129125                  !0 if no ice and 1 if yes 
    130                   zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) )  )  
     126                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 )  ) 
    131127                  !convert units ! very important that this line is here 
    132128                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    138134               DO ji = 1, jpi 
    139135                  !Energy of melting q(S,T) [J.m-3] 
    140                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s ) 
     136                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 
    141137                  !0 if no ice and 1 if yes 
    142                   zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) )  )  
     138                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 )  ) 
    143139                  !convert units ! very important that this line is here 
    144140                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    147143         END DO 
    148144      END DO 
    149  
    150       !----------------------------- 
    151       ! 1.3) Set some dummies to 0 
    152       !----------------------------- 
    153       !clem rdvosif(:,:) = 0.e0   ! variation of ice volume at surface 
    154       !clem rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom 
    155       !clem fdvolif(:,:) = 0.e0   ! total variation of ice volume 
    156       !clem rdvonif(:,:) = 0.e0   ! lateral variation of ice volume 
    157       !clem fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice 
    158       !clem ffltbif(:,:) = 0.e0   ! linked with fstric 
    159       !clem qfvbq  (:,:) = 0.e0   ! linked with fstric 
    160       !clem rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area 
    161       !clem rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area 
    162       !clem hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.  
    163       !clem sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean 
    164       !clem fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean 
    165       !clem sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay 
    166145 
    167146      !----------------------------------- 
    … …  
    182161!CDIR NOVERRCHK 
    183162         DO ji = 1, jpi 
    184             phicif(ji,jj)  = vt_i(ji,jj) 
    185             pfrld(ji,jj)   = 1.0 - at_i(ji,jj) 
    186             zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 
     163            zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) ) 
    187164            ! 
    188165            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    … …  
    195172 
    196173            ! here the drag will depend on ice thickness and type (0.006) 
    197             fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )  
     174            fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )  
    198175            ! also category dependent 
    199176            !           !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead  
    200             qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 
     177            qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice 
    201178            !                        
    202179            !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh)  
    … …  
    214191            zpareff        = 1.0 - zinda * zfntlat 
    215192            != 0 if ice and positive heat budget and 1 if one of those two is false 
    216             zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 
     193            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) ) 
    217194            ! 
    218195            ! Heat budget of the lead, energy transferred from ice to ocean 
    … …  
    221198            ! 
    222199            ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 
    223             qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 
     200            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
    224201            ! 
    225202            ! oceanic heat flux (limthd_dh) 
    226             fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 
     203            fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) ) 
    227204            ! 
    228205         END DO 
    … …  
    240217         ENDIF 
    241218 
    242          zareamin = 1.e-10 
     219         zareamin = epsi10 
    243220         nbpb = 0 
    244221         DO jj = 1, jpj 
    … …  
    248225                  npb(nbpb) = (jj - 1) * jpi + ji 
    249226               ENDIF 
    250                ! debug point to follow 
    251                IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 
    252                   jiindex_1d = nbpb 
    253                ENDIF 
    254             END DO 
    255          END DO 
     227            END DO 
     228         END DO 
     229 
     230         ! debug point to follow 
     231         jiindex_1d = 0 
     232         IF( ln_nicep ) THEN 
     233            DO ji = mi0(jiindx), mi1(jiindx) 
     234               DO jj = mj0(jjindx), mj1(jjindx) 
     235                  jiindex_1d = (jj - 1) * jpi + ji 
     236               END DO 
     237            END DO 
     238         ENDIF 
    256239 
    257240         !------------------------------------------------------------------------------! 
    … …  
    315298            !-------------------------------- 
    316299 
    317             IF( con_i )   CALL lim_thd_enmelt( 1, nbpb )   ! computes sea ice energy of melting 
    318             IF( con_i )   CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 
     300            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_enmelt( 1, nbpb )   ! computes sea ice energy of melting 
     301            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 
    319302 
    320303            !                                 !---------------------------------! 
    … …  
    324307            CALL lim_thd_enmelt( 1, nbpb )    ! computes sea ice energy of melting compulsory for limthd_dh 
    325308 
    326             IF( con_i )   CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
    327             IF( con_i )   CALL lim_thd_con_dif( 1 , nbpb , jl ) 
     309            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
     310            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_con_dif( 1 , nbpb , jl ) 
    328311 
    329312            !                                 !---------------------------------! 
    … …  
    340323 
    341324            !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting 
    342             IF( con_i )   CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
    343             IF( con_i )   CALL lim_thd_con_dh ( 1 , nbpb , jl ) 
     325            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
     326            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_con_dh ( 1 , nbpb , jl ) 
    344327 
    345328            !-------------------------------- 
    … …  
    431414!clem@mv-to-itd    dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    432415 
    433       IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
     416      IF( con_i .AND. jiindex_1d > 0 )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
    434417 
    435418      IF(ln_ctl) THEN            ! Control print 
    … …  
    522505      END DO 
    523506      ! 
    524       IF(lwp) WRITE(numout,*) ' lim_thd_glohec ' 
    525       IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 
    526       IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 
    527       IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 
     507      WRITE(numout,*) ' lim_thd_glohec ' 
     508      WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 
     509      WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 
     510      WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 
    528511      ! 
    529512   END SUBROUTINE lim_thd_glohec 
    … …  
    611594      WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit 
    612595 
    613       IF (jiindex_1D.GT.0) WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d) 
    614       IF (jiindex_1D.GT.0) WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl) 
    615       IF (jiindex_1D.GT.0) WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d) 
     596      WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d) 
     597      WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl) 
     598      WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d) 
    616599 
    617600      !--------------------------------------- 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4161 r4333  
    3232   PUBLIC   lim_thd_dh   ! called by lim_thd 
    3333 
    34    REAL(wp) ::   epsi20 = 1e-20   ! constant values 
    35    REAL(wp) ::   epsi13 = 1e-13   ! 
    36    REAL(wp) ::   epsi16 = 1e-16   ! 
    37    REAL(wp) ::   zzero  = 0.e0    ! 
    38    REAL(wp) ::   zone   = 1.e0    ! 
     34   REAL(wp) ::   epsi20 = 1.e-20   ! constant values 
     35   REAL(wp) ::   epsi10 = 1.e-10   ! 
     36   REAL(wp) ::   epsi13 = 1.e-13   ! 
     37   REAL(wp) ::   zzero  = 0._wp    ! 
     38   REAL(wp) ::   zone   = 1._wp    ! 
    3939 
    4040   !!---------------------------------------------------------------------- 
    … …  
    125125      INTEGER  ::   num_iter_max, numce_dh 
    126126      REAL(wp) ::   meance_dh 
     127      REAL(wp) ::   zinda  
    127128      REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt 
    128129      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i 
    … …  
    288289      END DO 
    289290 
    290       !                     !------------------- 
    291       IF( con_i ) THEN      ! Conservation test 
    292          !                  !------------------- 
     291      !                                          !------------------- 
     292      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test 
     293         !                                       !------------------- 
    293294         numce_dh  = 0 
    294295         meance_dh = 0._wp 
    … …  
    494495      END DO ! jk 
    495496 
    496       !                     !------------------- 
    497       IF( con_i ) THEN      ! Conservation test 
    498       !                     !------------------- 
     497      !                                          !------------------- 
     498      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test 
     499      !                                          !------------------- 
    499500         DO ji = kideb, kiut 
    500501            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN 
    … …  
    544545         sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 
    545546         !                     ! excessive energy is sent to lateral ablation 
    546          fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 
     547         zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) ) 
     548         fsup(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice 
    547549      END DO 
    548550 
    … …  
    577579 
    578580         ! energy of melting of remaining snow 
    579          zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 
     581         zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) ) 
     582         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda 
    580583         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 
    581584         zhnfi     =  zhni + zdhnm 
    … …  
    613616         ENDIF 
    614617         ! 
    615          ! diagnostic ( bottom ice growth ) 
     618         ! diagnostic  
    616619         ii = MOD( npb(ji) - 1, jpi ) + 1 
    617620         ij = ( npb(ji) - 1 ) / jpi + 1 
    … …  
    689692         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
    690693         IF (  num_sal == 2  ) THEN 
    691             i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 
     694            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) 
    692695            ! salinity dif due to snow-ice formation 
    693             dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch      
     696            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
    694697            ! salinity dif due to bottom growth  
    695698            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
    696                dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 
     699               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch 
    697700            ENDIF 
    698701         ENDIF 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r4205 r4333  
    3131   PUBLIC   lim_thd_dif   ! called by lim_thd 
    3232 
    33    REAL(wp) ::   epsi20 = 1e-20     ! constant values 
    34    REAL(wp) ::   epsi13 = 1e-13     ! constant values 
    35  
     33   REAL(wp) ::   epsi10      =  1.e-10_wp    ! 
    3634   !!---------------------------------------------------------------------- 
    3735   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    … …  
    107105      INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation 
    108106      INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    109       REAL(wp) ::   zeps      =  1.e-10_wp    ! 
    110107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    111108      REAL(wp) ::   zg1       =  2._wp        ! 
    … …  
    114111      REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow 
    115112      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
    116       REAL(wp) ::   zht_smin  =  1.e-4_wp     ! minimum snow depth 
    117113      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    118114      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    … …  
    312308         IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula 
    313309            DO ji = kideb , kiut 
    314                ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-zeps,t_i_b(ji,1)-rtt) 
     310               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt) 
    315311               ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin) 
    316312            END DO 
    … …  
    318314               DO ji = kideb , kiut 
    319315                  ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  & 
    320                      MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
     316                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
    321317                  ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 
    322318               END DO 
    … …  
    326322         IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
    327323            DO ji = kideb , kiut 
    328                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -zeps, t_i_b(ji,1)-rtt )   & 
     324               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   & 
    329325                  &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )   
    330326               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
    … …  
    333329               DO ji = kideb , kiut 
    334330                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
    335                      &                                  / MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
     331                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
    336332                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )   
    337333                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 
    … …  
    339335            END DO 
    340336            DO ji = kideb , kiut 
    341                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-zeps,t_bo_b(ji)-rtt)   & 
     337               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   & 
    342338                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )   
    343339               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
    … …  
    352348 
    353349            !-- Snow kappa factors 
    354             zkappa_s(ji,0)         = rcdsn / MAX(zeps,zh_s(ji)) 
    355             zkappa_s(ji,nlay_s)    = rcdsn / MAX(zeps,zh_s(ji)) 
     350            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji)) 
     351            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji)) 
    356352         END DO 
    357353 
    … …  
    359355            DO ji = kideb , kiut 
    360356               zkappa_s(ji,layer)  = 2.0 * rcdsn / & 
    361                   MAX(zeps,2.0*zh_s(ji)) 
     357                  MAX(epsi10,2.0*zh_s(ji)) 
    362358            END DO 
    363359         END DO 
    … …  
    367363               !-- Ice kappa factors 
    368364               zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ & 
    369                   MAX(zeps,2.0*zh_i(ji))  
    370             END DO 
    371          END DO 
    372  
    373          DO ji = kideb , kiut 
    374             zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) 
    375             zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) 
     365                  MAX(epsi10,2.0*zh_i(ji))  
     366            END DO 
     367         END DO 
     368 
     369         DO ji = kideb , kiut 
     370            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 
     371            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 
    376372            !-- Interface 
    377             zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 
     373            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 
    378374               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    379375            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
    … …  
    389385               ztitemp(ji,layer)   = t_i_b(ji,layer) 
    390386               zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 
    391                   MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps) 
     387                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10) 
    392388               zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 
    393                   zeps) 
     389                  epsi10) 
    394390            END DO 
    395391         END DO 
    … …  
    398394            DO ji = kideb , kiut 
    399395               ztstemp(ji,layer) = t_s_b(ji,layer) 
    400                zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps) 
     396               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
    401397            END DO 
    402398         END DO 
    … …  
    707703      END DO  ! End of the do while iterative procedure 
    708704 
    709       IF( ln_nicep ) THEN 
     705      IF( ln_nicep .AND. lwp ) THEN 
    710706         WRITE(numout,*) ' zerritmax : ', zerritmax 
    711707         WRITE(numout,*) ' nconv     : ', nconv 
    … …  
    732728      ! Heat conservation       ! 
    733729      !-------------------------! 
    734       IF( con_i ) THEN 
     730      IF( con_i .AND. jiindex_1d > 0 ) THEN 
    735731         DO ji = kideb, kiut 
    736732            ! Upper snow value 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r4161 r4333  
    3636   PUBLIC   lim_thd_ent         ! called by lim_thd 
    3737 
    38    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    39    REAL(wp) ::   epsi13 = 1e-13_wp   ! 
    40    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    41    REAL(wp) ::   epsi06 = 1e-06_wp   ! 
     38   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
     39   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    4240   REAL(wp) ::   zzero  = 0._wp      ! 
    4341   REAL(wp) ::   zone   = 1._wp      ! 
    … …  
    122120      REAL(wp), POINTER, DIMENSION(:,:) ::   zhl0      ! old and new layer thicknesses  
    123121      REAL(wp), POINTER, DIMENSION(:,:) ::   zrl01 
     122 
     123      REAL(wp) ::   zinda  
    124124      !!------------------------------------------------------------------- 
    125125 
    … …  
    391391         DO layer1 = ntop1, nbot1 
    392392            DO ji = kideb, kiut 
    393                zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   & 
     393               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 
     394               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   & 
    394395                  &                 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10))  
    395396               q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0)   & 
    … …  
    407408      END DO 
    408409 
    409       IF ( con_i ) THEN 
     410      IF ( con_i .AND. jiindex_1d > 0 ) THEN 
    410411         DO ji = kideb, kiut 
    411412            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    … …  
    429430      DO jk = 1, nlay_s 
    430431         DO ji = kideb, kiut 
    431             q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , epsi20 ) 
     432            zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) )         
     433            q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 ) 
    432434         END DO !ji 
    433435      END DO !jk   
    … …  
    555557         DO ji = kideb, kiut 
    556558            ! Heat conservation 
    557             zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 
     559            zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) & 
    558560               &                                   * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 
    559561         END DO 
    … …  
    601603         DO layer1 = ntop1, nbot1 
    602604            DO ji = kideb, kiut 
    603                zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 
     605               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 
     606               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 
    604607                  - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 
    605608               q_i_b(ji,layer1) = q_i_b(ji,layer1) &  
    606609                  + zrl01(layer1,layer0)*qm0(ji,layer0) & 
    607                   * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 
     610                  * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) & 
    608611                  * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 
    609612            END DO 
    … …  
    621624      END DO 
    622625      ! 
    623       IF ( con_i ) THEN 
     626      IF ( con_i .AND. jiindex_1d > 0 ) THEN 
    624627         DO ji = kideb, kiut 
    625628            IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    … …  
    646649      DO jk = 1, nlay_i 
    647650         DO ji = kideb, kiut 
    648             q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , epsi20 ) 
     651            zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) ) 
     652            q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 ) 
    649653         END DO !ji 
    650654      END DO !jk   

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r4161 r4333  
    3636   PUBLIC lim_thd_lac     ! called by lim_thd 
    3737 
    38    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    39    REAL(wp) ::   epsi13 = 1e-13_wp   ! 
    40    REAL(wp) ::   epsi11 = 1e-11_wp   ! 
    41    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    42    REAL(wp) ::   epsi06 = 1e-06_wp   ! 
    43    REAL(wp) ::   epsi03 = 1e-03_wp   ! 
     38   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    4439   REAL(wp) ::   zzero  = 0._wp      ! 
    4540   REAL(wp) ::   zone   = 1._wp      ! 
    … …  
    160155                  !Energy of melting q(S,T) [J.m-3] 
    161156                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i ) 
    162                   zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes 
     157                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes 
    163158                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
    164159               END DO 
    … …  
    286281               nbpac = nbpac + 1 
    287282               npac( nbpac ) = (jj - 1) * jpi + ji 
    288                IF( ji == jiindx .AND. jj == jjindx )   jiindex_1d = nbpac 
    289283            ENDIF 
    290284         END DO 
    291285      END DO 
    292286 
    293       IF( ln_nicep )   WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
     287      ! debug point to follow 
     288      jiindex_1d = 0 
     289      IF( ln_nicep ) THEN 
     290         DO ji = mi0(jiindx), mi1(jiindx) 
     291            DO jj = mj0(jjindx), mj1(jjindx) 
     292               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN 
     293                  jiindex_1d = (jj - 1) * jpi + ji 
     294               ENDIF 
     295            END DO 
     296         END DO 
     297      ENDIF 
     298    
     299      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    294300 
    295301      !------------------------------ 
    … …  
    499505         END DO 
    500506 
    501          IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl) 
     507         IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 
    502508         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    503509            DO ji = 1, nbpac 
    504510               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 
    505                zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem 
    506                zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 
    507             END DO 
    508          END DO 
    509          IF( ln_nicep )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl) 
     511               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) )  ! clem 
     512               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 ) 
     513            END DO 
     514         END DO 
     515         IF( ln_nicep .AND. jiindex_1d > 0 )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 
    510516 
    511517         !--------------------------------- 
    … …  
    649655         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)  
    650656         IF( ln_nicep ) THEN 
    651             WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx) 
    652             WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx) 
    653             WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx) 
    654             WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx) 
     657            DO ji = mi0(jiindx), mi1(jiindx) 
     658               DO jj = mj0(jjindx), mj1(jjindx) 
     659                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj) 
     660                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj) 
     661                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj) 
     662                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj) 
     663               END DO 
     664            END DO 
    655665         ENDIF 
    656666         ! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r4259 r4333  
    3636   PUBLIC   lim_trp    ! called by ice_step 
    3737 
    38    REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values 
    39    REAL(wp)  ::   epsi03 = 1.e-03_wp   
    4038   REAL(wp)  ::   epsi10 = 1.e-10_wp   
    41    REAL(wp)  ::   epsi16 = 1.e-16_wp 
    42    REAL(wp)  ::   epsi20 = 1.e-20_wp 
    4339   REAL(wp)  ::   rzero  = 0._wp    
    4440   REAL(wp)  ::   rone   = 1._wp 
    … …  
    446442 
    447443                  ! Switches and dummy variables 
    448                   zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 ) 
    449                   zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 ) 
     444                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi10 ) 
     445                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi10 ) 
    450446                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 
    451447                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
    … …  
    458454                  ENDIF 
    459455 
    460                   zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp  ) * a_i(ji,jj,jl) 
     456                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp  ) * a_i(ji,jj,jl) 
    461457                  oa_i (ji,jj,jl)  = zindic * zage  
    462458 
    … …  
    498494               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    499495               ! 
    500                zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) ) 
    501                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness 
     496               zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) ) 
     497               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness 
    502498            END DO 
    503499         END DO 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4293 r4333  
    4848   PUBLIC   lim_update1   ! routine called by ice_step 
    4949 
    50       REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants 
    51       REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       - 
    52       REAL(wp)  ::   epsi03 = 1.e-03_wp   !    -       - 
    5350      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    54       REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       - 
    55       REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       - 
    5651      REAL(wp)  ::   rzero  = 0._wp       !    -       - 
    5752      REAL(wp)  ::   rone   = 1._wp       !    -       - 
    … …  
    108103      !------------------------------------------------------------------------------ 
    109104 
    110       !--------------------- 
    111       ! Ice dynamics   
    112       !--------------------- 
    113       u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:) 
    114       v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:) 
    115   
    116       !----------------------------- 
    117       ! Update ice and snow volumes   
    118       !----------------------------- 
    119       DO jl = 1, jpl 
    120          v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_trp(:,:,jl)  
    121          v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_trp(:,:,jl) 
    122       END DO 
    123  
    124  
    125       !--------------------------------------------- 
    126       ! Ice concentration and ice heat content 
    127       !--------------------------------------------- 
    128       a_i (:,:,:)  = a_i (:,:,:) + d_a_i_trp(:,:,:) 
    129       e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)   
    130  
    131       !------------------------------ 
    132       ! Snow temperature and ice age 
    133       !------------------------------ 
    134       e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:) 
    135       oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_trp(:,:,:)   
    136  
    137       !-------------- 
    138       ! Ice salinity     
    139       !-------------- 
    140       IF(  num_sal == 2  ) THEN      ! general case 
    141          smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:) 
    142       ENDIF 
     105      !----------------- 
     106      !  Trend terms 
     107      !----------------- 
     108      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
     109      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
     110      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
     111      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
     112      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
     113      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
     114      d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:) 
     115      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
     116      d_smv_i_trp(:,:,:)   = 0._wp 
     117      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    143118 
    144119      ! mass and salt flux init (clem) 
    … …  
    160135      CALL lim_var_glo2eqv 
    161136 
    162       !--------------------------------- 
    163       ! Classify the pathological cases 
    164       !--------------------------------- 
    165       ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case) 
    166       ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation) 
    167       ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell) 
    168       ! 
    169       DO jl = 1, jpl 
    170          DO jj = 1, jpj 
    171             DO ji = 1, jpi 
    172                patho_case(ji,jj,jl) = 1 
    173                IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN 
    174                   patho_case(ji,jj,jl) = 3 
    175                ENDIF 
    176                IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 
    177                   patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell 
    178                ENDIF 
    179             END DO 
    180          END DO 
    181       END DO 
    182  
    183137      !-------------------------------------- 
    184138      ! 2. Review of all pathological cases 
    185139      !-------------------------------------- 
    186140 
     141! clem: useless now 
    187142      !------------------------------------------- 
    188143      ! 2.1) Advection of ice in an ice-free cell 
    189144      !------------------------------------------- 
    190145      ! should be removed since it is treated after dynamics now 
    191  
    192 !      !IF( ln_nicep ) THEN   
    193 !         WRITE(numout,*) ' limupdate1 - before h correction ' 
    194 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    195 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    196 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    197 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    198 !      !ENDIF 
     146!      zhimax = 5._wp 
     147!      ! first category 
     148!      DO jj = 1, jpj 
     149!         DO ji = 1, jpi 
     150!            !--- the thickness of such an ice is often out of bounds 
     151!            !--- thus we recompute a new area while conserving ice volume 
     152!            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
     153!            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) )  
     154!            IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) & 
     155!              &   .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
     156!              &   .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
     157!               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
     158!               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
     159!            ENDIF 
     160!         END DO 
     161!      END DO 
    199162! 
    200       zhimax = 5._wp 
    201       ! first category 
    202       DO jj = 1, jpj 
    203          DO ji = 1, jpi 
    204             !--- the thickness of such an ice is often out of bounds 
    205             !--- thus we recompute a new area while conserving ice volume 
    206             zat_i_old = SUM( old_a_i(ji,jj,:) ) 
    207             zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) )  
    208             IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) & 
    209               &   .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
    210               &   .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
    211                ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
    212                a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
    213             ENDIF 
    214          END DO 
    215       END DO 
    216  
    217  
    218 !      !IF( ln_nicep ) THEN   
    219 !         at_i(:,:) = 0._wp 
    220 !         DO jl = 1, jpl 
    221 !            at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    222 !         END DO 
    223 !         WRITE(numout,*) ' limupdate1 - after h correction 1 ' 
    224 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    225 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    226 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    227 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    228 !      !ENDIF 
    229  
    230       zhimax = 20._wp 
    231       ! other categories 
    232       DO jl = 2, jpl 
    233          jm = ice_types(jl) 
    234          DO jj = 1, jpj 
    235             DO ji = 1, jpi 
    236                zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) )  
    237                ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
    238                ! it makes problems when the advected volume and concentration do not seem to be  
    239                ! related with each other 
    240                ! the new thickness is sometimes very big! 
    241                ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
    242                ! which of course is plausible 
    243                ! but fuck! it fucks everything up :) 
    244                IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
    245                   &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
    246                   ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
    247                   a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
    248                ENDIF 
    249             END DO ! ji 
    250          END DO !jj 
    251       END DO !jl 
     163!      zhimax = 20._wp 
     164!      ! other categories 
     165!      DO jl = 2, jpl 
     166!         jm = ice_types(jl) 
     167!         DO jj = 1, jpj 
     168!            DO ji = 1, jpi 
     169!               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) )  
     170!               ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
     171!               ! it makes problems when the advected volume and concentration do not seem to be  
     172!               ! related with each other 
     173!               ! the new thickness is sometimes very big! 
     174!               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
     175!               ! which of course is plausible 
     176!               ! but fuck! it fucks everything up :) 
     177!               IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
     178!                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
     179!                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
     180!                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
     181!               ENDIF 
     182!            END DO ! ji 
     183!         END DO !jj 
     184!      END DO !jl 
    252185      
    253186      at_i(:,:) = 0._wp 
    … …  
    255188         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    256189      END DO 
    257  
    258 !      !IF( ln_nicep ) THEN   
    259 !         WRITE(numout,*) ' limupdate1 - after h correction 2 ' 
    260 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    261 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    262 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    263 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    264 !      !ENDIF 
    265190 
    266191      !---------------------------------------------------- 
    … …  
    285210 
    286211               !switches 
    287                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
     212               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
    288213               !switch = 1 if a_i > 1e-06 and 0 if not 
    289                zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not 
    290                zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not 
     214               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 
     215               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 
    291216               ! bug fix 25 avril 2007 
    292217               zindb         = zindb*zindic 
    … …  
    332257               !-------------------------------------------- 
    333258               ! if greater than 0, ice concentration cannot be smaller than 1e-10 
    334                !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 ) 
    335259               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 
    336260 
    … …  
    352276            DO jj = 1, jpj  
    353277               DO ji = 1, jpi 
    354                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )  
     278                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    355279                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 
    356280               END DO ! ji 
    … …  
    370294         DO ji = 1, jpi 
    371295            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
    372             zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) )  
     296            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) )  
    373297            DO jl  = 1, jpl 
    374                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 
     298               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 
    375299               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
    376300               ! 
    377                zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
    378                ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 
     301               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
     302               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 
    379303               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
    380304            END DO 
    … …  
    412336                     !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 
    413337                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
    414                      i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) ) 
    415                      smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
     338                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
     339                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
    416340                  END DO ! ji 
    417341               END DO ! jj 
    … …  
    512436            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ') 
    513437            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ') 
    514             CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update1  : Path. case  : ') 
    515438 
    516439            DO jk = 1, nlay_i 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r4293 r4333  
    4545   PUBLIC   lim_update2   ! routine called by ice_step 
    4646 
    47       REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants 
    48       REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       - 
    4947      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    50       REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       - 
    51       REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       - 
    5248      REAL(wp)  ::   rzero  = 0._wp       !    -       - 
    5349      REAL(wp)  ::   rone   = 1._wp       !    -       - 
    … …  
    10298      CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem 
    10399 
    104       !------------------------------------------------------------------------------ 
    105       ! 1. Update of Global variables                                               | 
    106       !------------------------------------------------------------------------------ 
    107  
    108       !----------------------------- 
    109       ! Update ice and snow volumes   
    110       !----------------------------- 
    111       DO jl = 1, jpl 
    112          v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_thd(:,:,jl)  
    113          v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_thd(:,:,jl) 
    114       END DO 
    115   
    116       !--------------------------------------------- 
    117       ! Ice concentration and ice heat content 
    118       !--------------------------------------------- 
    119       a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:) 
    120       e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)   
    121   
    122       !------------------------------ 
    123       ! Snow temperature and ice age 
    124       !------------------------------ 
    125       e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:) 
    126       oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_thd(:,:,:) 
    127  
    128       !-------------- 
    129       ! Ice salinity     
    130       !-------------- 
    131       IF(  num_sal == 2  ) THEN      ! general case 
    132          smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) 
    133       ENDIF 
     100      !---------------------------------------------------------------------------------------- 
     101      ! 1. Computation of trend terms       
     102      !---------------------------------------------------------------------------------------- 
     103      !- Trend terms 
     104      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:)  
     105      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:) 
     106      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)   
     107      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
     108      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
     109      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:) 
     110      d_smv_i_thd(:,:,:) = 0._wp 
     111      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     112      ! diag only (clem) 
     113      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    134114 
    135115      ! mass and salt flux init (clem) 
    … …  
    151131      CALL lim_var_glo2eqv 
    152132 
    153       !--------------------------------- 
    154       ! Classify the pathological cases 
    155       !--------------------------------- 
    156       ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case) 
    157       ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation) 
    158       ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation) 
    159       ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation  
    160       ! with negative advection, very pathological ) 
    161       ! 
    162       DO jl = 1, jpl 
    163          DO jj = 1, jpj 
    164             DO ji = 1, jpi 
    165                patho_case(ji,jj,jl) = 1 
    166                IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN 
    167                   IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN  
    168                      patho_case(ji,jj,jl) = 2 
    169                   ENDIF 
    170                ELSE 
    171                   patho_case(ji,jj,jl) = 3 
    172                   IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN  
    173                      patho_case(ji,jj,jl) = 4 
    174                   ENDIF 
    175                ENDIF 
    176              END DO 
    177          END DO 
    178       END DO 
    179  
    180133      !-------------------------------------- 
    181134      ! 2. Review of all pathological cases 
    182135      !-------------------------------------- 
     136 
     137! clem: useless now 
    183138      !------------------------------------------- 
    184139      ! 2.1) Advection of ice in an ice-free cell 
    185140      !------------------------------------------- 
    186141      ! should be removed since it is treated after dynamics now 
    187  
    188 !      !IF( ln_nicep ) THEN   
    189 !         WRITE(numout,*) ' limupdate2 - before h correction ' 
    190 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    191 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    192 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    193 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    194 !      !ENDIF 
    195  
    196       zhimax = 5._wp 
    197       ! first category 
    198       DO jj = 1, jpj 
    199          DO ji = 1, jpi 
    200             !--- the thickness of such an ice is often out of bounds 
    201             !--- thus we recompute a new area while conserving ice volume 
    202             zat_i_old = SUM( old_a_i(ji,jj,:) ) 
    203             zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) )  
    204             IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) & 
    205                &  .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
    206                &  .AND. ( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
    207                ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
    208                a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
    209             ENDIF 
    210          END DO 
    211       END DO 
    212  
    213 !      !IF( ln_nicep ) THEN   
    214 !         at_i(:,:) = 0._wp 
    215 !         DO jl = 1, jpl 
    216 !            at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     142!      zhimax = 5._wp 
     143!      ! first category 
     144!      DO jj = 1, jpj 
     145!         DO ji = 1, jpi 
     146!            !--- the thickness of such an ice is often out of bounds 
     147!            !--- thus we recompute a new area while conserving ice volume 
     148!            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
     149!            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) )  
     150!            IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) & 
     151!               &  .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
     152!               &  .AND. ( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
     153!               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
     154!               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
     155!            ENDIF 
    217156!         END DO 
    218 !         WRITE(numout,*) ' limupdate2 - after h correction 1 ' 
    219 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    220 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    221 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    222 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    223 !      !ENDIF 
    224 !  
    225       zhimax = 20._wp 
    226       ! other categories 
    227       DO jl = 2, jpl 
    228          jm = ice_types(jl) 
    229          DO jj = 1, jpj 
    230             DO ji = 1, jpi 
    231                zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) )  
    232               ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
    233                ! it makes problems when the advected volume and concentration do not seem to be  
    234                ! related with each other 
    235                ! the new thickness is sometimes very big! 
    236                ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
    237                ! which of course is plausible 
    238                ! but fuck! it fucks everything up :) 
    239                IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
    240                   &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
    241                   ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
    242                   a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
    243                ENDIF 
    244             END DO ! ji 
    245          END DO !jj 
    246       END DO !jl 
     157!      END DO 
     158 
     159!      zhimax = 20._wp 
     160!      ! other categories 
     161!      DO jl = 2, jpl 
     162!         jm = ice_types(jl) 
     163!         DO jj = 1, jpj 
     164!            DO ji = 1, jpi 
     165!               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) )  
     166!              ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
     167!               ! it makes problems when the advected volume and concentration do not seem to be  
     168!               ! related with each other 
     169!               ! the new thickness is sometimes very big! 
     170!               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
     171!               ! which of course is plausible 
     172!               ! but fuck! it fucks everything up :) 
     173!               IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
     174!                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
     175!                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
     176!                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
     177!               ENDIF 
     178!            END DO ! ji 
     179!         END DO !jj 
     180!      END DO !jl 
    247181      
    248182      at_i(:,:) = 0._wp 
    … …  
    250184         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    251185      END DO 
    252  
    253 !      !IF( ln_nicep ) THEN   
    254 !         WRITE(numout,*) ' limupdate2 - after h correction 2 ' 
    255 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    256 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    257 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    258 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    259 !      !ENDIF 
    260186 
    261187      !---------------------------------------------------- 
    … …  
    304230                     zesum = zesum + e_i(ji,jj,jk,jl) 
    305231                  END DO 
    306                   !clem IF (ind_im .LT. nlay_i ) THEN 
    307                   !clem   smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) ) 
    308                   !clem ENDIF 
    309232                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 
    310233                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 
    … …  
    368291            DO ji = 1, jpi 
    369292               ! snow energy of melting 
    370                ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), 1.0e-6 )  ! snow energy of melting 
     293               zinda   =  MAX( 0._wp, SIGN( 1._wp, v_s(ji,jj,jl) - epsi10 ) ) 
     294               ze_s = zinda * e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), epsi10 )  ! snow energy of melting 
    371295 
    372296               ! If snow energy of melting smaller then Lf 
    … …  
    401325 
    402326               !switches 
    403                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
     327               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
    404328               !switch = 1 if a_i > 1e-06 and 0 if not 
    405                zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not 
    406                zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not 
     329               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 
     330               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 
    407331               ! bug fix 25 avril 2007 
    408332               zindb         = zindb*zindic 
    … …  
    444368               !--- 2.7 Correction to ice concentrations  
    445369               !-------------------------------------------- 
    446                !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 ) 
    447370               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 
    448371 
    … …  
    464387            DO jj = 1, jpj  
    465388               DO ji = 1, jpi 
    466                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )  
     389                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    467390                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 
    468391               END DO ! ji 
    … …  
    478401               !--- 2.12 Constrain the thickness of the smallest category above 5 cm 
    479402               !---------------------------------------------------------------------- 
    480                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
    481                ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06) 
    482                zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) ) 
     403               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
     404               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi10) 
     405               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi10 ) ) 
    483406               ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh 
    484407               ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh 
    … …  
    502425         DO ji = 1, jpi 
    503426            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
    504             zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) )  
     427            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) )  
    505428            DO jl  = 1, jpl 
    506                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 
     429               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 
    507430               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
    508431               ! 
    509                zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
    510                ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 
     432               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
     433               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 
    511434               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
    512435            END DO 
    … …  
    542465                     !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 
    543466                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
    544                      i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) ) 
    545                      smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
     467                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
     468                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
    546469                  END DO ! ji 
    547470               END DO ! jj 
    … …  
    663586            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ') 
    664587            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ') 
    665             CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update2  : Path. case  : ') 
    666588 
    667589            DO jk = 1, nlay_i 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r4205 r4333  
    6666   PUBLIC   lim_var_salprof1d    ! 
    6767 
    68    REAL(wp) ::   epsi20 = 1.e-20_wp   ! module constants 
    69    REAL(wp) ::   epsi16 = 1.e-16_wp   !    -       - 
    70    REAL(wp) ::   epsi13 = 1.e-13_wp   !    -       - 
    7168   REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       - 
    72    REAL(wp) ::   epsi06 = 1.e-06_wp   !    -       - 
    7369   REAL(wp) ::   zzero = 0.e0        !    -       - 
    7470   REAL(wp) ::   zone  = 1.e0        !    -       - 
    … …  
    117113               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    118114               ! 
    119                zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) )  
    120                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness 
     115               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) )  
     116               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness 
    121117            END DO 
    122118         END DO 
    … …  
    138134            DO jj = 1, jpj 
    139135               DO ji = 1, jpi 
    140                   zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi16 ) )  
    141                   zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) )  
     136                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi10 ) )  
     137                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) )  
    142138                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content 
    143                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi16 ) * zinda   ! ice salinity 
    144                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi16 ) * zindb   ! ice age 
     139                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda   ! ice salinity 
     140                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * zindb   ! ice age 
    145141               END DO 
    146142            END DO 
    … …  
    209205               DO ji = 1, jpi 
    210206                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    211                   zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp)  
    212                   zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes 
     207                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)  
     208                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes 
    213209                  zq_i    = zq_i * unit_fac * zindb                              !convert units 
    214210                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature 
    … …  
    235231               DO ji = 1, jpi 
    236232                  !Energy of melting q(S,T) [J.m-3] 
    237                   zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp) 
    238                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes 
     233                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 
     234                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes 
    239235                  zq_s  = zq_s * unit_fac * zindb                                    ! convert units 
    240236                  ! 
    … …  
    341337                  zind01 = ( 1._wp - zind0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) )  
    342338                  ! If 2.sm_i GE sss_m then zindbal = 1 
     339                  ! this is to force a constant salinity profile in the Baltic Sea 
    343340                  zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
    344341                  zalpha(ji,jj,jl) = zind0  + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 
    … …  
    434431            DO jj = 1, jpj 
    435432               DO ji = 1, jpi 
    436                   zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi16 ) )  ) 
    437                   zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi16 ) )  ) 
    438                   zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi16 )   & 
     433                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  ) 
     434                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
     435                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   & 
    439436                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp) 
    440                   bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi16 ) 
     437                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi10 ) 
    441438               END DO 
    442439            END DO 
    … …  
    498495               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) )  
    499496               ! if 2.sm_i GE sss_m then zindbal = 1 
     497               ! this is to force a constant salinity profile in the Baltic Sea 
    500498               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) ) 
    501499               ! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r4161 r4333  
    5252   INTEGER            , DIMENSION(jpnoumax) ::   nc  , nca     ! switch for saving field ( = 1 ) or not ( = 0 ) 
    5353 
    54    REAL(wp)  ::   epsi06 = 1e-6_wp 
     54   REAL(wp)  ::   epsi06 = 1.e-6_wp 
    5555   REAL(wp)  ::   zzero  = 0._wp 
    5656   REAL(wp)  ::   zone   = 1._wp       
    … …  
    192192 
    193193      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    194       IF( ln_nicep ) THEN 
    195          WRITE(numout,*) 
    196          WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit 
    197          WRITE(numout,*) '~~~~~~~ ' 
    198          WRITE(numout,*) ' kindic = ', kindic 
    199       ENDIF 
     194      !IF( ln_nicep ) THEN 
     195      !   WRITE(numout,*) 
     196      !   WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit 
     197      !   WRITE(numout,*) '~~~~~~~ ' 
     198      !   WRITE(numout,*) ' kindic = ', kindic 
     199      !ENDIF 
    200200      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    201201 
    … …  
    570570      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    571571      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     572      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     573      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     574      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     575      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     576      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     577      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     578      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     579      CALL histdef( kid, "iisfxthd", "Salt flux from thermo"   , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     580      CALL histdef( kid, "iisfxmec", "Salt flux from dynmics"  , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     581      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     582 
    572583 
    573584      !CALL histdef( kid, "iice_itd", "Ice concentration by cat", "%"      , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt )  
    … …  
    592603      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    593604 
     605      CALL histwrite( kid, "iicebopr", kt, diag_bot_gr        , jpi*jpj, (/1/) ) 
     606      CALL histwrite( kid, "iicedypr", kt, diag_dyn_gr        , jpi*jpj, (/1/) ) 
     607      CALL histwrite( kid, "iicelapr", kt, diag_lat_gr        , jpi*jpj, (/1/) ) 
     608      CALL histwrite( kid, "iicesipr", kt, diag_sni_gr        , jpi*jpj, (/1/) ) 
     609      CALL histwrite( kid, "iicerepr", kt, diag_res_pr        , jpi*jpj, (/1/) ) 
     610      CALL histwrite( kid, "iicebome", kt, diag_bot_me        , jpi*jpj, (/1/) ) 
     611      CALL histwrite( kid, "iicesume", kt, diag_sur_me        , jpi*jpj, (/1/) ) 
     612      CALL histwrite( kid, "iisfxthd", kt, sfx_thd        , jpi*jpj, (/1/) ) 
     613      CALL histwrite( kid, "iisfxmec", kt, sfx_mec        , jpi*jpj, (/1/) ) 
     614      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) ) 
     615 
    594616      !CALL histwrite( kid, "iice_itd", kt, a_i  , jpi*jpj*jpl, (/1/)  )   ! area 
    595617      !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/)  )   ! thickness 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r4292 r4333  
    104104   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points 
    105105 
    106 #if defined key_lim2 
    107    CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim2      ! Choice of boundary condition for sea ice variables  
    108    INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim2_dta  !: = 0 use the initial state as bdy dta ;  
     106#if ( defined key_lim2 || defined key_lim3 ) 
     107   CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim       ! Choice of boundary condition for sea ice variables  
     108   INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim_dta   !: = 0 use the initial state as bdy dta ;  
    109109                                                              !: = 1 read it in a NetCDF file 
    110110#endif 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4292 r4333  
    197197 
    198198#if defined key_lim2 
    199             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
     199            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN  
    200200               ilen1(:) = nblen(:) 
    201201               IF( dta%ll_frld ) THEN 
    … …  
    649649#if defined key_lim2 
    650650            ! sea ice 
    651             IF( nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
     651            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 
    652652 
    653653               IF( dta%ll_frld ) THEN 
    … …  
    724724               ENDIF 
    725725 
     726            ENDIF 
    726727#endif 
    727728            ! Recalculate field counts 
    … …  
    850851#if defined key_lim2 
    851852         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
    852             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
     853            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
    853854               ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) ) 
    854855               ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) ) 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r4292 r4333  
    4545 
    4646   REAL(wp) ::   epsi20 = 1.e-20_wp  ! module constants 
     47   REAL(wp) ::   epsi10 = 1.e-10_wp  ! min area allowed by ice model 
    4748   !!---------------------------------------------------------------------- 
    4849   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    9192      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index      !! 
    9293 
     94      INTEGER  ::   jpbound            ! 0 = incoming ice 
     95                                       ! 1 = outgoing ice 
    9396      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices 
    94       INTEGER  ::   ji, jj             ! local scalar 
     97      INTEGER  ::   ji, jj, ii, ij     ! local scalar 
    9598      REAL(wp) ::   zwgt, zwgt1        ! local scalar 
    96       REAL(wp) ::   zinda, ztmelts 
     99      REAL(wp) ::   zinda, ztmelts, zdh 
     100 
     101      REAL(wp), PARAMETER  ::   zsal = 6.3    ! arbitrary salinity    for incoming ice 
     102      REAL(wp), PARAMETER  ::   ztem = 270.0  ! arbitrary temperature for incoming ice 
     103      REAL(wp), PARAMETER  ::   zage = 30.0   ! arbitrary age         for incoming ice 
    97104      !!------------------------------------------------------------------------------ 
    98105      ! 
    99106      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs') 
    100107      ! 
    101       !IF( kt /= nit000 ) THEN 
    102108      jgrd = 1      ! Everything is at T-points here 
    103109      ! 
    … …  
    170176            ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth  
    171177            ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth 
     178 
     179            ! ----------------- 
     180            ! Pathological case 
     181            ! ----------------- 
     182            ! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to  
     183            ! very large transformation from snow to ice (see limthd_dh.F90) 
     184 
     185            ! Then, a) transfer the snow excess into the ice (different from limthd_dh) 
     186            zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
     187            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment) 
     188            !zdh = MAX( 0._wp, ht_s(ji,jj,jl) * rhosn / rhoic ) 
     189 
     190            ! recompute ht_i, ht_s 
     191            ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
     192            ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn )  
     193 
    172194         ENDDO 
    173195         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy ) 
    174196         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
    175197         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy ) 
    176  
     198      ENDDO 
     199      ! retrieve at_i 
     200      at_i(:,:) = 0._wp 
     201      DO jl = 1, jpl 
     202         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     203      END DO 
     204 
     205      DO jl = 1, jpl 
    177206         DO jb = 1, idx%nblen(jgrd) 
    178207            ji    = idx%nbi(jb,jgrd) 
    179208            jj    = idx%nbj(jb,jgrd) 
    180             ! clem : condition on ice thickness depends on the ice velocity 
    181             ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values  
    182             IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) THEN 
    183                ht_i(ji,jj,jl) = ht_i(ji+1,jj  ,jl) 
    184                a_i (ji,jj,jl) = a_i (ji+1,jj  ,jl) 
    185                ht_s(ji,jj,jl) = ht_s(ji+1,jj  ,jl) 
    186             ENDIF 
    187             IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) THEN 
    188                ht_i(ji,jj,jl) = ht_i(ji-1,jj  ,jl) 
    189                a_i (ji,jj,jl) = a_i (ji-1,jj  ,jl) 
    190                ht_s(ji,jj,jl) = ht_s(ji-1,jj  ,jl) 
    191             ENDIF 
    192             IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) THEN 
    193                ht_i(ji,jj,jl) = ht_i(ji  ,jj+1,jl) 
    194                a_i (ji,jj,jl) = a_i (ji  ,jj+1,jl) 
    195                ht_s(ji,jj,jl) = ht_s(ji  ,jj+1,jl) 
    196             ENDIF 
    197             IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) THEN 
    198                ht_i(ji,jj,jl) = ht_i(ji  ,jj-1,jl) 
    199                a_i (ji,jj,jl) = a_i (ji  ,jj-1,jl) 
    200                ht_s(ji,jj,jl) = ht_s(ji  ,jj-1,jl) 
    201             ENDIF 
    202  
    203             zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - a_i(ji,jj,jl) ) ) ! 0 if no ice 
    204             ! -------------------- 
     209 
     210            ! condition on ice thickness depends on the ice velocity 
     211            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 
     212            jpbound = 0; ii = ji; ij = jj; 
     213 
     214            IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj 
     215            IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj 
     216            IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1 
     217            IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1 
     218 
     219            zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ii,ij) + 0.01 ) ) ! 0 if no ice 
     220 
     221            ! concentration and thickness 
     222            a_i (ji,jj,jl) = a_i (ii,ij,jl) * zinda 
     223            ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * zinda 
     224            ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * zinda 
     225 
    205226            ! Ice and snow volumes 
    206             ! -------------------- 
    207227            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 
    208228            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) 
    209229 
    210             ! ------------- 
    211             ! Ice salinity 
    212             !--------------- 
    213             sm_i(ji,jj,jl)   = zinda * bulk_sal  + ( 1.0 - zinda ) * s_i_min 
     230            SELECT CASE( jpbound ) 
     231 
     232            CASE( 0 ) ! velocity is inward 
     233 
     234               ! Ice salinity, age, temperature 
     235               sm_i(ji,jj,jl)   = zinda * zsal  + ( 1.0 - zinda ) * s_i_min 
     236               o_i(ji,jj,jl)    = zinda * zage  + ( 1.0 - zinda ) 
     237               t_su(ji,jj,jl)   = zinda * ztem  + ( 1.0 - zinda ) * ztem 
     238               DO jk = 1, nlay_s 
     239                  t_s(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt 
     240               END DO 
     241               DO jk = 1, nlay_i 
     242                  t_i(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt  
     243                  s_i(ji,jj,jk,jl) = zinda * zsal + ( 1.0 - zinda ) * s_i_min 
     244               END DO 
     245                
     246            CASE( 1 ) ! velocity is outward 
     247  
     248               ! Ice salinity, age, temperature 
     249               sm_i(ji,jj,jl)   = zinda * sm_i(ii,ij,jl)  + ( 1.0 - zinda ) * s_i_min 
     250               o_i(ji,jj,jl)    = zinda * o_i(ii,ij,jl)   + ( 1.0 - zinda ) 
     251               t_su(ji,jj,jl)   = zinda * t_su(ii,ij,jl)  + ( 1.0 - zinda ) * rtt 
     252               DO jk = 1, nlay_s 
     253                  t_s(ji,jj,jk,jl) = zinda * t_s(ii,ij,jk,jl) + ( 1.0 - zinda ) * rtt 
     254               END DO 
     255               DO jk = 1, nlay_i 
     256                  t_i(ji,jj,jk,jl) = zinda * t_i(ii,ij,jk,jl) + ( 1.0 - zinda ) * rtt 
     257                  s_i(ji,jj,jk,jl) = zinda * s_i(ii,ij,jk,jl) + ( 1.0 - zinda ) * s_i_min 
     258               END DO 
     259 
     260            END SELECT 
     261 
     262            ! contents 
    214263            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    215              
    216             !---------- 
    217             ! Ice age 
    218             !---------- 
    219             o_i(ji,jj,jl)    = zinda * 1.0   + ( 1.0 - zinda ) 
    220264            oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    221              
    222             !------------------------------ 
    223             ! Sea ice surface temperature 
    224             !------------------------------ 
    225             t_su(ji,jj,jl)   = zinda * 270.0 + ( 1.0 - zinda ) * t_bo(ji,jj) 
    226              
    227             !------------------------------------ 
    228             ! Snow temperature and heat content 
    229             !------------------------------------ 
    230265            DO jk = 1, nlay_s 
    231                t_s(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt 
    232266               ! Snow energy of melting 
    233267               e_s(ji,jj,jk,jl) = zinda * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    … …  
    235269               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    236270               ! Multiply by volume, so that heat content in 10^9 Joules 
    237                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    238                     v_s(ji,jj,jl)  / nlay_s 
    239             END DO !jk 
    240              
    241             !----------------------------------------------- 
    242             ! Ice salinities, temperature and heat content  
    243             !----------------------------------------------- 
     271               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
     272            END DO 
    244273            DO jk = 1, nlay_i 
    245                t_i(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt  
    246                s_i(ji,jj,jk,jl) = zinda * bulk_sal + ( 1.0 - zinda ) * s_i_min 
    247                ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    248                 
     274               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K                   
    249275               ! heat content per unit volume 
    250276               e_i(ji,jj,jk,jl) = zinda * rhoic * & 
    251                     (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    252                     +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
    253                     - rcp      * ( ztmelts - rtt ) ) 
    254                 
     277                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     278                  +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
     279                  - rcp      * ( ztmelts - rtt ) ) 
    255280               ! Correct dimensions to avoid big values 
    256281               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    257                 
    258282               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    259                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    260                     area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 
    261             END DO ! jk 
     283               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 
     284            END DO 
     285 
    262286 
    263287         END DO !jb 
    264288  
    265  
    266289         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions 
    267290         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
    … …  
    288311#endif 
    289312      !       
    290       !ENDIF !nit000/=0 
    291313      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs') 
    292314      ! 
    … …  
    294316 
    295317 
    296    SUBROUTINE bdy_ice_lim_dyn( kn, pvar1, pvar2, pvar12 ) 
     318   SUBROUTINE bdy_ice_lim_dyn( cd_type ) 
    297319      !!------------------------------------------------------------------------------ 
    298320      !!                 ***  SUBROUTINE bdy_ice_lim_dyn  *** 
    299321      !!                     
    300322      !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries. 
    301       !!      kn = 1: u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free 
     323      !!              u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free 
    302324      !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities 
    303       !!      kn = 2: the stress tensor is set to 0 (i.e. pvar1, pvar2, pvar12) 
    304325      !! 
    305326      !! 2013-06 : C. Rousset 
    306327      !!------------------------------------------------------------------------------ 
    307328      !! 
    308       INTEGER,  INTENT( in    )                           ::   kn     ! set up of ice vel (kn=1) or stress tensor (kn=2) 
    309       REAL(wp), INTENT( inout ), DIMENSION(:,:), OPTIONAL ::   pvar1, pvar2, pvar12  ! stress tensor components (zs1,zs2,zs12)  
     329      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points 
    310330      INTEGER  ::   jb, jgrd   ! dummy loop indices 
    311331      INTEGER  ::   ji, jj             ! local scalar 
    … …  
    317337      ! 
    318338      DO ib_bdy=1, nb_bdy 
    319       ! 
    320       SELECT CASE( nn_ice_lim(ib_bdy) ) 
    321  
    322       CASE(jp_none) 
    323          CYCLE 
    324  
    325       CASE(jp_frs) 
    326  
    327          IF ( kn == 1 ) THEN ! set up of u_ice and v_ice 
    328             
    329             jgrd = 2      ! u velocity 
    330             DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
    331                ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
    332                jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    333                zflag = idx_bdy(ib_bdy)%flagu(jb) 
    334  
    335                IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries 
    336                   ! one of the two zmsk is always 0 (because of zflag) 
    337                   zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 
    338                   zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 
     339         ! 
     340         SELECT CASE( nn_ice_lim(ib_bdy) ) 
     341 
     342         CASE('none') 
     343 
     344            CYCLE 
     345             
     346         CASE('frs') 
     347             
     348 
     349            SELECT CASE ( cd_type ) 
     350 
     351            CASE ( 'U' ) 
     352                
     353               jgrd = 2      ! u velocity 
     354               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
     355                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
     356                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
     357                  zflag = idx_bdy(ib_bdy)%flagu(jb) 
    339358                   
    340                   ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    341                   u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
    342                      &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
    343                      &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    344                ELSE                             ! everywhere else 
    345                   u_ice(ji,jj) = u_oce(ji,jj) 
    346                ENDIF 
    347                ! mask ice velocities 
    348                zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice 
    349                u_ice(ji,jj) = zinda * u_ice(ji,jj) 
    350  
    351             ENDDO 
     359                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries 
     360                     ! one of the two zmsk is always 0 (because of zflag) 
     361                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 
     362                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 
     363                      
     364                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
     365                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
     366                        &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
     367                        &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
     368                  ELSE                             ! everywhere else 
     369                     !u_ice(ji,jj) = u_oce(ji,jj) 
     370                     u_ice(ji,jj) = 0._wp 
     371                  ENDIF 
     372                  ! mask ice velocities 
     373                  zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice 
     374                  u_ice(ji,jj) = zinda * u_ice(ji,jj) 
     375                   
     376               ENDDO 
     377 
     378               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy ) 
     379                
     380            CASE ( 'V' ) 
     381                
     382               jgrd = 3      ! v velocity 
     383               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
     384                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
     385                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
     386                  zflag = idx_bdy(ib_bdy)%flagv(jb) 
     387                   
     388                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries 
     389                     ! one of the two zmsk is always 0 (because of zflag) 
     390                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice 
     391                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice 
     392                      
     393                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
     394                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
     395                        &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
     396                        &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
     397                  ELSE                             ! everywhere else 
     398                     !v_ice(ji,jj) = v_oce(ji,jj) 
     399                     v_ice(ji,jj) = 0._wp 
     400                  ENDIF 
     401                  ! mask ice velocities 
     402                  zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice 
     403                  v_ice(ji,jj) = zinda * v_ice(ji,jj) 
     404                   
     405               ENDDO 
     406                
     407               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy ) 
     408                
     409            END SELECT 
    352410             
    353             jgrd = 3      ! v velocity 
    354             DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
    355                ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
    356                jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    357                zflag = idx_bdy(ib_bdy)%flagv(jb) 
    358                
    359                IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries 
    360                   ! one of the two zmsk is always 0 (because of zflag) 
    361                   zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice 
    362                   zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice 
    363                    
    364                   ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    365                   v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
    366                      &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
    367                      &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    368                ELSE                             ! everywhere else 
    369                   v_ice(ji,jj) = v_oce(ji,jj) 
    370                ENDIF 
    371                ! mask ice velocities 
    372                zinda = 1.0 - MAX( 0.0_wp, SIGN( 1.0_wp , - at_i(ji,jj) ) ) ! 0 if no ice 
    373                v_ice(ji,jj) = zinda * v_ice(ji,jj) 
    374  
    375             ENDDO 
    376              
    377          CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy ) 
    378          CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy ) 
    379   
    380          ELSEIF ( kn == 2 ) THEN ! set up of stress tensor (not sure it works) 
     411         CASE DEFAULT 
     412            CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' ) 
     413         END SELECT 
    381414          
    382             jgrd = 1      ! T grid 
    383             DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
    384                ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
    385                jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    386  
    387                pvar1 (ji,jj) = 0._wp 
    388                pvar2 (ji,jj) = 0._wp 
    389                pvar12(ji,jj) = 0._wp 
    390  
    391             ENDDO 
    392  
    393          ENDIF 
    394  
    395       CASE DEFAULT 
    396          CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' ) 
    397       END SELECT 
    398  
    399415      ENDDO 
    400416 
    401417      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn') 
    402   
     418       
    403419    END SUBROUTINE bdy_ice_lim_dyn 
    404420 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r4294 r4333  
    488488            DO igrd = 1, jpbgrd 
    489489               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    490                nblendta(igrd,ib_bdy) = kdimsz(1) 
    491                jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     490               !clem nblendta(igrd,ib_bdy) = kdimsz(1) 
     491               !clem jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     492               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 
     493               jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz)) 
    492494            ENDDO 
    493495            CALL iom_close( inum ) 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r4328 r4333  
    3636   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets 
    3737 
    38    REAL(dp), SAVE                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends 
    39    REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ssh_ini              ! 
    40    REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
    41    REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcssh_loc_ini, scssh_loc_ini     ! 
     38   REAL(wp), SAVE                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends 
     39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ssh_ini              ! 
     40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcssh_loc_ini, scssh_loc_ini     ! 
    4242 
    4343   !! * Substitutions 
    … …  
    6767      !! 
    6868      INTEGER    ::   jk                          ! dummy loop indice 
    69       REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations 
    70       REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation 
    71       REAL(dp)   ::   z_hc        , z_sc          ! heat and salt content 
    72       REAL(dp)   ::   z_v1        , z_v2          ! volume 
    73       REAL(dp)   ::   zdeltat                     !    -     - 
    74       REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     - 
    75       REAL(dp)   ::   z_frc_trd_v                 !    -     - 
    76       REAL(dp), POINTER, DIMENSION(:,:)   ::   zsurf              ! 
     69      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations 
     70      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation 
     71      REAL(wp)   ::   z_hc        , z_sc          ! heat and salt content 
     72      REAL(wp)   ::   z_v1        , z_v2          ! volume 
     73      REAL(wp)   ::   zdeltat                     !    -     - 
     74      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     - 
     75      REAL(wp)   ::   z_frc_trd_v                 !    -     - 
     76      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsurf              ! 
    7777      !!--------------------------------------------------------------------------- 
    7878      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')       
    … …  
    104104      ! 2a -  Content variations ! 
    105105      ! ------------------------ ! 
    106       zdiff_v2 = 0._dp 
    107       zdiff_hc = 0._dp 
    108       zdiff_sc = 0._dp 
     106      zdiff_v2 = 0._wp 
     107      zdiff_hc = 0._wp 
     108      zdiff_sc = 0._wp 
    109109      ! volume variation (calculated with ssh) 
    110110      zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) ) 
    111111      DO jk = 1, jpkm1 
    112112         ! volume variation (calculated with scale factors) 
    113          zdiff_v2 = zdiff_v2 & 
    114             &  + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 
     113         zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 
    115114         ! heat content variation 
    116          zdiff_hc = zdiff_hc & 
    117             &  + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   & 
     115         zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   & 
    118116            &                           - hc_loc_ini(:,:,jk) ) ) 
    119117         ! salt content variation 
    120          zdiff_sc = zdiff_sc & 
    121             &  + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   & 
     118         zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   & 
    122119            &                           - sc_loc_ini(:,:,jk) ) ) 
    123120      ENDDO 
    … …  
    140137      ! 2b -  Content           ! 
    141138      ! ----------------------- ! 
    142       z_v2 = 0._dp 
    143       z_hc = 0._dp 
    144       z_sc = 0._dp 
     139      z_v2 = 0._wp 
     140      z_hc = 0._wp 
     141      z_sc = 0._wp 
    145142      ! volume (calculated with ssh) 
    146143      z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) ) 
    147144      DO jk = 1, jpkm1 
    148145         ! volume (calculated with scale factors) 
    149          z_v2 = z_v2 & 
    150             &     + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
     146         z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
    151147         ! heat content 
    152          z_hc = z_hc & 
    153             &     + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) ) 
     148         z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) ) 
    154149         ! salt content 
    155          z_sc = z_sc & 
    156             &     + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) ) 
     150         z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) ) 
    157151      ENDDO 
    158152      ! add ssh if not vvl 
    … …  
    170164      CALL iom_put( 'bgtemper' , z_hc / z_v2 )                      ! Temperature (C)  
    171165      CALL iom_put( 'bgsaline' , z_sc / z_v2 )                      ! Salinity (psu) 
    172       CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_dp ) ! Heat content variation (10^9 J) 
     166      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J) 
    173167      CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9 )                 ! Salt content variation (psu*km3)  
    174168      CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 )                    ! volume ssh (km3)   
    … …  
    176170      CALL iom_put( 'bgvoltot' , zdiff_v2 * 1.e-9 )                 ! volume total (km3)  
    177171      CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 )                     ! vol - surface forcing (volume)  
    178       CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_dp ) ! hc  - surface forcing (heat content)  
     172      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content)  
    179173      CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 )                     ! sc  - surface forcing (salt content)  
    180174      ! 
    … …  
    224218      IF(lwp) THEN                   ! Control print 
    225219         WRITE(numout,*) 
     220         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
     221         WRITE(numout,*) '~~~~~~~~~~~~' 
    226222         WRITE(numout,*) '   Namelist namhsb : set hsb parameters' 
    227223         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb 
    … …  
    308304          hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
    309305          scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
    310           frc_v = 0._dp                                            
    311           frc_t = 0._dp                                            
    312           frc_s = 0._dp                                                   
     306          frc_v = 0._wp                                            
     307          frc_t = 0._wp                                            
     308          frc_s = 0._wp                                                   
    313309       ENDIF 
    314310 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r4313 r4333  
    2323   USE eosbn2          ! equation of state            (eos bn2 routine) 
    2424   USE trdmld_oce      ! ocean active mixed layer tracers trends variables 
     25   USE domvvl          ! variable volume 
    2526   USE divcur          ! hor. divergence and curl      (div & cur routines) 
    2627   USE sbc_ice, ONLY : lk_lim3 
    … …  
    210211         CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   ) 
    211212         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    ) 
     213         IF( lk_vvl )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    212214      ELSE 
    213215         neuler = 0 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r4306 r4333  
    793793        Cd_n10(:,:) =   cdn_wave 
    794794      ELSE 
    795         Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
     795        Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
    796796      ENDIF 
    797797      sqrt_Cd_n10 = sqrt(Cd_n10) 
    798       Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b) 
    799       Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) !   L & Y eq. (6c), (6d) 
     798      Ce_n10  = 1.e-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b) 
     799      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) !   L & Y eq. (6c), (6d) 
    800800      !! 
    801801      !! Initializing transfert coefficients with their first guess neutral equivalents : 
    … …  
    824824 
    825825           !! Updating the neutral 10m transfer coefficients : 
    826            Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
     826           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
    827827           sqrt_Cd_n10 = sqrt(Cd_n10) 
    828            Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
     828           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
    829829           stab    = 0.5 + sign(0.5,zeta) 
    830            Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d) 
     830           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab))           !  L & Y eq. (6c), (6d) 
    831831 
    832832           !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) : 
    … …  
    927927        Cd_n10(:,:) =   cdn_wave 
    928928      ELSE 
    929         Cd_n10  = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
     929        Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
    930930      ENDIF 
    931931      sqrt_Cd_n10 = sqrt(Cd_n10) 
    932       Ce_n10  = 1E-3*( 34.6 * sqrt_Cd_n10 ) 
    933       Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab)) 
     932      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
     933      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    934934 
    935935      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    … …  
    973973         ELSE 
    974974           !! Updating the neutral 10m transfer coefficients : 
    975            Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a) 
     975           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a) 
    976976           sqrt_Cd_n10 = sqrt(Cd_n10) 
    977            Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b) 
     977           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b) 
    978978           stab    = 0.5 + sign(0.5,zeta_u) 
    979            Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d) 
     979           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c-6d) 
    980980           !! 
    981981           !! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4269 r4333  
    5858   USE in_out_manager  ! I/O manager 
    5959   USE prtctl          ! Print control 
     60   USE lib_fortran     !  
    6061 
    6162#if defined key_bdy  
    … …  
    167168         ! 
    168169         IF( ln_nicep ) THEN      ! control print at a given point 
    169             jiindx = 6   ;   jjindx = 47 
    170             WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
     170            jiindx = 177   ;   jjindx = 112 
     171            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    171172         ENDIF 
    172173      ENDIF 
    … …  
    303304         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp 
    304305         fhmec  (:,:) = 0._wp   ;    
    305          fmmec  (:,:) = 0._wp      
     306         fmmec  (:,:) = 0._wp 
     307         fmmflx (:,:) = 0._wp      
    306308         focea2D(:,:) = 0._wp 
    307309         fsup2D (:,:) = 0._wp 
    … …  
    328330                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
    329331         ! 
    330          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    331          ! 
    332          IF( .NOT. lk_c1d ) THEN                     ! Ice dynamics & transport (except in 1D case) 
     332         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
     333         ! ---------------------------------------------- 
     334         ! ice dynamics and transport (except in 1D case) 
     335         ! ---------------------------------------------- 
     336         IF( .NOT. lk_c1d ) THEN 
    333337                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    334338                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    335339                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    336          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
     340         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
    337341                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    338342                          CALL lim_var_agg( 1 )  
     343#if defined key_bdy 
     344                          ! bdy ice thermo  
     345                          CALL lim_var_glo2eqv            ! equivalent variables 
     346                          CALL bdy_ice_lim( kt ) 
     347                          CALL lim_itd_me_zapsmall 
     348                          CALL lim_var_agg(1) 
     349         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
     350#endif 
    339351                          CALL lim_update1 
    340352         ENDIF 
    … …  
    349361                          old_oa_i(:,:,:)  = oa_i(:,:,:) 
    350362                          old_smv_i(:,:,:) = smv_i (:,:,:) 
    351          !                                           ! Ice thermodynamics  
     363  
     364         ! ---------------------------------------------- 
     365         ! ice thermodynamic 
     366         ! ---------------------------------------------- 
    352367                          CALL lim_var_glo2eqv            ! equivalent variables 
    353368                          CALL lim_var_agg(1)             ! aggregate ice categories 
     369                          ! previous lead fraction and ice volume for flux calculations 
     370                          pfrld(:,:)   = 1._wp - at_i(:,:) 
     371                          phicif(:,:)  = vt_i(:,:) 
     372                          ! 
    354373                          CALL lim_var_bv                 ! bulk brine volume (diag) 
    355374                          CALL lim_thd( kt )              ! Ice thermodynamics  
    … …  
    357376                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    358377                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin) 
    359          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
     378         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    360379                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
    361          ! 
    362          !                                           ! Global variables update 
    363380                          CALL lim_var_agg( 1 )           ! requested by limupdate 
    364                           CALL lim_update2                 ! Global variables update 
    365 #if defined key_bdy 
    366                           CALL lim_var_glo2eqv            ! 
    367                           CALL bdy_ice_lim( kt )          ! clem modif: bdy ice 
    368                           CALL bdy_ice_lim_dyn( 1 )       !?? repeated from limrhg.F90 
    369 #endif 
     381                          CALL lim_update2                ! Global variables update 
     382 
    370383                          CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
    371384                          CALL lim_var_agg(2)             ! aggregate ice thickness categories 
    372          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
     385         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    373386         ! 
    374387                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes 
    375388         ! 
    376          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
     389         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
    377390         ! 
    378391         !                                           ! Diagnostics and outputs  
    … …  
    387400                          CALL lim_var_glo2eqv            ! ??? 
    388401         ! 
    389          IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash 
     402         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash 
    390403         ! 
    391404      ENDIF                                    ! End sea-ice time step only 
    … …  
    414427 
    415428 
    416    SUBROUTINE lim_ctl 
     429   SUBROUTINE lim_ctl( kt ) 
    417430      !!----------------------------------------------------------------------- 
    418431      !!                   ***  ROUTINE lim_ctl ***  
    … …  
    420433      !! ** Purpose :   Alerts in case of model crash 
    421434      !!------------------------------------------------------------------- 
     435      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    422436      INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices 
    423437      INTEGER  ::   inb_altests       ! number of alert tests (max 20) 
    … …  
    439453            DO ji = 1, jpi 
    440454               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    441                   WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    442                   WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    443                   WRITE(numout,*) ' Point - category', ji, jj, jl 
    444                   WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl) 
    445                   WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl) 
    446                   WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    447                   WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
     455                  !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
     456                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
     457                  !WRITE(numout,*) ' Point - category', ji, jj, jl 
     458                  !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl) 
     459                  !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl) 
     460                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
     461                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
    448462                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    449463               ENDIF 
    … …  
    459473         DO ji = 1, jpi 
    460474            IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN 
    461                CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
     475               !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    462476               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    463477            ENDIF 
    … …  
    470484      DO jj = 1, jpj 
    471485         DO ji = 1, jpi 
    472             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5  .AND.  & 
     486            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
    473487               &  at_i(ji,jj) > 0._wp   ) THEN 
    474                CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    475                WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    476                WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    477                WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    478                WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    479                WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    480                WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
    481                WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
    482                WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    483                WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    484                WRITE(numout,*)  
     488               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
     489               !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
     490               !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
     491               !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
     492               !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
     493               !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
     494               !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
     495               !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
     496               !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
     497               !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
     498               !WRITE(numout,*)  
    485499               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    486500            ENDIF 
    … …  
    494508         DO ji = 1, jpi 
    495509            IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    496                CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    497                WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
    498                WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    499                WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    500                WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    501                WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    502                WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    503                WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    504                WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
     510               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
     511               !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
     512               !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
     513               !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
     514               !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
     515               !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
     516               !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
     517               !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
     518               !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    505519               ! 
    506520               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    … …  
    516530         DO jj = 1, jpj 
    517531            DO ji = 1, jpi 
    518 !!gm  test twice sm_i ...  ????  bug? 
    519                IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 )   .OR.  & 
    520                      ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. & 
    521                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    522 !                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
     532               IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     533!                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    523534!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    524535!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    … …  
    541552                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    542553                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    543                   CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
     554                  !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    544555                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    545556               ENDIF 
    … …  
    553564      DO jj = 1, jpj 
    554565         DO ji = 1, jpi 
    555             IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN 
    556                CALL lim_prt_state( ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    557                DO jl = 1, jpl 
    558                   WRITE(numout,*) ' Category no: ', jl 
    559                   WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)    
    560                   WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    561                   WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)    
    562                   WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    563                   WRITE(numout,*) ' ' 
    564                END DO 
     566            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
     567               !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
     568               !DO jl = 1, jpl 
     569                  !WRITE(numout,*) ' Category no: ', jl 
     570                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)    
     571                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
     572                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)    
     573                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
     574                  !WRITE(numout,*) ' ' 
     575               !END DO 
    565576               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    566577            ENDIF 
    … …  
    573584      DO jj = 1, jpj 
    574585         DO ji = 1, jpi 
    575             IF(   ABS( qns(ji,jj) ) > 1500._wp  .AND.  at_i(ji,jj) > 0._wp  ) THEN 
     586            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    576587               ! 
    577                WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    578                WRITE(numout,*) ' ji, jj    : ', ji, jj 
    579                WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    580                WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    581                WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    582                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
    583                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
    584                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
    585                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
    586                WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
    587                WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
    588                WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
    589                WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
    590                WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
    591                WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
    592                WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
    593                WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
    594                WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
     588               !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
     589               !WRITE(numout,*) ' ji, jj    : ', ji, jj 
     590               !WRITE(numout,*) ' qns       : ', qns(ji,jj) 
     591               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
     592               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
     593               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
     594               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
     595               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
     596               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
     597               !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
     598               !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
     599               !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
     600               !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
     601               !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
     602               !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
     603               !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
     604               !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
     605               !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
    595606               ! 
    596                CALL lim_prt_state( ji, jj, 2, '   ') 
     607               !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
    597608               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    598609               ! 
    … …  
    611622               DO ji = 1, jpi 
    612623                  ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt 
    613                   IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-6   & 
     624                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    614625                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    615                      WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    616                      WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
    617                      WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
    618                      WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
    619                      WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 
    620                      WRITE(numout,*) ' ztmelts : ', ztmelts 
     626                     !WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
     627                     !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
     628                     !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
     629                     !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
     630                     !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 
     631                     !WRITE(numout,*) ' ztmelts : ', ztmelts 
    621632                     inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    622633                  ENDIF 
    … …  
    626637      END DO 
    627638 
    628       ialert_id = 1                                 ! reference number of this alert 
    629       cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
    630       WRITE(numout,*) 
    631       WRITE(numout,*) ' All alerts at the end of ice model ' 
    632       DO ialert_id = 1, inb_altests 
    633          WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
    634       END DO 
    635       ! 
     639      ! sum of the alerts on all processors 
     640      IF( lk_mpp ) THEN 
     641         DO ialert_id = 1, inb_altests 
     642            CALL mpp_sum(inb_alp(ialert_id)) 
     643         END DO 
     644      ENDIF 
     645 
     646      ! print alerts 
     647      IF( lwp ) THEN 
     648         ialert_id = 1                                 ! reference number of this alert 
     649         cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
     650         WRITE(numout,*) ' time step ',kt 
     651         WRITE(numout,*) ' All alerts at the end of ice model ' 
     652         DO ialert_id = 1, inb_altests 
     653            WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
     654         END DO 
     655      ENDIF 
     656     ! 
    636657   END SUBROUTINE lim_ctl 
    637658  
    638659    
    639    SUBROUTINE lim_prt_state( ki, kj, kn, cd1 ) 
     660   SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 ) 
    640661      !!----------------------------------------------------------------------- 
    641662      !!                   ***  ROUTINE lim_prt_state ***  
    … …  
    651672      !!                n : number of the option 
    652673      !!------------------------------------------------------------------- 
     674      INTEGER         , INTENT(in) ::   kt      ! ocean time step 
    653675      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    654676      CHARACTER(len=*), INTENT(in) ::   cd1           ! 
    655677      !! 
    656       INTEGER :: jl 
     678      INTEGER :: jl, ji, jj 
    657679      !!------------------------------------------------------------------- 
    658680 
    659       WRITE(numout,*) cd1             ! print title 
    660  
    661       !---------------- 
    662       !  Simple state 
    663       !---------------- 
    664  
    665       IF ( kn == 1 .OR. kn == -1 ) THEN 
    666          WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 
    667          WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    668          WRITE(numout,*) ' Simple state ' 
    669          WRITE(numout,*) ' masks s,u,v   : ', tms(ki,kj), tmu(ki,kj), tmv(ki,kj) 
    670          WRITE(numout,*) ' lat - long    : ', gphit(ki,kj), glamt(ki,kj) 
    671          WRITE(numout,*) ' Time step     : ', numit 
    672          WRITE(numout,*) ' - Ice drift   ' 
    673          WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    674          WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj) 
    675          WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj) 
    676          WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1) 
    677          WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj) 
    678          WRITE(numout,*) ' strength      : ', strength(ki,kj) 
    679          WRITE(numout,*) 
    680          WRITE(numout,*) ' - Cell values ' 
    681          WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    682          WRITE(numout,*) ' cell area     : ', area(ki,kj) 
    683          WRITE(numout,*) ' at_i          : ', at_i(ki,kj)        
    684          WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)        
    685          WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)        
    686          DO jl = 1, jpl 
    687             WRITE(numout,*) ' - Category (', jl,')' 
    688             WRITE(numout,*) ' a_i           : ', a_i(ki,kj,jl) 
    689             WRITE(numout,*) ' ht_i          : ', ht_i(ki,kj,jl) 
    690             WRITE(numout,*) ' ht_s          : ', ht_s(ki,kj,jl) 
    691             WRITE(numout,*) ' v_i           : ', v_i(ki,kj,jl) 
    692             WRITE(numout,*) ' v_s           : ', v_s(ki,kj,jl) 
    693             WRITE(numout,*) ' e_s           : ', e_s(ki,kj,1,jl)/1.0e9 
    694             WRITE(numout,*) ' e_i           : ', e_i(ki,kj,1:nlay_i,jl)/1.0e9 
    695             WRITE(numout,*) ' t_su          : ', t_su(ki,kj,jl) 
    696             WRITE(numout,*) ' t_snow        : ', t_s(ki,kj,1,jl) 
    697             WRITE(numout,*) ' t_i           : ', t_i(ki,kj,1:nlay_i,jl) 
    698             WRITE(numout,*) ' sm_i          : ', sm_i(ki,kj,jl) 
    699             WRITE(numout,*) ' smv_i         : ', smv_i(ki,kj,jl) 
    700             WRITE(numout,*) 
    701             WRITE(numout,*) ' Pathological case : ', patho_case(ki,kj,jl) 
    702          END DO 
    703       ENDIF 
    704       IF( kn == -1 ) THEN 
    705          WRITE(numout,*) ' Mechanical Check ************** ' 
    706          WRITE(numout,*) ' Check what means ice divergence ' 
    707          WRITE(numout,*) ' Total ice concentration ', at_i (ki,kj) 
    708          WRITE(numout,*) ' Total lead fraction     ', ato_i(ki,kj) 
    709          WRITE(numout,*) ' Sum of both             ', ato_i(ki,kj) + at_i(ki,kj) 
    710          WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ki,kj) + at_i(ki,kj) - 1.00 
    711       ENDIF 
    712  
    713  
    714       !-------------------- 
    715       !  Exhaustive state 
    716       !-------------------- 
    717  
    718       IF ( kn .EQ. 2 ) THEN 
    719          WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 
    720          WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    721          WRITE(numout,*) ' Exhaustive state ' 
    722          WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj) 
    723          WRITE(numout,*) ' Time step ', numit 
    724          WRITE(numout,*)  
    725          WRITE(numout,*) ' - Cell values ' 
    726          WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    727          WRITE(numout,*) ' cell area     : ', area(ki,kj) 
    728          WRITE(numout,*) ' at_i          : ', at_i(ki,kj)        
    729          WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)        
    730          WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)        
    731          WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj) 
    732          WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj) 
    733          WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1) 
    734          WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj) 
    735          WRITE(numout,*) ' strength      : ', strength(ki,kj) 
    736          WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ki,kj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ki,kj) 
    737          WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ki,kj)  , ' old_v_ice     : ', old_v_ice(ki,kj)   
    738          WRITE(numout,*) 
    739  
    740          DO jl = 1, jpl 
    741               WRITE(numout,*) ' - Category (',jl,')' 
    742               WRITE(numout,*) '   ~~~~~~~~         '  
    743               WRITE(numout,*) ' ht_i       : ', ht_i(ki,kj,jl)             , ' ht_s       : ', ht_s(ki,kj,jl) 
    744               WRITE(numout,*) ' t_i        : ', t_i(ki,kj,1:nlay_i,jl) 
    745               WRITE(numout,*) ' t_su       : ', t_su(ki,kj,jl)             , ' t_s        : ', t_s(ki,kj,1,jl) 
    746               WRITE(numout,*) ' sm_i       : ', sm_i(ki,kj,jl)             , ' o_i        : ', o_i(ki,kj,jl) 
    747               WRITE(numout,*) ' a_i        : ', a_i(ki,kj,jl)              , ' old_a_i    : ', old_a_i(ki,kj,jl)    
    748               WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ki,kj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ki,kj,jl)  
    749               WRITE(numout,*) ' v_i        : ', v_i(ki,kj,jl)              , ' old_v_i    : ', old_v_i(ki,kj,jl)    
    750               WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ki,kj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ki,kj,jl)  
    751               WRITE(numout,*) ' v_s        : ', v_s(ki,kj,jl)              , ' old_v_s    : ', old_v_s(ki,kj,jl)   
    752               WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ki,kj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ki,kj,jl) 
    753               WRITE(numout,*) ' e_i1       : ', e_i(ki,kj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ki,kj,1,jl)/1.0e9  
    754               WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ki,kj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ki,kj,1,jl)/1.0e9 
    755               WRITE(numout,*) ' e_i2       : ', e_i(ki,kj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ki,kj,2,jl)/1.0e9   
    756               WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ki,kj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ki,kj,2,jl)/1.0e9 
    757               WRITE(numout,*) ' e_snow     : ', e_s(ki,kj,1,jl)            , ' old_e_snow : ', old_e_s(ki,kj,1,jl)  
    758               WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ki,kj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ki,kj,1,jl) 
    759               WRITE(numout,*) ' smv_i      : ', smv_i(ki,kj,jl)            , ' old_smv_i  : ', old_smv_i(ki,kj,jl)    
    760               WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ki,kj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ki,kj,jl)  
    761               WRITE(numout,*) ' oa_i       : ', oa_i(ki,kj,jl)             , ' old_oa_i   : ', old_oa_i(ki,kj,jl) 
    762               WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ki,kj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ki,kj,jl) 
    763               WRITE(numout,*) ' Path. case : ', patho_case(ki,kj,jl) 
    764         END DO !jl 
    765  
    766         WRITE(numout,*) 
    767         WRITE(numout,*) ' - Heat / FW fluxes ' 
    768         WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    769 !       WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ki,kj) 
    770 !       WRITE(numout,*) ' sfx        : ', sfx      (ki,kj) 
    771 !       WRITE(numout,*) ' sfx_res  : ', sfx_res(ki,kj) 
    772         WRITE(numout,*) ' fmmec      : ', fmmec    (ki,kj) 
    773         WRITE(numout,*) ' fhmec      : ', fhmec    (ki,kj) 
    774         WRITE(numout,*) ' fhbri      : ', fhbri    (ki,kj) 
    775         WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ki,kj) 
    776         WRITE(numout,*)  
    777         WRITE(numout,*) ' sst        : ', sst_m(ki,kj)   
    778         WRITE(numout,*) ' sss        : ', sss_m(ki,kj)   
    779         WRITE(numout,*)  
    780         WRITE(numout,*) ' - Stresses ' 
    781         WRITE(numout,*) '   ~~~~~~~~ ' 
    782         WRITE(numout,*) ' utau_ice   : ', utau_ice(ki,kj)  
    783         WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ki,kj) 
    784         WRITE(numout,*) ' utau       : ', utau    (ki,kj)  
    785         WRITE(numout,*) ' vtau       : ', vtau    (ki,kj) 
    786         WRITE(numout,*) ' oc. vel. u : ', u_oce   (ki,kj) 
    787         WRITE(numout,*) ' oc. vel. v : ', v_oce   (ki,kj) 
    788      ENDIF 
    789  
    790      !--------------------- 
    791      ! Salt / heat fluxes 
    792      !--------------------- 
    793  
    794      IF ( kn .EQ. 3 ) THEN 
    795         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 
    796         WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    797         WRITE(numout,*) ' - Salt / Heat Fluxes ' 
    798         WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    799         WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj) 
    800         WRITE(numout,*) ' Time step ', numit 
    801         WRITE(numout,*) 
    802         WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    803         WRITE(numout,*) ' qsr       : ', qsr(ki,kj) 
    804         WRITE(numout,*) ' qns       : ', qns(ki,kj) 
    805         WRITE(numout,*) 
    806         WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    807         WRITE(numout,*) ' emp       : ', emp    (ki,kj) 
    808         WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ki,kj) 
    809         WRITE(numout,*) ' sfx       : ', sfx    (ki,kj) 
    810         WRITE(numout,*) ' sfx_res   : ', sfx_res(ki,kj) 
    811         WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ki,kj) 
    812         WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    813         WRITE(numout,*) ' fheat_res : ', fheat_res(ki,kj) 
    814         WRITE(numout,*) 
    815         WRITE(numout,*) ' - Momentum fluxes ' 
    816         WRITE(numout,*) ' utau      : ', utau(ki,kj)  
    817         WRITE(numout,*) ' vtau      : ', vtau(ki,kj) 
    818       ENDIF 
    819       WRITE(numout,*) ' ' 
    820       ! 
     681      DO ji = mi0(ki), mi1(ki) 
     682         DO jj = mj0(kj), mj1(kj) 
     683 
     684            WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title 
     685 
     686            !---------------- 
     687            !  Simple state 
     688            !---------------- 
     689             
     690            IF ( kn == 1 .OR. kn == -1 ) THEN 
     691               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
     692               WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     693               WRITE(numout,*) ' Simple state ' 
     694               WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
     695               WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj) 
     696               WRITE(numout,*) ' Time step     : ', numit 
     697               WRITE(numout,*) ' - Ice drift   ' 
     698               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     699               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
     700               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
     701               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
     702               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
     703               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
     704               WRITE(numout,*) 
     705               WRITE(numout,*) ' - Cell values ' 
     706               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     707               WRITE(numout,*) ' cell area     : ', area(ji,jj) 
     708               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     709               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     710               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
     711               DO jl = 1, jpl 
     712                  WRITE(numout,*) ' - Category (', jl,')' 
     713                  WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl) 
     714                  WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl) 
     715                  WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl) 
     716                  WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl) 
     717                  WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl) 
     718                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9 
     719                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 
     720                  WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl) 
     721                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl) 
     722                  WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl) 
     723                  WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl) 
     724                  WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl) 
     725                  WRITE(numout,*) 
     726               END DO 
     727            ENDIF 
     728            IF( kn == -1 ) THEN 
     729               WRITE(numout,*) ' Mechanical Check ************** ' 
     730               WRITE(numout,*) ' Check what means ice divergence ' 
     731               WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
     732               WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
     733               WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
     734               WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
     735            ENDIF 
     736             
     737 
     738            !-------------------- 
     739            !  Exhaustive state 
     740            !-------------------- 
     741             
     742            IF ( kn .EQ. 2 ) THEN 
     743               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
     744               WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     745               WRITE(numout,*) ' Exhaustive state ' 
     746               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
     747               WRITE(numout,*) ' Time step ', numit 
     748               WRITE(numout,*)  
     749               WRITE(numout,*) ' - Cell values ' 
     750               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     751               WRITE(numout,*) ' cell area     : ', area(ji,jj) 
     752               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     753               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     754               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
     755               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
     756               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
     757               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
     758               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
     759               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
     760               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
     761               WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj)   
     762               WRITE(numout,*) 
     763                
     764               DO jl = 1, jpl 
     765                  WRITE(numout,*) ' - Category (',jl,')' 
     766                  WRITE(numout,*) '   ~~~~~~~~         '  
     767                  WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl) 
     768                  WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl) 
     769                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
     770                  WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
     771                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)    
     772                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
     773                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)    
     774                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
     775                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl)   
     776                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
     777                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9  
     778                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9 
     779                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9   
     780                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9 
     781                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl)  
     782                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl) 
     783                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)    
     784                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
     785                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl) 
     786                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
     787               END DO !jl 
     788                
     789               WRITE(numout,*) 
     790               WRITE(numout,*) ' - Heat / FW fluxes ' 
     791               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
     792               WRITE(numout,*) ' emp        : ', emp      (ji,jj) 
     793               WRITE(numout,*) ' sfx        : ', sfx      (ji,jj) 
     794               WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj) 
     795               WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj) 
     796               WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj) 
     797               WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj) 
     798               WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj) 
     799               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj) 
     800               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj) 
     801               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj) 
     802               WRITE(numout,*)  
     803               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
     804               WRITE(numout,*) ' sss        : ', sss_m(ji,jj)   
     805               WRITE(numout,*)  
     806               WRITE(numout,*) ' - Stresses ' 
     807               WRITE(numout,*) '   ~~~~~~~~ ' 
     808               WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj)  
     809               WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj) 
     810               WRITE(numout,*) ' utau       : ', utau    (ji,jj)  
     811               WRITE(numout,*) ' vtau       : ', vtau    (ji,jj) 
     812               WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj) 
     813               WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj) 
     814            ENDIF 
     815             
     816            !--------------------- 
     817            ! Salt / heat fluxes 
     818            !--------------------- 
     819             
     820            IF ( kn .EQ. 3 ) THEN 
     821               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
     822               WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     823               WRITE(numout,*) ' - Salt / Heat Fluxes ' 
     824               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
     825               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
     826               WRITE(numout,*) ' Time step ', numit 
     827               WRITE(numout,*) 
     828               WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
     829               WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
     830               WRITE(numout,*) ' qns       : ', qns(ji,jj) 
     831               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
     832               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice 
     833               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice 
     834               WRITE(numout,*) 
     835               WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
     836               WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
     837               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
     838               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
     839               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
     840               WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj) 
     841               WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
     842               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
     843               WRITE(numout,*) 
     844               WRITE(numout,*) ' - Momentum fluxes ' 
     845               WRITE(numout,*) ' utau      : ', utau(ji,jj)  
     846               WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
     847            ENDIF 
     848            WRITE(numout,*) ' ' 
     849            ! 
     850         END DO 
     851      END DO 
     852 
    821853   END SUBROUTINE lim_prt_state 
    822854 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4230 r4333  
    343343      !                                                           ! (update freshwater fluxes) 
    344344!RBbug do not understand why see ticket 667 
    345       CALL lbc_lnk( emp, 'T', 1. ) 
     345      !clem-bugsal CALL lbc_lnk( emp, 'T', 1. ) 
    346346      ! 
    347347      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 

  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4292 r4333  
    3333   USE timing         ! Timing 
    3434   USE sbc_ice, ONLY : lk_lim3 
    35  
    3635 
    3736   IMPLICIT NONE 
Note: See TracChangeset for help on using the changeset viewer.