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 6515 for branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2016-05-09T16:42:28+02:00 (8 years ago)
Author:
clem
Message:

implement several developments for LIM3: new advection scheme (ultimate-macho, not yet perfect) ; lateral ice melt ; enabling/disabling thermo and dyn with namelist options ; simplifications (including a clarified namelist)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6399 r6515  
    140140         !----------------------------------------------------------------- 
    141141         SELECT CASE( kblk ) 
    142          CASE( jp_clio    )   ;   CALL blk_ice_clio_tau                         ! CLIO bulk formulation             
    143          CASE( jp_core    )   ;   CALL blk_ice_core_tau                         ! CORE bulk formulation 
    144          CASE( jp_purecpl )   ;   CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation 
     142            CASE( jp_clio    )   ;    CALL blk_ice_clio_tau                         ! CLIO bulk formulation             
     143            CASE( jp_core    )   ;    CALL blk_ice_core_tau                         ! CORE bulk formulation 
     144            CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation 
    145145         END SELECT 
    146146          
    147          IF( ln_mixcpl) THEN   ! Case of a mixed Bulk/Coupled formulation 
     147         IF( ln_mixcpl) THEN                                                       ! Case of a mixed Bulk/Coupled formulation 
    148148            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice) 
    149             CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 
     149                                      CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 
    150150            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
    151151            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
     
    158158         numit = numit + nn_fsbc                  ! Ice model time step 
    159159         !                                                    
    160          CALL sbc_lim_bef                         ! Store previous ice values 
    161          CALL sbc_lim_diag0                       ! set diag of mass, heat and salt fluxes to 0 
    162          CALL lim_rst_opn( kt )                   ! Open Ice restart file 
    163          ! 
    164          IF( .NOT. lk_c1d ) THEN 
     160                                      CALL sbc_lim_bef         ! Store previous ice values 
     161                                      CALL sbc_lim_diag0       ! set diag of mass, heat and salt fluxes to 0 
     162                                      CALL lim_rst_opn( kt )   ! Open Ice restart file 
     163         ! 
     164         ! --- zap this if no ice dynamics --- ! 
     165         IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN 
    165166            ! 
    166             CALL lim_dyn( kt )                    ! Ice dynamics    ( rheology/dynamics )    
    167             ! 
    168             CALL lim_trp( kt )                    ! Ice transport   ( Advection/diffusion ) 
    169             ! 
    170             IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 
    171             ! 
    172 #if defined key_bdy 
    173             CALL bdy_ice_lim( kt )                ! bdy ice thermo  
    174             IF( ln_icectl )       CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 
    175 #endif 
    176             ! 
    177             CALL lim_update1( kt )                ! Corrections 
     167            IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics 
     168                                      CALL lim_dyn( kt )       !     rheology   
     169            ELSE 
     170               u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity 
     171               v_ice(:,:) = rn_vice * vmask(:,:,1) 
     172            ENDIF 
     173                                      CALL lim_trp( kt )       ! -- Ice transport (Advection/diffusion) 
     174            IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- Mechanical redistribution (ridging/rafting) 
     175               &                      CALL lim_itd_me          
     176            IF( nn_limdyn == 2 )      CALL lim_update1( kt )   ! -- Corrections 
    178177            ! 
    179178         ENDIF 
    180           
     179         ! --- 
     180#if defined key_bdy 
     181         IF( ln_limthd )              CALL bdy_ice_lim( kt )   ! bdy ice thermo  
     182         IF( ln_icectl )              CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 
     183#endif 
     184 
    181185         ! previous lead fraction and ice volume for flux calculations 
    182          CALL sbc_lim_bef                         
    183          CALL lim_var_glo2eqv                     ! ht_i and ht_s for ice albedo calculation 
    184          CALL lim_var_agg(1)                      ! at_i for coupling (via pfrld)  
     186                                      CALL sbc_lim_bef                         
     187                                      CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation 
     188                                      CALL lim_var_agg(1)      ! at_i for coupling (via pfrld)  
     189         ! 
    185190         pfrld(:,:)   = 1._wp - at_i(:,:) 
    186191         phicif(:,:)  = vt_i(:,:) 
     
    197202         !---------------------------------------------------------------------------------------- 
    198203         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    199          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    200  
     204          
     205                                      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    201206         SELECT CASE( kblk ) 
    202          CASE( jp_clio )                                       ! CLIO bulk formulation 
    203             ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
    204             ! (alb_ice) is computed within the bulk routine 
    205                                  CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
    206             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    207             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    208          CASE( jp_core )                                       ! CORE bulk formulation 
    209             ! albedo depends on cloud fraction because of non-linear spectral effects 
    210             alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    211                                  CALL blk_ice_core_flx( t_su, alb_ice ) 
    212             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    213             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    214          CASE ( jp_purecpl ) 
    215             ! albedo depends on cloud fraction because of non-linear spectral effects 
    216             alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    217                                  CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    218             IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     207 
     208            CASE( jp_clio )                                       ! CLIO bulk formulation 
     209               ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     210               ! (alb_ice) is computed within the bulk routine 
     211                                      CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
     212               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     213               IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     214                
     215            CASE( jp_core )                                       ! CORE bulk formulation 
     216               ! albedo depends on cloud fraction because of non-linear spectral effects 
     217               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     218                                      CALL blk_ice_core_flx( t_su, alb_ice ) 
     219               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     220               IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     221                
     222            CASE ( jp_purecpl )                                    ! Coupled formulation 
     223               ! albedo depends on cloud fraction because of non-linear spectral effects 
     224               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     225                                      CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     226               IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     227 
    219228         END SELECT 
     229 
    220230         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    221231 
     
    223233         ! --- ice thermodynamics --- ! 
    224234         !----------------------------! 
    225          CALL lim_thd( kt )                         ! Ice thermodynamics       
    226          ! 
    227          CALL lim_update2( kt )                     ! Corrections 
    228          ! 
    229          CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes 
    230          ! 
    231          IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs  
    232          ! 
    233          CALL lim_wri( 1 )                          ! Ice outputs  
     235         ! --- zap this if no ice thermo --- ! 
     236         IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics       
     237         IF( ln_limthd )              CALL lim_update2( kt )    ! -- Corrections 
     238         ! --- 
     239                                      CALL lim_var_glo2eqv      ! necessary calls (at least for coupling) 
     240                                      CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling) 
     241                                      ! 
     242                                      CALL lim_sbc_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes 
     243                                      ! 
     244         IF(ln_limdiaout)             CALL lim_diahsb           ! -- Diagnostics and outputs  
     245         ! 
     246                                      CALL lim_wri( 1 )         ! -- Ice outputs  
    234247         ! 
    235248         IF( kt == nit000 .AND. ln_rstart )   & 
    236             &             CALL iom_close( numrir )  ! close input ice restart file 
    237          ! 
    238          IF( lrst_ice )   CALL lim_rst_write( kt )  ! Ice restart file  
    239          ! 
    240          IF( ln_icectl )  CALL lim_ctl( kt )        ! alerts in case of model crash 
     249            &                         CALL iom_close( numrir )  ! close input ice restart file 
     250         ! 
     251         IF( lrst_ice )               CALL lim_rst_write( kt )  ! -- Ice restart file  
     252         ! 
     253         IF( ln_icectl )              CALL lim_ctl( kt )        ! alerts in case of model crash 
    241254         ! 
    242255      ENDIF   ! End sea-ice time step only 
     
    246259      !-------------------------! 
    247260      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 
    248       IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
     261      !    using before instantaneous surf. currents 
     262      IF( ln_limdyn )                 CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) ) 
    249263!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    250264      ! 
     
    264278      !!---------------------------------------------------------------------- 
    265279      IF(lwp) WRITE(numout,*) 
    266       IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
     280      IF(lwp) WRITE(numout,*) 'sbc_lim_init : update ocean surface boudary condition'  
    267281      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
    268282      ! 
     
    279293      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
    280294      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics 
    281       ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
     295      IF( ln_limdyn )   ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
    282296      ! 
    283297      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     
    300314      CALL lim_msh                     ! ice mesh initialization 
    301315      ! 
    302       CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation 
     316      IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation 
    303317      !                                ! Initial sea-ice state 
    304318      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
     
    347361      !!------------------------------------------------------------------- 
    348362      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    349       NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  & 
    350          &                ln_limdyn, rn_amax_n, rn_amax_s, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt   
     363      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   & 
     364         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice   
     365      NAMELIST/namicediag/ ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt   
    351366      !!------------------------------------------------------------------- 
    352367      !                     
     
    360375      IF(lwm) WRITE ( numoni, namicerun ) 
    361376      ! 
     377      REWIND( numnam_ice_ref )              ! Namelist namicediag in reference namelist : Parameters for ice 
     378      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903) 
     379903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp ) 
     380 
     381      REWIND( numnam_ice_cfg )              ! Namelist namicediag in configuration namelist : Parameters for ice 
     382      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 ) 
     383904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp ) 
     384      IF(lwm) WRITE ( numoni, namicediag ) 
    362385      ! 
    363386      IF(lwp) THEN                        ! control print 
    364387         WRITE(numout,*) 
    365          WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
     388         WRITE(numout,*) 'ice_run : ice share~d parameters for dynamics/advection/thermo of sea-ice' 
    366389         WRITE(numout,*) ' ~~~~~~' 
    367390         WRITE(numout,*) '   number of ice  categories                               = ', jpl 
    368391         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i 
    369392         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s 
    370          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
    371393         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n  
    372394         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s 
    373          WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
    374          WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
     395         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)            ln_limthd  = ', ln_limthd 
     396         WRITE(numout,*) '   Ice dynamics       (T) or not (F)            ln_limdyn  = ', ln_limdyn 
     397         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch          nn_limdyn  = ', nn_limdyn 
     398         WRITE(numout,*) '       2: total' 
     399         WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)' 
     400         WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)' 
     401         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)     = ', rn_uice 
     402         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)     = ', rn_vice 
     403         WRITE(numout,*) 
     404         WRITE(numout,*) '...and ice diagnostics' 
     405         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~' 
     406         WRITE(numout,*) '   Diagnose heat/mass/salt budget or not     ln_limdiahsb  = ', ln_limdiahsb 
     407         WRITE(numout,*) '   Output   heat/mass/salt budget or not     ln_limdiaout  = ', ln_limdiaout 
    375408         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_icectl 
    376409         WRITE(numout,*) '   i-index for control prints (ln_icectl=true)             = ', iiceprt 
     
    410443      ! 
    411444      REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice 
    412       READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 903) 
    413 903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 
     445      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905) 
     446905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 
    414447 
    415448      REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice 
    416       READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 904 ) 
    417 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 
     449      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 ) 
     450906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 
    418451      IF(lwm) WRITE ( numoni, namiceitd ) 
    419       ! 
    420452      ! 
    421453      IF(lwp) THEN                        ! control print 
    422454         WRITE(numout,*) 
    423          WRITE(numout,*) 'ice_itd : ice cat distribution' 
    424          WRITE(numout,*) ' ~~~~~~' 
     455         WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution ' 
     456         WRITE(numout,*) '~~~~~~~~~~~~' 
    425457         WRITE(numout,*) '   shape of ice categories distribution                          nn_catbnd = ', nn_catbnd 
    426458         WRITE(numout,*) '   mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean 
     
    430462      !- Thickness categories boundaries  
    431463      !---------------------------------- 
    432       IF(lwp) WRITE(numout,*) 
    433       IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution ' 
    434       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    435  
    436464      hi_max(:) = 0._wp 
    437465 
     
    581609      !!---------------------------------------------------------------------- 
    582610      sfx    (:,:) = 0._wp   ; 
    583       sfx_bri(:,:) = 0._wp   ;  
     611      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp 
    584612      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    585613      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     
    592620      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    593621      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    594       wfx_spr(:,:) = 0._wp   ;    
     622      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp   
    595623       
    596624      hfx_thd(:,:) = 0._wp   ;    
Note: See TracChangeset for help on using the changeset viewer.