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 7309 for branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Ignore:
Timestamp:
2016-11-22T18:43:11+01:00 (8 years ago)
Author:
clem
Message:

first implementations

Location:
branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90

    r6140 r7309  
    44   !! Ocean forcing:  analytical momentum, heat and freshwater forcings 
    55   !!===================================================================== 
    6    !! History :  3.0   ! 2006-06  (G. Madec)  Original code 
    7    !!            3.2   ! 2009-07  (G. Madec)  Style only 
     6   !! History :  3.0   ! 2006-06  (G. Madec)    Original code 
     7   !!            3.2   ! 2009-07  (G. Madec)    Style only 
     8   !!            3.7   ! 2016-10  (C. Rousset)  Add analytic for LIM3 (ana_ice) 
    89   !!---------------------------------------------------------------------- 
    910 
     
    1516   USE dom_oce         ! ocean space and time domain 
    1617   USE sbc_oce         ! Surface boundary condition: ocean fields 
     18   USE sbc_ice         ! Surface boundary condition: ice   fields 
    1719   USE phycst          ! physical constants 
    1820   USE in_out_manager  ! I/O manager 
     
    2022   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2123   USE lib_fortran 
    22  
     24   USE wrk_nemo 
     25#if defined key_lim3 
     26   USE ice, ONLY       : pfrld, a_i_b 
     27   USE limthd_dh       ! for CALL lim_thd_snwblow 
     28#endif 
     29    
    2330   IMPLICIT NONE 
    2431   PRIVATE 
    2532 
    26    PUBLIC   sbc_ana    ! routine called in sbcmod module 
    27    PUBLIC   sbc_gyre   ! routine called in sbcmod module 
     33   PUBLIC   sbc_ana         ! routine called in sbcmod module 
     34   PUBLIC   sbc_gyre        ! routine called in sbcmod module 
     35#if defined key_lim3 
     36   PUBLIC   ana_ice_tau     ! routine called in sbc_ice_lim module 
     37   PUBLIC   ana_ice_flx     ! routine called in sbc_ice_lim module 
     38#endif 
    2839 
    2940   !                       !!* Namelist namsbc_ana * 
    30    INTEGER  ::   nn_tau000  ! nb of time-step during which the surface stress 
    31    !                        ! increase from 0 to its nominal value  
    32    REAL(wp) ::   rn_utau0   ! constant wind stress value in i-direction 
    33    REAL(wp) ::   rn_vtau0   ! constant wind stress value in j-direction 
    34    REAL(wp) ::   rn_qns0    ! non solar heat flux 
    35    REAL(wp) ::   rn_qsr0    !     solar heat flux 
    36    REAL(wp) ::   rn_emp0    ! net freshwater flux 
     41   ! --- oce variables --- ! 
     42   INTEGER  ::   nn_tau000 ! nb of time-step during which the surface stress 
     43   !                       ! increase from 0 to its nominal value  
     44   REAL(wp) ::   rn_utau0  ! constant wind stress value in i-direction 
     45   REAL(wp) ::   rn_vtau0  ! constant wind stress value in j-direction 
     46   REAL(wp) ::   rn_qns0   ! non solar heat flux 
     47   REAL(wp) ::   rn_qsr0   !     solar heat flux 
     48   REAL(wp) ::   rn_emp0   ! net freshwater flux 
     49   ! --- ice variables --- ! 
     50   REAL(wp) ::   rn_iutau0 ! constant wind stress value in i-direction over ice 
     51   REAL(wp) ::   rn_ivtau0 ! constant wind stress value in j-direction over ice 
     52   REAL(wp) ::   rn_iqns0  ! non solar heat flux over ice 
     53   REAL(wp) ::   rn_iqsr0  !     solar heat flux over ice 
     54   REAL(wp) ::   rn_sprec0 ! snow precip 
     55   REAL(wp) ::   rn_ievap0 ! sublimation 
    3756    
    3857   !! * Substitutions 
     
    6887      REAL(wp) ::   zcoef, zty, zmod      !   -      - 
    6988      !! 
    70       NAMELIST/namsbc_ana/ nn_tau000, rn_utau0, rn_vtau0, rn_qns0, rn_qsr0, rn_emp0 
     89      NAMELIST/namsbc_ana/ nn_tau000, rn_utau0, rn_vtau0, rn_qns0, rn_qsr0, rn_emp0,  & 
     90         &                 rn_iutau0, rn_ivtau0, rn_iqsr0, rn_iqns0, rn_sprec0, rn_ievap0 
    7191      !!--------------------------------------------------------------------- 
    7292      ! 
     
    85105         IF(lwp) WRITE(numout,*)' sbc_ana : Constant surface fluxes read in namsbc_ana namelist' 
    86106         IF(lwp) WRITE(numout,*)' ~~~~~~~ ' 
    87          IF(lwp) WRITE(numout,*)'              spin up of the stress  nn_tau000 = ', nn_tau000, ' time-steps' 
    88          IF(lwp) WRITE(numout,*)'              constant i-stress      rn_utau0  = ', rn_utau0 , ' N/m2' 
    89          IF(lwp) WRITE(numout,*)'              constant j-stress      rn_vtau0  = ', rn_vtau0 , ' N/m2' 
    90          IF(lwp) WRITE(numout,*)'              non solar heat flux    rn_qns0   = ', rn_qns0  , ' W/m2' 
    91          IF(lwp) WRITE(numout,*)'              solar heat flux        rn_qsr0   = ', rn_qsr0  , ' W/m2' 
    92          IF(lwp) WRITE(numout,*)'              net heat flux          rn_emp0   = ', rn_emp0  , ' Kg/m2/s' 
     107         IF(lwp) WRITE(numout,*)'              spin up of the stress         nn_tau000 = ', nn_tau000 , ' time-steps' 
     108         IF(lwp) WRITE(numout,*)'              constant i-stress             rn_utau0  = ', rn_utau0  , ' N/m2' 
     109         IF(lwp) WRITE(numout,*)'              constant j-stress             rn_vtau0  = ', rn_vtau0  , ' N/m2' 
     110         IF(lwp) WRITE(numout,*)'              non solar heat flux           rn_qns0   = ', rn_qns0   , ' W/m2' 
     111         IF(lwp) WRITE(numout,*)'              solar heat flux               rn_qsr0   = ', rn_qsr0   , ' W/m2' 
     112         IF(lwp) WRITE(numout,*)'              net freshwater flux           rn_emp0   = ', rn_emp0   , ' Kg/m2/s' 
     113         IF(lwp) WRITE(numout,*)'              constant ice-atm stress       rn_iutau0 = ', rn_iutau0 , ' N/m2' 
     114         IF(lwp) WRITE(numout,*)'              constant ice-atm stress       rn_ivtau0 = ', rn_ivtau0 , ' N/m2' 
     115         IF(lwp) WRITE(numout,*)'              solar heat flux over ice      rn_iqsr0  = ', rn_iqsr0  , ' W/m2' 
     116         IF(lwp) WRITE(numout,*)'              non solar heat flux over ice  rn_iqns0  = ', rn_iqns0  , ' W/m2' 
     117         IF(lwp) WRITE(numout,*)'              snow precip                   rn_sprec0 = ', rn_sprec0 , ' Kg/m2/s' 
     118         IF(lwp) WRITE(numout,*)'              sublimation                   rn_ievap0 = ', rn_ievap0 , ' Kg/m2/s' 
    93119         ! 
    94120         nn_tau000 = MAX( nn_tau000, 1 )     ! must be >= 1 
     
    132158   END SUBROUTINE sbc_ana 
    133159 
    134  
     160#if defined key_lim3 
     161   SUBROUTINE ana_ice_tau 
     162      !!--------------------------------------------------------------------- 
     163      !!                     ***  ROUTINE ana_ice_tau  *** 
     164      !! 
     165      !! ** Purpose :   provide the surface boundary (momentum) condition over sea-ice 
     166      !!--------------------------------------------------------------------- 
     167      utau_ice(:,:) = rn_iutau0 
     168      vtau_ice(:,:) = rn_ivtau0 
     169      
     170   END SUBROUTINE ana_ice_tau 
     171    
     172   SUBROUTINE ana_ice_flx 
     173      !!--------------------------------------------------------------------- 
     174      !!                     ***  ROUTINE ana_ice_flx  *** 
     175      !! 
     176      !! ** Purpose :   provide the surface boundary (flux) condition over sea-ice 
     177      !!--------------------------------------------------------------------- 
     178      REAL(wp), DIMENSION(:,:), POINTER ::   zsnw       ! snw distribution after wind blowing 
     179      !!--------------------------------------------------------------------- 
     180      CALL wrk_alloc( jpi,jpj, zsnw )  
     181 
     182      ! ocean variables (renaming) 
     183      emp_oce (:,:)   = rn_emp0 
     184      qsr_oce (:,:)   = rn_qsr0 
     185      qns_oce (:,:)   = rn_qns0 
     186       
     187      ! ice variables 
     188      alb_ice (:,:,:) = 0.7_wp ! useless 
     189      qsr_ice (:,:,:) = rn_iqsr0 
     190      qns_ice (:,:,:) = rn_iqns0 
     191      sprecip (:,:)   = rn_sprec0 
     192      evap_ice(:,:,:) = rn_ievap0 
     193 
     194      ! ice variables deduced from above 
     195      zsnw(:,:) = 1._wp 
     196      !!CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     197      emp_ice  (:,:)   = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw(:,:) 
     198      emp_oce  (:,:)   = emp_oce(:,:) - sprecip(:,:) * (1._wp - zsnw(:,:) ) 
     199      qevap_ice(:,:,:) =   0._wp 
     200      qprec_ice(:,:)   =   rhosn * ( sst_m(:,:) * cpic - lfus ) * tmask(:,:,1) ! in J/m3 
     201      qemp_oce (:,:)   = - emp_oce(:,:) * sst_m(:,:) * rcp 
     202      qemp_ice (:,:)   =   sprecip(:,:) * zsnw * ( sst_m(:,:) * cpic - lfus ) * tmask(:,:,1) ! solid precip (only) 
     203 
     204      ! total fluxes 
     205      emp_tot (:,:) = emp_ice  + emp_oce  
     206      qns_tot (:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 
     207      qsr_tot (:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
     208 
     209      !-------------------------------------------------------------------- 
     210      ! FRACTIONs of net shortwave radiation which is not absorbed in the 
     211      ! thin surface layer and penetrates inside the ice cover 
     212      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
     213      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     214      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     215 
     216      CALL wrk_dealloc( jpi,jpj, zsnw )  
     217       
     218   END SUBROUTINE ana_ice_flx 
     219#endif 
     220 
     221    
    135222   SUBROUTINE sbc_gyre( kt ) 
    136223      !!--------------------------------------------------------------------- 
  • branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r6813 r7309  
    1616   !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE 
    1717   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk 
     18   !!            4.0  !  2016-06  (C. Rousset) Add new param of drags with sea-ice (Lupkes at al 2012) 
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    3839   USE lib_fortran    ! to use key_nosignedzero 
    3940#if defined key_lim3 
    40    USE ice     , ONLY :   u_ice, v_ice, jpl, pfrld, a_i_b 
     41   USE ice, ONLY      : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 
    4142   USE limthd_dh      ! for CALL lim_thd_snwblow 
    4243#elif defined key_lim2 
    43    USE ice_2   , ONLY :  u_ice, v_ice 
    44    USE par_ice_2      ! LIM-2 parameters 
     44   USE ice_2, ONLY    : u_ice, v_ice 
     45   USE par_ice_2 
    4546#endif 
    4647   ! 
     
    6162   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module 
    6263#endif 
    63    PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
     64   PUBLIC   turb_core_2z         ! routine called in sbcblk_mfs module 
    6465 
    6566   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
     
    7778 
    7879   !                                             !!! CORE bulk parameters 
    79    REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density 
    80    REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air 
    81    REAL(wp), PARAMETER ::   Lv   =    2.5e6       ! latent heat of vaporization 
    82    REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation 
    83    REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant 
    84    REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
    85    REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant 
     80   REAL(wp), PARAMETER ::   rhoa   =    1.22        ! air density 
     81   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air 
     82   REAL(wp), PARAMETER ::   Lv     =    2.5e6       ! latent heat of vaporization 
     83   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
     84   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
     85   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
     86   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    8687 
    8788   !                        !!* Namelist namsbc_core : CORE bulk parameters 
     
    9293   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9394   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
    94  
     95   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     96 
     97   ! 
     98   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     99    
    95100   !! * Substitutions 
    96101#  include "vectopt_loop_substitute.h90" 
     
    102107CONTAINS 
    103108 
     109   INTEGER FUNCTION sbc_blk_core_alloc() 
     110      !!------------------------------------------------------------------- 
     111      !!             ***  ROUTINE sbc_blk_core_alloc (clem) *** 
     112      !!------------------------------------------------------------------- 
     113      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_core_alloc ) 
     114      ! 
     115      IF( lk_mpp                  )   CALL mpp_sum( sbc_blk_core_alloc ) 
     116      IF( sbc_blk_core_alloc /= 0 )   CALL ctl_warn('sbc_blk_core_alloc: failed to allocate arrays') 
     117   END FUNCTION sbc_blk_core_alloc 
     118 
     119    
    104120   SUBROUTINE sbc_blk_core( kt ) 
    105121      !!--------------------------------------------------------------------- 
     
    149165      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 " 
    150166      NAMELIST/namsbc_core/ cn_dir , ln_taudif, rn_pfac, rn_efac, rn_vfac,  & 
    151          &                  sn_wndi, sn_wndj  , sn_humi, sn_qsr ,           & 
    152          &                  sn_qlw , sn_tair  , sn_prec, sn_snow,           & 
    153          &                  sn_tdif, rn_zqt   ,  rn_zu 
     167         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
     168         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
     169         &                  sn_tdif, rn_zqt,  rn_zu    , ln_Cd_L12 
    154170      !!--------------------------------------------------------------------- 
    155171      ! 
     
    157173      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
    158174         !                                      ! ====================== ! 
     175         ! 
     176         !                                      ! allocate sbc_blk_core array (clem) 
     177         IF( sbc_blk_core_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core : unable to allocate standard arrays' ) 
    159178         ! 
    160179         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters 
     
    321340         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    322341     
     342      Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem) 
     343      
    323344      ! ... tau module, i and j component 
    324345      DO jj = 1, jpj 
     
    439460      !!--------------------------------------------------------------------- 
    440461      INTEGER  ::   ji, jj    ! dummy loop indices 
    441       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2 
    442462      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
    443463      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     464      REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau) 
    444465      !!--------------------------------------------------------------------- 
    445466      ! 
    446467      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau') 
    447468      ! 
    448       ! local scalars ( place there for vector optimisation purposes) 
    449       zcoef_wnorm  = rhoa * Cice 
    450       zcoef_wnorm2 = rhoa * Cice * 0.5 
     469      CALL wrk_alloc( jpi,jpj, Cd ) 
     470 
     471      Cd(:,:) = Cd_ice 
     472       
     473      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     474#if defined key_lim3 
     475      IF( ln_Cd_L12 ) THEN 
     476         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     477      ENDIF 
     478#endif 
    451479 
    452480!!gm brutal.... 
     
    469497               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    470498                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    471                zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     499               zwnorm_f = rhoa * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    472500               ! ... ice stress at I-point 
    473501               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     
    478506               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    479507                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    480                wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     508               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    481509            END DO 
    482510         END DO 
     
    495523         DO jj = 2, jpjm1 
    496524            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    497                utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     525               utau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    498526                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    499                vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     527               vtau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    500528                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    501529            END DO 
     
    511539         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice_core: wndm_ice : ') 
    512540      ENDIF 
     541 
     542      CALL wrk_dealloc( jpi,jpj, Cd ) 
    513543 
    514544      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
     
    542572      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
    543573      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3) 
     574      REAL(wp), DIMENSION(:,:)  , POINTER ::   Cd                ! transfer coefficient for momentum      (tau) 
    544575      !!--------------------------------------------------------------------- 
    545576      ! 
     
    547578      ! 
    548579      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     580      CALL wrk_alloc( jpi,jpj, Cd ) 
     581 
     582      Cd(:,:) = Cd_ice 
     583 
     584      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     585#if defined key_lim3 
     586      IF( ln_Cd_L12 ) THEN 
     587         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     588      ENDIF 
     589#endif 
    549590 
    550591      ! local scalars ( place there for vector optimisation purposes) 
    551592      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    552       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    553       zcoef_dqsb   = rhoa * cpa * Cice 
     593      zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
     594      zcoef_dqsb   = rhoa * cpa 
    554595 
    555596      zztmp = 1. / ( 1. - albo ) 
     
    577618               ! ... turbulent heat fluxes 
    578619               ! Sensible Heat 
    579                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     620               z_qsb(ji,jj,jl) = rhoa * cpa * Cd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    580621               ! Latent Heat 
    581                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     622               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cd(ji,jj) * wndm_ice(ji,jj)   &                            
    582623                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    583624              ! Latent heat sensitivity for ice (Dqla/Dt) 
    584625               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    585                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
     626                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Cd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
    586627               ELSE 
    587628                  dqla_ice(ji,jj,jl) = 0._wp 
     
    589630 
    590631               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    591                z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj) 
     632               z_dqsb(ji,jj,jl) = zcoef_dqsb * Cd(ji,jj) * wndm_ice(ji,jj) 
    592633 
    593634               ! ----------------------------! 
     
    668709 
    669710      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
     711      CALL wrk_dealloc( jpi,jpj, Cd ) 
    670712      ! 
    671713      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
     
    905947   END FUNCTION psi_h 
    906948 
     949 
     950#if defined key_lim3 
     951   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     952      !!---------------------------------------------------------------------- 
     953      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
     954      !! 
     955      !! ** Purpose :    Recompute the ice-atm drag at 10m height to make 
     956      !!                 it dependent on edges at leads, melt ponds and flows. 
     957      !!                 After some approximations, this can be resumed to a dependency 
     958      !!                 on ice concentration. 
     959      !!                 
     960      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
     961      !!                 with the highest level of approximation: level4, eq.(59) 
     962      !!                 The generic drag over a cell partly covered by ice can be re-written as follows: 
     963      !! 
     964      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu 
     965      !! 
     966      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59) 
     967      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59) 
     968      !!                    A is the concentration of ice minus melt ponds (if any) 
     969      !! 
     970      !!                 This new drag has a parabolic shape (as a function of A) starting at 
     971      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     972      !!                 and going down to Cdi(say 1.4e-3) for A=1 
     973      !! 
     974      !!                 It is theoretically applicable to all ice conditions (not only MIZ) 
     975      !!                 => see Lupkes et al (2013) 
     976      !! 
     977      !! ** References : Lupkes et al. JGR 2012 (theory) 
     978      !!                 Lupkes et al. GRL 2013 (application to GCM) 
     979      !! 
     980      !!---------------------------------------------------------------------- 
     981      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     982      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
     983      REAL(wp), PARAMETER ::   znu   = 1._wp 
     984      REAL(wp), PARAMETER ::   zmu   = 1._wp 
     985      REAL(wp), PARAMETER ::   zbeta = 1._wp 
     986      REAL(wp)            ::   zcoef 
     987      !!---------------------------------------------------------------------- 
     988      zcoef = znu + 1._wp / ( 10._wp * zbeta ) 
     989 
     990      ! generic drag over a cell partly covered by ice 
     991      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
     992      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
     993      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
     994 
     995      ! ice-atm drag 
     996      Cd(:,:) = Cd_ice +  &                                                          ! pure ice drag 
     997         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)   ! change due to sea-ice morphology 
     998       
     999   END SUBROUTINE Cdn10_Lupkes2012 
     1000#endif 
     1001    
    9071002   !!====================================================================== 
    9081003END MODULE sbcblk_core 
  • branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6722 r7309  
    168168#  include "vectopt_loop_substitute.h90" 
    169169   !!---------------------------------------------------------------------- 
    170    !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
     170   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    171171   !! $Id$ 
    172172   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    15991599      ENDIF 
    16001600 
    1601       !! clem: we should output qemp_oce and qemp_ice (at least) 
    1602       IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average) 
    1603       !! these diags are not outputed yet 
    1604 !!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average) 
    1605 !!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average) 
    1606 !!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average) 
     1601      ! some more outputs 
     1602      IF( iom_use('hflx_snow_cea') )    CALL iom_put('hflx_snow_cea',   sprecip(:,:) * ( zcptn(:,:) - Lfus ) )                       ! heat flux from snow (cell average) 
     1603      IF( iom_use('hflx_rain_cea') )    CALL iom_put('hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )                 ! heat flux from rain (cell average) 
     1604      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea',sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average) 
     1605      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea',sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) )           ! heat flux from snow (cell average) 
    16071606 
    16081607#else 
  • branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6416 r7309  
    2424   USE ice             ! LIM-3: ice variables 
    2525   USE thd_ice         ! LIM-3: thermodynamical variables 
    26    USE dom_ice         ! LIM-3: ice domain 
    2726   ! 
    2827   USE sbc_oce         ! Surface boundary condition: ocean fields 
     
    3130   USE sbcblk_clio     ! Surface boundary condition: CLIO bulk 
    3231   USE sbccpl          ! Surface boundary condition: coupled interface 
     32   USE sbcana          ! Surface boundary condition: analytic formulation 
    3333   USE albedo          ! ocean & ice albedo 
    3434   ! 
     
    4848   USE limvar          ! Ice variables switch 
    4949   USE limctl          !  
    50    USE limmsh          ! LIM mesh 
    5150   USE limistate       ! LIM initial state 
    5251   USE limthd_sal      ! LIM ice thermodynamics: salinity 
     
    6564   USE bdyice_lim       ! unstructured open boundary data  (bdy_ice_lim routine) 
    6665#endif 
     66# if defined key_agrif 
     67   USE agrif_ice 
     68   USE agrif_lim3_update 
     69   USE agrif_lim3_interp 
     70# endif 
    6771 
    6872   IMPLICIT NONE 
     
    102106      !!--------------------------------------------------------------------- 
    103107      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    104       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED) 
    105       !! 
    106       INTEGER  ::    jl                 ! dummy loop index 
     108      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=1 ANALYTIC, =3 CLIO, =4 CORE, =5 COUPLED) 
     109      !! 
     110      INTEGER  ::   jl                 ! dummy loop index 
    107111      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    108112      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice  
     
    110114 
    111115      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
     116 
     117      ! clem: it is important to initialize agrif_lim3 variables here and not in sbc_lim_init 
     118# if defined key_agrif 
     119      IF( kt == nit000 ) THEN 
     120         IF( .NOT. Agrif_Root() )   CALL Agrif_InitValues_cont_lim3 
     121      ENDIF 
     122# endif 
    112123 
    113124      !-----------------------! 
     
    115126      !-----------------------! 
    116127      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
     128 
     129# if defined key_agrif 
     130         IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1 
     131# endif 
    117132 
    118133         ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean) 
     
    136151         !----------------------------------------------------------------- 
    137152         SELECT CASE( kblk ) 
    138          CASE( jp_clio    )   ;   CALL blk_ice_clio_tau                         ! CLIO bulk formulation             
    139          CASE( jp_core    )   ;   CALL blk_ice_core_tau                         ! CORE bulk formulation 
    140          CASE( jp_purecpl )   ;   CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation 
     153            CASE( jp_ana     )   ;    CALL ana_ice_tau                              ! analytic formulation             
     154            CASE( jp_clio    )   ;    CALL blk_ice_clio_tau                         ! CLIO bulk formulation             
     155            CASE( jp_core    )   ;    CALL blk_ice_core_tau                         ! CORE bulk formulation 
     156            CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation 
    141157         END SELECT 
    142158          
    143          IF( ln_mixcpl) THEN   ! Case of a mixed Bulk/Coupled formulation 
     159         IF( ln_mixcpl) THEN                                                       ! Case of a mixed Bulk/Coupled formulation 
    144160            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice) 
    145             CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 
     161                                      CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 
    146162            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
    147163            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
     
    154170         numit = numit + nn_fsbc                  ! Ice model time step 
    155171         !                                                    
    156          CALL sbc_lim_bef                         ! Store previous ice values 
    157          CALL sbc_lim_diag0                       ! set diag of mass, heat and salt fluxes to 0 
    158          CALL lim_rst_opn( kt )                   ! Open Ice restart file 
    159          ! 
    160          IF( .NOT. lk_c1d ) THEN 
     172                                      CALL sbc_lim_bef         ! Store previous ice values 
     173                                      CALL sbc_lim_diag0       ! set diag of mass, heat and salt fluxes to 0 
     174                                      CALL lim_rst_opn( kt )   ! Open Ice restart file 
     175         ! 
     176         ! --- zap this if no ice dynamics --- ! 
     177         IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN 
    161178            ! 
    162             CALL lim_dyn( kt )                    ! Ice dynamics    ( rheology/dynamics )    
    163             ! 
    164             CALL lim_trp( kt )                    ! Ice transport   ( Advection/diffusion ) 
    165             ! 
    166             IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 
    167             ! 
    168 #if defined key_bdy 
    169             CALL bdy_ice_lim( kt )                ! bdy ice thermo  
    170             IF( ln_icectl )       CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 
    171 #endif 
    172             ! 
    173             CALL lim_update1( kt )                ! Corrections 
     179            IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics 
     180                                      CALL lim_dyn( kt )       !     rheology   
     181            ELSE 
     182               u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity 
     183               v_ice(:,:) = rn_vice * vmask(:,:,1) 
     184            ENDIF 
     185                                      CALL lim_trp( kt )       ! -- Ice transport (Advection/diffusion) 
     186            IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- Mechanical redistribution (ridging/rafting) 
     187               &                      CALL lim_itd_me          
     188            IF( nn_limdyn == 2 )      CALL lim_update1( kt )   ! -- Corrections 
    174189            ! 
    175190         ENDIF 
    176           
     191 
     192         ! --- 
     193#if defined key_agrif 
     194         IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T') 
     195#endif 
     196#if defined key_bdy 
     197         IF( ln_limthd )              CALL bdy_ice_lim( kt )   ! -- bdy ice thermo  
     198#endif 
     199 
    177200         ! previous lead fraction and ice volume for flux calculations 
    178          CALL sbc_lim_bef                         
    179          CALL lim_var_glo2eqv                     ! ht_i and ht_s for ice albedo calculation 
    180          CALL lim_var_agg(1)                      ! at_i for coupling (via pfrld)  
     201                                      CALL sbc_lim_bef                         
     202                                      CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation 
     203                                      CALL lim_var_agg(1)      ! at_i for coupling (via pfrld)  
     204         ! 
    181205         pfrld(:,:)   = 1._wp - at_i(:,:) 
    182206         phicif(:,:)  = vt_i(:,:) 
     
    193217         !---------------------------------------------------------------------------------------- 
    194218         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    195          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    196  
     219          
     220                                      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    197221         SELECT CASE( kblk ) 
    198          CASE( jp_clio )                                       ! CLIO bulk formulation 
    199             ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
    200             ! (alb_ice) is computed within the bulk routine 
    201                                  CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
    202             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    203             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    204          CASE( jp_core )                                       ! CORE bulk formulation 
    205             ! albedo depends on cloud fraction because of non-linear spectral effects 
    206             alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    207                                  CALL blk_ice_core_flx( t_su, alb_ice ) 
    208             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    209             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    210          CASE ( jp_purecpl ) 
    211             ! albedo depends on cloud fraction because of non-linear spectral effects 
    212             alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    213                                  CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    214             IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     222 
     223            CASE( jp_ana )                                        ! analytic formulation 
     224                                      CALL ana_ice_flx 
     225                
     226            CASE( jp_clio )                                       ! CLIO bulk formulation 
     227               ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     228               ! (alb_ice) is computed within the bulk routine 
     229                                      CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
     230               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     231               IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     232                
     233            CASE( jp_core )                                       ! CORE bulk formulation 
     234               ! albedo depends on cloud fraction because of non-linear spectral effects 
     235               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     236                                      CALL blk_ice_core_flx( t_su, alb_ice ) 
     237               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     238               IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     239                
     240            CASE ( jp_purecpl )                                    ! Coupled formulation 
     241               ! albedo depends on cloud fraction because of non-linear spectral effects 
     242               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     243                                      CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     244               IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     245 
    215246         END SELECT 
     247 
    216248         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    217249 
     
    219251         ! --- ice thermodynamics --- ! 
    220252         !----------------------------! 
    221          CALL lim_thd( kt )                         ! Ice thermodynamics       
    222          ! 
    223          CALL lim_update2( kt )                     ! Corrections 
    224          ! 
    225          CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes 
    226          ! 
    227          IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs  
    228          ! 
    229          CALL lim_wri( 1 )                          ! Ice outputs  
     253         ! --- zap this if no ice thermo --- ! 
     254         IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics       
     255         IF( ln_limthd )              CALL lim_update2( kt )    ! -- Corrections 
     256         ! --- 
     257# if defined key_agrif 
     258         IF( .NOT. Agrif_Root() )     CALL agrif_update_lim3( kt ) 
     259# endif 
     260                                      CALL lim_var_glo2eqv      ! necessary calls (at least for coupling) 
     261                                      CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling) 
     262                                      ! 
     263# if defined key_agrif 
     264!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()  ! clem: should be called at the update frequency only (cf agrif_lim3_update) 
     265# endif 
     266                                      CALL lim_sbc_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes 
     267# if defined key_agrif 
     268!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()  ! clem: should be called at the update frequency only (cf agrif_lim3_update) 
     269# endif 
     270         IF( ln_limdiahsb )           CALL lim_diahsb( kt )     ! -- Diagnostics and outputs  
     271         ! 
     272                                      CALL lim_wri( 1 )         ! -- Ice outputs  
    230273         ! 
    231274         IF( kt == nit000 .AND. ln_rstart )   & 
    232             &             CALL iom_close( numrir )  ! close input ice restart file 
    233          ! 
    234          IF( lrst_ice )   CALL lim_rst_write( kt )  ! Ice restart file  
    235          ! 
    236          IF( ln_icectl )  CALL lim_ctl( kt )        ! alerts in case of model crash 
     275            &                         CALL iom_close( numrir )  ! close input ice restart file 
     276         ! 
     277         IF( lrst_ice )               CALL lim_rst_write( kt )  ! -- Ice restart file  
     278         ! 
     279         IF( ln_limctl )              CALL lim_ctl( kt )        ! alerts in case of model crash 
    237280         ! 
    238281      ENDIF   ! End sea-ice time step only 
     
    242285      !-------------------------! 
    243286      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 
    244       IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
     287      !    using before instantaneous surf. currents 
     288      IF( ln_limdyn )                 CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) ) 
    245289!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    246290      ! 
     
    259303      !!---------------------------------------------------------------------- 
    260304      IF(lwp) WRITE(numout,*) 
    261       IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
     305      IF(lwp) WRITE(numout,*) 'sbc_lim_init : update ocean surface boudary condition'  
    262306      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
    263307      ! 
     
    271315      !                                ! Allocate the ice arrays 
    272316      ierr =        ice_alloc        ()      ! ice variables 
    273       ierr = ierr + dom_ice_alloc    ()      ! domain 
    274317      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
    275318      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics 
    276       ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
     319      IF( ln_limdyn )   ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
    277320      ! 
    278321      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    279322      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays') 
    280323      ! 
    281       !                                ! adequation jpk versus ice/snow layers/categories 
    282       IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk )   & 
    283          &      CALL ctl_stop( 'STOP',                          & 
    284          &     'sbc_lim_init: the 3rd dimension of workspace arrays is too small.',   & 
    285          &     'use more ocean levels or less ice/snow layers/categories.' ) 
     324      CALL lim_dyn_init                ! set ice dynamics parameters 
    286325      ! 
    287326      CALL lim_itd_init                ! ice thickness distribution initialization 
     
    293332      CALL lim_thd_sal_init            ! set ice salinity parameters 
    294333      ! 
    295       CALL lim_msh                     ! ice mesh initialization 
    296       ! 
    297       CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation 
     334      IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation 
    298335      !                                ! Initial sea-ice state 
    299336      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
     
    305342         numit = nit000 - 1 
    306343      ENDIF 
    307       CALL lim_var_agg(1) 
     344      CALL lim_var_agg(2) 
    308345      CALL lim_var_glo2eqv 
    309346      ! 
    310347      CALL lim_sbc_init                 ! ice surface boundary condition    
     348      ! 
     349      IF( ln_limdiahsb) CALL lim_diahsb_init  ! initialization for diags 
    311350      ! 
    312351      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
     
    342381      !!------------------------------------------------------------------- 
    343382      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    344       NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  & 
    345          &                ln_limdyn, rn_amax_n, rn_amax_s, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt   
     383      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   & 
     384         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice   
     385      NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt   
    346386      !!------------------------------------------------------------------- 
    347387      !                     
    348388      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice 
    349389      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 
    350 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 
    351       ! 
     390901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 
     391 
    352392      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice 
    353393      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 
    354 902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 
    355       IF(lwm) WRITE( numoni, namicerun ) 
     394902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 
     395      IF(lwm) WRITE ( numoni, namicerun ) 
     396      ! 
     397      REWIND( numnam_ice_ref )              ! Namelist namicediag in reference namelist : Parameters for ice 
     398      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903) 
     399903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp ) 
     400 
     401      REWIND( numnam_ice_cfg )              ! Namelist namicediag in configuration namelist : Parameters for ice 
     402      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 ) 
     403904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp ) 
     404      IF(lwm) WRITE ( numoni, namicediag ) 
    356405      ! 
    357406      IF(lwp) THEN                        ! control print 
     
    362411         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i 
    363412         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s 
    364          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
    365413         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n  
    366414         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s 
    367          WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
    368          WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
    369          WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_icectl 
    370          WRITE(numout,*) '   i-index for control prints (ln_icectl=true)             = ', iiceprt 
    371          WRITE(numout,*) '   j-index for control prints (ln_icectl=true)             = ', jiceprt 
     415         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)            ln_limthd  = ', ln_limthd 
     416         WRITE(numout,*) '   Ice dynamics       (T) or not (F)            ln_limdyn  = ', ln_limdyn 
     417         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch          nn_limdyn  = ', nn_limdyn 
     418         WRITE(numout,*) '       2: total' 
     419         WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)' 
     420         WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)' 
     421         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)     = ', rn_uice 
     422         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)     = ', rn_vice 
     423         WRITE(numout,*) 
     424         WRITE(numout,*) '...and ice diagnostics' 
     425         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~' 
     426         WRITE(numout,*) '   Diagnose online heat/mass/salt budget     ln_limdiachk  = ', ln_limdiachk 
     427         WRITE(numout,*) '   Output          heat/mass/salt budget     ln_limdiahsb  = ', ln_limdiahsb 
     428         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl 
     429         WRITE(numout,*) '   i-index for control prints (ln_limctl=true)             = ', iiceprt 
     430         WRITE(numout,*) '   j-index for control prints (ln_limctl=true)             = ', jiceprt 
    372431      ENDIF 
    373432      ! 
    374433      ! sea-ice timestep and inverse 
    375       rdt_ice   = nn_fsbc * rdt   
     434      rdt_ice   = REAL(nn_fsbc) * rdt   
    376435      r1_rdtice = 1._wp / rdt_ice  
    377436 
     
    381440      ! 
    382441#if defined key_bdy 
    383       IF( lwp .AND. ln_limdiahsb )  CALL ctl_warn('online conservation check activated but it does not work with BDY') 
     442      IF( lwp .AND. ln_limdiachk )  CALL ctl_warn('online conservation check activated but it does not work with BDY') 
    384443#endif 
     444      ! 
     445      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice 
    385446      ! 
    386447   END SUBROUTINE ice_run 
     
    404465      ! 
    405466      REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice 
    406       READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 903) 
    407 903   IF( ios /= 0 )  CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 
    408       ! 
     467      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905) 
     468905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 
     469 
    409470      REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice 
    410       READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 904 ) 
    411 904   IF( ios /= 0 )  CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 
    412       IF(lwm) WRITE( numoni, namiceitd ) 
     471      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 ) 
     472906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 
     473      IF(lwm) WRITE ( numoni, namiceitd ) 
    413474      ! 
    414475      IF(lwp) THEN                        ! control print 
    415476         WRITE(numout,*) 
    416          WRITE(numout,*) 'ice_itd : ice cat distribution' 
    417          WRITE(numout,*) ' ~~~~~~' 
    418          WRITE(numout,*) '   shape of ice categories distribution                     nn_catbnd = ', nn_catbnd 
    419          WRITE(numout,*) '   mean ice thickness in the domain (used if nn_catbnd=2)  rn_himean = ', rn_himean 
     477         WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution ' 
     478         WRITE(numout,*) '~~~~~~~~~~~~' 
     479         WRITE(numout,*) '   shape of ice categories distribution                          nn_catbnd = ', nn_catbnd 
     480         WRITE(numout,*) '   mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean 
    420481      ENDIF 
    421482      ! 
     
    423484      !- Thickness categories boundaries  
    424485      !---------------------------------- 
    425       IF(lwp) WRITE(numout,*) 
    426       IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution ' 
    427       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    428       ! 
    429486      hi_max(:) = 0._wp 
    430487      ! 
     
    463520 
    464521    
    465    SUBROUTINE ice_lim_flx( ptn_ice , palb_ice, pqns_ice ,    & 
    466       &                    pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 
     522   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 
    467523      !!--------------------------------------------------------------------- 
    468524      !!                  ***  ROUTINE ice_lim_flx  *** 
     
    557613      u_ice_b(:,:)     = u_ice(:,:) 
    558614      v_ice_b(:,:)     = v_ice(:,:) 
    559       !       
     615      at_i_b (:,:)     = SUM( a_i_b(:,:,:), dim=3 ) 
     616       
    560617   END SUBROUTINE sbc_lim_bef 
    561618 
     
    569626      !!---------------------------------------------------------------------- 
    570627      sfx    (:,:) = 0._wp   ; 
    571       sfx_bri(:,:) = 0._wp   ;  
     628      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp 
    572629      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    573630      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     
    580637      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    581638      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    582       wfx_spr(:,:) = 0._wp   ;    
    583       ! 
     639      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp   
     640       
    584641      hfx_thd(:,:) = 0._wp   ;    
    585642      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     
    597654      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp ; 
    598655      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp ; 
    599       ! 
     656 
     657      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here) 
     658       
    600659   END SUBROUTINE sbc_lim_diag0 
    601660 
  • branches/2016/dev_CNRS_AGRIF_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r6460 r7309  
    115115      ! 
    116116      !                          ! overwrite namelist parameter using CPP key information 
    117       IF( Agrif_Root() ) THEN                ! AGRIF zoom 
    118         IF( lk_lim2 )   nn_ice      = 2 
    119         IF( lk_lim3 )   nn_ice      = 3 
    120         IF( lk_cice )   nn_ice      = 4 
    121       ENDIF 
    122       IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration 
    123           ln_ana      = .TRUE.    
    124           nn_ice      =   0 
    125       ENDIF 
    126       ! 
     117#if defined key_agrif 
     118      IF( Agrif_Root() ) THEN                ! AGRIF zoom (cf r1242: possibility to run without ice in fine grid) 
     119         IF( lk_lim2 )   nn_ice      = 2 
     120         IF( lk_lim3 )   nn_ice      = 3 
     121         IF( lk_cice )   nn_ice      = 4 
     122      ENDIF 
     123#else 
     124      IF( lk_lim2 )   nn_ice      = 2 
     125      IF( lk_lim3 )   nn_ice      = 3 
     126      IF( lk_cice )   nn_ice      = 4      
     127#endif 
     128 
     129      IF( cp_cfg == 'gyre' )   ln_ana = .TRUE.          ! GYRE configuration 
     130              
    127131      IF(lwp) THEN               ! Control print 
    128132         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)' 
     
    200204 
    201205      !                                            ! restartability    
    202       IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. ln_cpl ) )   & 
    203          &   CALL ctl_stop( 'LIM sea-ice model requires a bulk formulation or coupled configuration' ) 
    204206      IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. ln_cpl ) )   & 
    205207         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or ln_cpl' ) 
Note: See TracChangeset for help on using the changeset viewer.