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 5051 for branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2015-02-02T18:31:34+01:00 (9 years ago)
Author:
clem
Message:

LIM3 initialization is now called at the same time as other sbc fields

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5048 r5051  
    1919   !!---------------------------------------------------------------------- 
    2020   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area 
    21    !!   lim_ctl       : alerts in case of ice model crash 
    22    !!   lim_prt_state : ice control print at a given grid point 
    2321   !!---------------------------------------------------------------------- 
    2422   USE oce             ! ocean dynamics and tracers 
     
    2624   USE par_ice         ! sea-ice parameters 
    2725   USE ice             ! LIM-3: ice variables 
    28    USE iceini          ! LIM-3: ice initialisation 
     26   USE thd_ice         ! LIM-3: thermodynamical variables 
    2927   USE dom_ice         ! LIM-3: ice domain 
    3028 
     
    4745   USE limwri          ! Ice outputs 
    4846   USE limrst          ! Ice restarts 
    49    USE limupdate1       ! update of global variables 
    50    USE limupdate2       ! update of global variables 
     47   USE limupdate1      ! update of global variables 
     48   USE limupdate2      ! update of global variables 
    5149   USE limvar          ! Ice variables switch 
     50 
     51   USE limmsh          ! LIM mesh 
     52   USE limistate       ! LIM initial state 
     53   USE limthd_sal      ! LIM ice thermodynamics: salinity 
    5254 
    5355   USE c1d             ! 1D vertical configuration 
     
    6062   USE prtctl          ! Print control 
    6163   USE lib_fortran     !  
     64   USE limctl 
    6265 
    6366#if defined key_bdy  
     
    6972 
    7073   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
    71    PUBLIC lim_prt_state 
     74   PUBLIC sbc_lim_init ! routine called by sbcmod.F90 
    7275    
    7376   !! * Substitutions 
     
    114117      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    115118 
    116       IF( kt == nit000 ) THEN 
    117          IF(lwp) WRITE(numout,*) 
    118          IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
    119          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
    120          ! 
    121          CALL ice_init 
    122          ! 
    123          IF( ln_nicep ) THEN      ! control print at a given point 
    124             jiindx = 15    ;   jjindx =  44 
    125             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    126          ENDIF 
    127       ENDIF 
    128  
    129119      !                                        !----------------------! 
    130120      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
     
    133123         !                                           !----------------! 
    134124         ! 
    135          u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point 
    136          v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)                    ! (C-grid dynamics :  U- & V-points as the ocean) 
    137          ! 
    138          t_bo(:,:) = ( eos_fzp( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )  ! masked sea surface freezing temperature [Kelvin] 
    139          !                                                                                  ! (set to rt0 over land) 
    140          !                                           ! Ice albedo 
    141          CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )       
    142  
     125         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)      ! mean surface ocean current at ice velocity point 
     126         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)      ! (C-grid dynamics :  U- & V-points as the ocean) 
     127          
     128         ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
     129         t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )   
     130         !                                                                                       
     131         ! Ice albedo 
     132         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
    143133         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
    144134 
     
    153143         END SELECT 
    154144          
    155          !                                           ! Mask sea ice surface temperature 
     145         ! Mask sea ice surface temperature (set to rt0 over land) 
    156146         DO jl = 1, jpl 
    157             t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) ) 
     147            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
    158148         END DO 
    159149      
     
    210200         v_ice_b(:,:)     = v_ice(:,:) 
    211201 
    212          ! salt, heat and mass fluxes 
    213          sfx    (:,:) = 0._wp   ; 
    214          sfx_bri(:,:) = 0._wp   ;  
    215          sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    216          sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    217          sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    218          sfx_res(:,:) = 0._wp 
    219  
    220          wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    221          wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
    222          wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
    223          wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    224          wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    225          wfx_spr(:,:) = 0._wp   ;    
    226  
    227          afx_tot(:,:) = at_i(:,:) ;  afx_dyn(:,:) = 0._wp 
    228          afx_thd(:,:) = 0._wp 
    229  
    230          hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
    231          hfx_thd(:,:) = 0._wp   ;    
    232          hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
    233          hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
    234          hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
    235          hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
    236          hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
    237          hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     202                          CALL sbc_lim_flx0          ! set diag of mass, heat and salt fluxes to 0 
    238203 
    239204                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
    240205         ! 
    241          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
     206         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    242207         ! ---------------------------------------------- 
    243208         ! ice dynamics and transport (except in 1D case) 
     
    245210         IF( .NOT. lk_c1d ) THEN 
    246211 
    247          IF ( ln_limdyn ) afx_dyn(:,:)     = at_i(:,:) 
    248  
    249212                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    250213                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    251214                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    252          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
    253          IF( nn_monocat /= 2 ) CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
     215         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
     216         IF( nn_monocat /= 2 )   & 
     217            &             CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    254218                          CALL lim_var_agg( 1 )  
    255219#if defined key_bdy 
     
    257221                          CALL lim_var_glo2eqv            ! equivalent variables 
    258222                          CALL bdy_ice_lim( kt ) 
    259                           CALL lim_itd_me_zapsmall 
     223                          CALL lim_var_zapsmall 
    260224                          CALL lim_var_agg(1) 
    261          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
     225         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
    262226#endif 
    263227                          CALL lim_update1 
     
    273237                          oa_i_b (:,:,:)   = oa_i (:,:,:) 
    274238                          smv_i_b(:,:,:)   = smv_i(:,:,:) 
    275  
    276          IF ( ln_limdyn ) afx_dyn(:,:)     = ( at_i(:,:) - afx_dyn(:,:) ) * r1_rdtice 
    277                           afx_thd(:,:)     = at_i(:,:) 
    278239  
    279240         ! ---------------------------------------------- 
     
    282243                          CALL lim_var_glo2eqv            ! equivalent variables 
    283244                          CALL lim_var_agg(1)             ! aggregate ice categories 
     245 
    284246                          ! previous lead fraction and ice volume for flux calculations 
    285247                          pfrld(:,:)   = 1._wp - at_i(:,:) 
     
    301263                          zcoef = rdt_ice /rday           !  Ice natural aging 
    302264                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    303          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    304                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
    305                           CALL lim_var_agg( 1 )           ! requested by limupdate 
     265         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
     266                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion 
    306267                          CALL lim_update2                ! Global variables update 
    307268 
    308                           CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
    309                           CALL lim_var_agg(2)             ! aggregate ice thickness categories 
    310          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
     269         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    311270         ! 
    312271                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes 
    313272         ! 
    314          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
    315          ! 
    316          !                                           ! Diagnostics and outputs  
    317          IF (ln_limdiaout) CALL lim_diahsb 
    318  
    319                           afx_thd(:,:)     = ( at_i(:,:) - afx_thd(:,:) ) * r1_rdtice 
    320                           afx_tot(:,:)     = ( at_i(:,:) - afx_tot(:,:) ) * r1_rdtice 
     273         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
     274         ! 
     275         IF(ln_limdiaout) CALL lim_diahsb                 ! Diagnostics and outputs  
    321276 
    322277                          CALL lim_wri( 1  )              ! Ice outputs  
     
    342297      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
    343298!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    344  
    345299      ! 
    346300      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
     
    348302   END SUBROUTINE sbc_ice_lim 
    349303    
     304 
     305   SUBROUTINE sbc_lim_init 
     306      !!---------------------------------------------------------------------- 
     307      !!                  ***  ROUTINE sbc_lim_init  *** 
     308      !! 
     309      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules 
     310      !!---------------------------------------------------------------------- 
     311      INTEGER :: ierr 
     312      !!---------------------------------------------------------------------- 
     313      IF(lwp) WRITE(numout,*) 
     314      IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
     315      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
     316      ! 
     317      ! 
     318      !                                ! Allocate the ice arrays 
     319      ierr =        ice_alloc        ()      ! ice variables 
     320      ierr = ierr + dom_ice_alloc    ()      ! domain 
     321      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
     322      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics 
     323      ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
     324      ! 
     325      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     326      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays') 
     327      ! 
     328      !                                ! adequation jpk versus ice/snow layers/categories 
     329      IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk )   & 
     330         &      CALL ctl_stop( 'STOP',                     & 
     331         &     'sbc_lim_init: the 3rd dimension of workspace arrays is too small.',   & 
     332         &     'use more ocean levels or less ice/snow layers/categories.' ) 
     333 
     334                                       ! Open the reference and configuration namelist files and namelist output file  
     335      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )  
     336      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     337      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) 
     338      ! 
     339      CALL ice_run                     ! set some ice run parameters 
     340      ! 
     341      CALL lim_thd_init                ! set ice thermodynics parameters 
     342      ! 
     343      CALL lim_thd_sal_init            ! set ice salinity parameters 
     344      ! 
     345      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep 
     346      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse 
     347      ! 
     348      CALL lim_msh                     ! ice mesh initialization 
     349      ! 
     350      CALL lim_itd_init                 ! ice thickness distribution initialization 
     351      ! 
     352      CALL lim_itd_me_init             ! ice thickness distribution initialization 
     353      !                                ! Initial sea-ice state 
     354      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
     355         numit = 0 
     356         numit = nit000 - 1 
     357         CALL lim_istate 
     358         CALL lim_var_agg(1) 
     359         CALL lim_var_glo2eqv 
     360      ELSE                                    ! start from a restart file 
     361         CALL lim_rst_read 
     362         numit = nit000 - 1 
     363         CALL lim_var_agg(1) 
     364         CALL lim_var_glo2eqv 
     365      ENDIF 
     366      ! 
     367      CALL lim_sbc_init                 ! ice surface boundary condition    
     368      ! 
     369      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
     370      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
     371      ! 
     372      nstart = numit  + nn_fsbc       
     373      nitrun = nitend - nit000 + 1  
     374      nlast  = numit  + nitrun  
     375      ! 
     376      IF( nstock == 0 )   nstock = nlast + 1 
     377      ! 
     378   END SUBROUTINE sbc_lim_init 
     379 
     380 
     381   SUBROUTINE ice_run 
     382      !!------------------------------------------------------------------- 
     383      !!                  ***  ROUTINE ice_run *** 
     384      !!                  
     385      !! ** Purpose :   Definition some run parameter for ice model 
     386      !! 
     387      !! ** Method  :   Read the namicerun namelist and check the parameter  
     388      !!              values called at the first timestep (nit000) 
     389      !! 
     390      !! ** input   :   Namelist namicerun 
     391      !!------------------------------------------------------------------- 
     392      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout 
     393      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     394      !!------------------------------------------------------------------- 
     395      !                     
     396      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice 
     397      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 
     398901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 
     399 
     400      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice 
     401      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 
     402902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 
     403      IF(lwm) WRITE ( numoni, namicerun ) 
     404      ! 
     405      ! 
     406      IF(lwp) THEN                        ! control print 
     407         WRITE(numout,*) 
     408         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
     409         WRITE(numout,*) ' ~~~~~~' 
     410         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
     411         WRITE(numout,*) '   maximum ice concentration                               = ', amax  
     412         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep 
     413         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
     414         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
     415      ENDIF 
     416      ! 
     417      !IF( lk_mpp .AND. ln_nicep ) THEN 
     418      !   ln_nicep = .FALSE. 
     419      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
     420      !ENDIF 
     421      IF( ln_nicep ) THEN      ! control print at a given point 
     422         jiindx = 15    ;   jjindx =  44 
     423         IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
     424      ENDIF 
     425      ! 
     426   END SUBROUTINE ice_run 
     427 
     428 
     429   SUBROUTINE lim_itd_init 
     430      !!------------------------------------------------------------------ 
     431      !!                ***  ROUTINE lim_itd_init *** 
     432      !! 
     433      !! ** Purpose :   Initializes the ice thickness distribution 
     434      !! ** Method  :   ... 
     435      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for 
     436      !!              limistate (only) and is changed to 99 m in sbc_lim_init 
     437      !!------------------------------------------------------------------ 
     438      INTEGER  ::   jl                   ! dummy loop index 
     439      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars 
     440      REAL(wp) ::   zhmax, znum, zden, zalpha ! 
     441      !!------------------------------------------------------------------ 
     442 
     443      IF(lwp) WRITE(numout,*) 
     444      IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice thickness distribution ' 
     445      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     446 
     447      !------------------------------------------------------------------------------! 
     448      ! 1) Ice thickness distribution parameters initialization     
     449      !------------------------------------------------------------------------------! 
     450      IF(lwp) WRITE(numout,*) ' Number of ice categories jpl = ', jpl 
     451 
     452      !- Thickness categories boundaries  
     453      !---------------------------------- 
     454      !  Clem - je sais pas encore dans quelle namelist les mettre, ca depend des chgts liés à bdy 
     455      nn_hibnd  = 2                !  thickness category boundaries: tanh function (1) h^(-alpha) (2) 
     456      rn_hibnd  = 2.5              !  mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only) 
     457 
     458      hi_max(:) = 0._wp 
     459 
     460      SELECT CASE ( nn_hibnd  )        
     461                                   !---------------------- 
     462         CASE (1)                  ! tanh function (CICE) 
     463                                   !---------------------- 
     464         zc1 =  3._wp / REAL( jpl, wp ) 
     465         zc2 = 10._wp * zc1 
     466         zc3 =  3._wp 
     467 
     468         DO jl = 1, jpl 
     469            zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 
     470            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 
     471         END DO 
     472 
     473                                   !---------------------- 
     474         CASE (2)                  ! h^(-alpha) function 
     475                                   !---------------------- 
     476         zalpha = 0.05             ! exponent of the transform function 
     477 
     478         zhmax  = 3.*rn_hibnd 
     479 
     480         DO jl = 1, jpl  
     481            znum = jpl * ( zhmax+1 )**zalpha 
     482            zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl 
     483            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
     484         END DO 
     485 
     486      END SELECT 
     487 
     488      DO jl = 1, jpl 
     489         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     490      END DO 
     491      ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 
     492      hi_max(jpl) = 99._wp 
     493 
     494      IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' 
     495      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl) 
     496      ! 
     497   END SUBROUTINE lim_itd_init 
     498 
    350499    
    351       SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
     500   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
    352501         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 
    353502      !!--------------------------------------------------------------------- 
     
    427576      ! 
    428577   END SUBROUTINE ice_lim_flx 
    429     
    430     
    431    SUBROUTINE lim_ctl( kt ) 
    432       !!----------------------------------------------------------------------- 
    433       !!                   ***  ROUTINE lim_ctl ***  
    434       !!                  
    435       !! ** Purpose :   Alerts in case of model crash 
    436       !!------------------------------------------------------------------- 
    437       INTEGER, INTENT(in) ::   kt      ! ocean time step 
    438       INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices 
    439       INTEGER  ::   inb_altests       ! number of alert tests (max 20) 
    440       INTEGER  ::   ialert_id         ! number of the current alert 
    441       REAL(wp) ::   ztmelts           ! ice layer melting point 
    442       CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert 
    443       INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive 
    444       !!------------------------------------------------------------------- 
    445  
    446       inb_altests = 10 
    447       inb_alp(:)  =  0 
    448  
    449       ! Alert if incompatible volume and concentration 
    450       ialert_id = 2 ! reference number of this alert 
    451       cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
    452  
    453       DO jl = 1, jpl 
    454          DO jj = 1, jpj 
    455             DO ji = 1, jpi 
    456                IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    457                   !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    458                   !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    459                   !WRITE(numout,*) ' Point - category', ji, jj, jl 
    460                   !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl) 
    461                   !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl) 
    462                   !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    463                   !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
    464                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    465                ENDIF 
    466             END DO 
    467          END DO 
    468       END DO 
    469  
    470       ! Alerte if very thick ice 
    471       ialert_id = 3 ! reference number of this alert 
    472       cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert 
    473       jl = jpl  
    474       DO jj = 1, jpj 
    475          DO ji = 1, jpi 
    476             IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN 
    477                !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    478                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    479             ENDIF 
    480          END DO 
    481       END DO 
    482  
    483       ! Alert if very fast ice 
    484       ialert_id = 4 ! reference number of this alert 
    485       cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert 
    486       DO jj = 1, jpj 
    487          DO ji = 1, jpi 
    488             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
    489                &  at_i(ji,jj) > 0._wp   ) THEN 
    490                !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    491                !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    492                !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    493                !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    494                !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    495                !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    496                !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
    497                !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
    498                !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    499                !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    500                !WRITE(numout,*)  
    501                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    502             ENDIF 
    503          END DO 
    504       END DO 
    505  
    506       ! Alert if there is ice on continents 
    507       ialert_id = 6 ! reference number of this alert 
    508       cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert 
    509       DO jj = 1, jpj 
    510          DO ji = 1, jpi 
    511             IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    512                !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    513                !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
    514                !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    515                !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    516                !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    517                !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    518                !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    519                !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    520                !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    521                ! 
    522                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    523             ENDIF 
    524          END DO 
    525       END DO 
    526  
    527 ! 
    528 !     ! Alert if very fresh ice 
    529       ialert_id = 7 ! reference number of this alert 
    530       cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert 
    531       DO jl = 1, jpl 
    532          DO jj = 1, jpj 
    533             DO ji = 1, jpi 
    534                IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    535 !                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    536 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    537 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    538 !                 WRITE(numout,*)  
    539                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    540                ENDIF 
    541             END DO 
    542          END DO 
    543       END DO 
    544 ! 
    545  
    546 !     ! Alert if too old ice 
    547       ialert_id = 9 ! reference number of this alert 
    548       cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert 
    549       DO jl = 1, jpl 
    550          DO jj = 1, jpj 
    551             DO ji = 1, jpi 
    552                IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. & 
    553                       ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    554                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    555                   !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    556                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    557                ENDIF 
    558             END DO 
    559          END DO 
    560       END DO 
    561   
    562       ! Alert on salt flux 
    563       ialert_id = 5 ! reference number of this alert 
    564       cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
    565       DO jj = 1, jpj 
    566          DO ji = 1, jpi 
    567             IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
    568                !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    569                !DO jl = 1, jpl 
    570                   !WRITE(numout,*) ' Category no: ', jl 
    571                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)    
    572                   !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    573                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)    
    574                   !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    575                   !WRITE(numout,*) ' ' 
    576                !END DO 
    577                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    578             ENDIF 
    579          END DO 
    580       END DO 
    581  
    582       ! Alert if qns very big 
    583       ialert_id = 8 ! reference number of this alert 
    584       cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
    585       DO jj = 1, jpj 
    586          DO ji = 1, jpi 
    587             IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    588                ! 
    589                !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    590                !WRITE(numout,*) ' ji, jj    : ', ji, jj 
    591                !WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    592                !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    593                !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    594                ! 
    595                !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
    596                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    597                ! 
    598             ENDIF 
    599          END DO 
    600       END DO 
    601       !+++++ 
    602   
    603       ! Alert if very warm ice 
    604       ialert_id = 10 ! reference number of this alert 
    605       cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert 
    606       inb_alp(ialert_id) = 0 
    607       DO jl = 1, jpl 
    608          DO jk = 1, nlay_i 
    609             DO jj = 1, jpj 
    610                DO ji = 1, jpi 
    611                   ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt 
    612                   IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    613                      &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    614                      !WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    615                      !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
    616                      !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
    617                      !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
    618                      !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 
    619                      !WRITE(numout,*) ' ztmelts : ', ztmelts 
    620                      inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    621                   ENDIF 
    622                END DO 
    623             END DO 
    624          END DO 
    625       END DO 
    626  
    627       ! sum of the alerts on all processors 
    628       IF( lk_mpp ) THEN 
    629          DO ialert_id = 1, inb_altests 
    630             CALL mpp_sum(inb_alp(ialert_id)) 
    631          END DO 
    632       ENDIF 
    633  
    634       ! print alerts 
    635       IF( lwp ) THEN 
    636          ialert_id = 1                                 ! reference number of this alert 
    637          cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
    638          WRITE(numout,*) ' time step ',kt 
    639          WRITE(numout,*) ' All alerts at the end of ice model ' 
    640          DO ialert_id = 1, inb_altests 
    641             WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
    642          END DO 
    643       ENDIF 
    644      ! 
    645    END SUBROUTINE lim_ctl 
    646   
    647     
    648    SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 ) 
    649       !!----------------------------------------------------------------------- 
    650       !!                   ***  ROUTINE lim_prt_state ***  
    651       !!                  
    652       !! ** Purpose :   Writes global ice state on the (i,j) point  
    653       !!                in ocean.ouput  
    654       !!                3 possibilities exist  
    655       !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1) 
    656       !!                n = 2    -> exhaustive state 
    657       !!                n = 3    -> ice/ocean salt fluxes 
    658       !! 
    659       !! ** input   :   point coordinates (i,j)  
    660       !!                n : number of the option 
    661       !!------------------------------------------------------------------- 
    662       INTEGER         , INTENT(in) ::   kt            ! ocean time step 
    663       INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    664       CHARACTER(len=*), INTENT(in) ::   cd1           ! 
    665       !! 
    666       INTEGER :: jl, ji, jj 
    667       !!------------------------------------------------------------------- 
    668  
    669       DO ji = mi0(ki), mi1(ki) 
    670          DO jj = mj0(kj), mj1(kj) 
    671  
    672             WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title 
    673  
    674             !---------------- 
    675             !  Simple state 
    676             !---------------- 
    677              
    678             IF ( kn == 1 .OR. kn == -1 ) THEN 
    679                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    680                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    681                WRITE(numout,*) ' Simple state ' 
    682                WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
    683                WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj) 
    684                WRITE(numout,*) ' Time step     : ', numit 
    685                WRITE(numout,*) ' - Ice drift   ' 
    686                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    687                WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
    688                WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
    689                WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
    690                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    691                WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    692                WRITE(numout,*) 
    693                WRITE(numout,*) ' - Cell values ' 
    694                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    695                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
    696                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    697                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    698                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    699                DO jl = 1, jpl 
    700                   WRITE(numout,*) ' - Category (', jl,')' 
    701                   WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl) 
    702                   WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl) 
    703                   WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl) 
    704                   WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl) 
    705                   WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl) 
    706                   WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9 
    707                   WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 
    708                   WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl) 
    709                   WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl) 
    710                   WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl) 
    711                   WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl) 
    712                   WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl) 
    713                   WRITE(numout,*) 
    714                END DO 
    715             ENDIF 
    716             IF( kn == -1 ) THEN 
    717                WRITE(numout,*) ' Mechanical Check ************** ' 
    718                WRITE(numout,*) ' Check what means ice divergence ' 
    719                WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
    720                WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
    721                WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
    722                WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
    723             ENDIF 
    724              
    725  
    726             !-------------------- 
    727             !  Exhaustive state 
    728             !-------------------- 
    729              
    730             IF ( kn .EQ. 2 ) THEN 
    731                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    732                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    733                WRITE(numout,*) ' Exhaustive state ' 
    734                WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
    735                WRITE(numout,*) ' Time step ', numit 
    736                WRITE(numout,*)  
    737                WRITE(numout,*) ' - Cell values ' 
    738                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    739                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
    740                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    741                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    742                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    743                WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
    744                WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
    745                WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
    746                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    747                WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    748                WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
    749                WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)   
    750                WRITE(numout,*) 
    751                 
    752                DO jl = 1, jpl 
    753                   WRITE(numout,*) ' - Category (',jl,')' 
    754                   WRITE(numout,*) '   ~~~~~~~~         '  
    755                   WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl) 
    756                   WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl) 
    757                   WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
    758                   WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
    759                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)    
    760                   WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    761                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)    
    762                   WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    763                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)   
    764                   WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
    765                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9  
    766                   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 
    767                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9   
    768                   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 
    769                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    770                   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) 
    771                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)    
    772                   WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
    773                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl) 
    774                   WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
    775                END DO !jl 
    776                 
    777                WRITE(numout,*) 
    778                WRITE(numout,*) ' - Heat / FW fluxes ' 
    779                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    780                WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 
    781                WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) ) 
    782                WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) ) 
    783                WRITE(numout,*) 
    784                WRITE(numout,*)  
    785                WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
    786                WRITE(numout,*) ' sss        : ', sss_m(ji,jj)   
    787                WRITE(numout,*)  
    788                WRITE(numout,*) ' - Stresses ' 
    789                WRITE(numout,*) '   ~~~~~~~~ ' 
    790                WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj)  
    791                WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj) 
    792                WRITE(numout,*) ' utau       : ', utau    (ji,jj)  
    793                WRITE(numout,*) ' vtau       : ', vtau    (ji,jj) 
    794                WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj) 
    795                WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj) 
    796             ENDIF 
    797              
    798             !--------------------- 
    799             ! Salt / heat fluxes 
    800             !--------------------- 
    801              
    802             IF ( kn .EQ. 3 ) THEN 
    803                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    804                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    805                WRITE(numout,*) ' - Salt / Heat Fluxes ' 
    806                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    807                WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
    808                WRITE(numout,*) ' Time step ', numit 
    809                WRITE(numout,*) 
    810                WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    811                WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
    812                WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    813                WRITE(numout,*) 
    814                WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj) 
    815                WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
    816                WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
    817                WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
    818                WRITE(numout,*) 
    819                WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 
    820                WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj) 
    821                WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj) 
    822                WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj)  
    823                WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice 
    824                WRITE(numout,*) 
    825                WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    826                WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
    827                WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
    828                WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
    829                WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
    830                WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj) 
    831                WRITE(numout,*) 
    832                WRITE(numout,*) ' - Momentum fluxes ' 
    833                WRITE(numout,*) ' utau      : ', utau(ji,jj)  
    834                WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
    835             ENDIF  
    836             WRITE(numout,*) ' ' 
    837             ! 
    838          END DO 
    839       END DO 
    840       ! 
    841    END SUBROUTINE lim_prt_state 
    842     
    843       
     578 
     579   SUBROUTINE sbc_lim_flx0 
     580      !!---------------------------------------------------------------------- 
     581      !!                  ***  ROUTINE sbc_lim_flx0  *** 
     582      !! 
     583      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining 
     584      !!               of the time step 
     585      !!---------------------------------------------------------------------- 
     586      sfx    (:,:) = 0._wp   ; 
     587      sfx_bri(:,:) = 0._wp   ;  
     588      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     589      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     590      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     591      sfx_res(:,:) = 0._wp 
     592       
     593      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     594      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     595      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     596      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     597      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     598      wfx_spr(:,:) = 0._wp   ;    
     599       
     600      hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
     601      hfx_thd(:,:) = 0._wp   ;    
     602      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     603      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     604      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     605      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     606      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     607      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     608 
     609      afx_tot(:,:) = 0._wp   ; 
     610      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp 
     611       
     612   END SUBROUTINE sbc_lim_flx0 
     613       
    844614   FUNCTION fice_cell_ave ( ptab ) 
    845615      !!-------------------------------------------------------------------------- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4990 r5051  
    271271      IF( ln_ssr           )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation 
    272272      ! 
     273      IF( nn_ice == 3      )   CALL sbc_lim_init               ! LIM3 initialisation 
     274 
    273275      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
    274276      ! 
Note: See TracChangeset for help on using the changeset viewer.