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 13284 for NEMO/releases/r4.0/r4.0-HEAD/src/OCE – NEMO

Ignore:
Timestamp:
2020-07-09T17:12:23+02:00 (4 years ago)
Author:
smasson
Message:

4.0-HEAD: merge 4.0-HEAD_r12713_clem_dan_fixcpl into 4.0-HEAD

Location:
NEMO/releases/r4.0/r4.0-HEAD/src/OCE
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/BDY/bdy_oce.F90

    r11536 r13284  
    6363      REAL(wp), POINTER, DIMENSION(:,:) ::  aip    !: now ice  pond concentration 
    6464      REAL(wp), POINTER, DIMENSION(:,:) ::  hip    !: now ice  pond depth 
     65      REAL(wp), POINTER, DIMENSION(:,:) ::  hil    !: now ice  pond lid depth 
    6566#if defined key_top 
    6667      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply 
    … …  
    115116   REAL(wp), DIMENSION(jp_bdy) ::   rice_apnd               !: pond conc.  of incoming sea ice 
    116117   REAL(wp), DIMENSION(jp_bdy) ::   rice_hpnd               !: pond thick. of incoming sea ice 
     118   REAL(wp), DIMENSION(jp_bdy) ::   rice_hlid               !: pond lid thick. of incoming sea ice 
    117119   ! 
    118120   !!---------------------------------------------------------------------- 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/BDY/bdydta.F90

    r13255 r13284  
    4343   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90 
    4444 
    45    INTEGER , PARAMETER ::   jpbdyfld  = 16    ! maximum number of files to read  
     45   INTEGER , PARAMETER ::   jpbdyfld  = 17    ! maximum number of files to read  
    4646   INTEGER , PARAMETER ::   jp_bdyssh = 1     !  
    4747   INTEGER , PARAMETER ::   jp_bdyu2d = 2     !  
    … …  
    6060   INTEGER , PARAMETER ::   jp_bdyaip = 15    !  
    6161   INTEGER , PARAMETER ::   jp_bdyhip = 16    !  
     62   INTEGER , PARAMETER ::   jp_bdyhil = 17    !  
    6263#if ! defined key_si3 
    6364   INTEGER , PARAMETER ::   jpl = 1 
    … …  
    190191                        dta_bdy(jbdy)%aip(ib,jl) =  a_ip(ii,ij,jl) * tmask(ii,ij,1)  
    191192                        dta_bdy(jbdy)%hip(ib,jl) =  h_ip(ii,ij,jl) * tmask(ii,ij,1)  
     193                        dta_bdy(jbdy)%hil(ib,jl) =  h_il(ii,ij,jl) * tmask(ii,ij,1)  
    192194                     END DO 
    193195                  END DO 
    … …  
    299301               &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i ) 
    300302            IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy) 
     303            IF( TRIM(bf_alias(jp_bdyhil)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhil)%fnow(:,1,:) = rice_hlid(jbdy) 
    301304 
    302305            ! if T_i is read and not T_su, set T_su = T_i 
    … …  
    323326               bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp 
    324327               bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp 
     328               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp 
     329            ENDIF 
     330            IF ( .NOT.ln_pnd_lids ) THEN 
     331               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp 
    325332            ENDIF 
    326333             
    … …  
    328335            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)             
    329336            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output) 
    330                CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & 
    331                   &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & 
    332                   &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), & 
    333                   &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), & 
    334                   &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), & 
    335                   &              dta_alias%t_i                  , dta_alias%t_s                  , & 
    336                   &              dta_alias%tsu                  , dta_alias%s_i                  , & 
    337                   &              dta_alias%aip                  , dta_alias%hip ) 
     337               CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & ! in 
     338                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & ! out 
     339                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &                                  ! in (optional) 
     340                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &                                  ! in     - 
     341                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), bf_alias(jp_bdyhil)%fnow(:,1,:), & ! in     - 
     342                  &              dta_alias%t_i                  , dta_alias%t_s                  , &                                  ! out    - 
     343                  &              dta_alias%tsu                  , dta_alias%s_i                  , &                                  ! out    - 
     344                  &              dta_alias%aip                  , dta_alias%hip                  , dta_alias%hil )                    ! out    - 
    338345            ENDIF 
    339346         ENDIF 
    … …  
    379386      !                                                         ! =F => baroclinic velocities in 3D boundary data 
    380387      LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta 
    381       REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd  
     388      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid 
    382389      INTEGER                                ::   ipk,ipl       ! 
    383390      INTEGER                                ::   idvar         ! variable ID 
    … …  
    392399      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill 
    393400      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
    394       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip        
     401      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil        
    395402      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill 
    396403      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias 
    397404      ! 
    398       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d,         & 
    399                          & bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, & 
    400                          & rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd,           & 
     405      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d,                 & 
     406                         & bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip, bn_hil, & 
     407                         & rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid,      & 
    401408                         & ln_full_vel, ln_zinterp 
    402409      !!--------------------------------------------------------------------------- 
    … …  
    455462#if defined key_si3 
    456463         IF( .NOT.ln_pnd ) THEN 
    457             rn_ice_apnd = 0. ; rn_ice_hpnd = 0. 
    458             CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' ) 
     464            rn_ice_apnd = 0. ; rn_ice_hpnd = 0. ; rn_ice_hlid = 0. 
     465            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 & rn_ice_hlid = 0 when no ponds' ) 
     466         ENDIF 
     467         IF( .NOT.ln_pnd_lids ) THEN 
     468            rn_ice_hlid = 0. 
    459469         ENDIF 
    460470#endif 
    … …  
    466476         rice_apnd(jbdy) = rn_ice_apnd 
    467477         rice_hpnd(jbdy) = rn_ice_hpnd 
    468           
     478         rice_hlid(jbdy) = rn_ice_hlid 
     479 
    469480          
    470481         DO jfld = 1, jpbdyfld 
    … …  
    567578            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. & 
    568579               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. & 
    569                & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip     ) THEN 
     580               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip .OR. jfld == jp_bdyhil ) THEN 
    570581               igrd = 1                                                    ! T point 
    571582               ipk = ipl                                                   ! jpl-cat data 
    … …  
    618629               bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy 
    619630               bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta  
     631            ENDIF 
     632            IF( jfld == jp_bdyhil ) THEN 
     633               cl3 = 'hil' 
     634               bf_alias => bf(jp_bdyhil,jbdy:jbdy)                         ! alias for hil structure of bdy number jbdy 
     635               bn_alias => bn_hil                                          ! alias for hil structure of nambdy_dta  
    620636            ENDIF 
    621637 
    … …  
    687703                  ENDIF 
    688704               ENDIF 
     705               IF( jfld == jp_bdyhil ) THEN 
     706                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hil => bf_alias(1)%fnow(:,1,:) 
     707                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hil(iszdim,jpl) ) 
     708                  ENDIF 
     709               ENDIF 
    689710            ENDIF 
    690711 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/BDY/bdyice.F90

    r12520 r13284  
    9494         IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction 
    9595            ! exchange 3d arrays 
    96             CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. & 
    97                  &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. & 
    98                  &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                & 
    99                  &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
     96            CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1.                 & 
     97               &                        , s_i , 'T', 1., t_su, 'T', 1., v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1. & 
     98               &                        , a_ip, 'T', 1., v_ip, 'T', 1., v_il, 'T', 1.                                & 
     99               &                        , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    100100            ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk 
    101101            CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    … …  
    163163            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration 
    164164            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth 
     165            h_il(ji,jj,  jl) = ( h_il(ji,jj,  jl) * zwgt1 + dta%hil(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid depth 
    165166            ! 
    166167            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) 
    … …  
    170171               a_ip(ji,jj,jl) = 0._wp 
    171172               h_ip(ji,jj,jl) = 0._wp 
     173               h_il(ji,jj,jl) = 0._wp 
     174            ENDIF 
     175 
     176            IF( .NOT.ln_pnd_lids ) THEN 
     177               h_il(ji,jj,jl) = 0._wp 
    172178            ENDIF 
    173179            ! 
    … …  
    231237               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) 
    232238               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl) 
     239               h_il(ji,jj,  jl) = h_il(ib,jb,  jl) 
    233240               ! 
    234241               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) 
    … …  
    268275               ! 
    269276               ! melt ponds 
    270                IF( a_i(ji,jj,jl) > epsi10 ) THEN 
    271                   a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl) 
    272                ELSE 
    273                   a_ip_frac(ji,jj,jl) = 0._wp 
    274                ENDIF 
    275277               v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl) 
     278               v_il(ji,jj,jl) = h_il(ji,jj,jl) * a_ip(ji,jj,jl) 
    276279               ! 
    277280            ELSE   ! no ice at the boundary 
    … …  
    281284               h_s (ji,jj,  jl) = 0._wp 
    282285               oa_i(ji,jj,  jl) = 0._wp 
    283                a_ip(ji,jj,  jl) = 0._wp 
    284                v_ip(ji,jj,  jl) = 0._wp 
    285286               t_su(ji,jj,  jl) = rt0 
    286287               t_s (ji,jj,:,jl) = rt0 
    287288               t_i (ji,jj,:,jl) = rt0  
    288289 
    289                a_ip_frac(ji,jj,jl) = 0._wp 
    290                h_ip     (ji,jj,jl) = 0._wp 
    291                a_ip     (ji,jj,jl) = 0._wp 
    292                v_ip     (ji,jj,jl) = 0._wp 
     290               a_ip(ji,jj,jl) = 0._wp 
     291               h_ip(ji,jj,jl) = 0._wp 
     292               h_il(ji,jj,jl) = 0._wp 
    293293                
    294294               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
    … …  
    306306               e_s (ji,jj,:,jl) = 0._wp 
    307307               e_i (ji,jj,:,jl) = 0._wp 
     308               v_ip(ji,jj,  jl) = 0._wp 
     309               v_il(ji,jj,  jl) = 0._wp 
    308310 
    309311            ENDIF 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynnxt.F90

    r12366 r13284  
    3434   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme 
    3535   USE domvvl         ! variable volume 
    36    USE bdy_oce   , ONLY: ln_bdy 
     36   USE bdy_oce , ONLY : ln_bdy 
    3737   USE bdydta         ! ocean open boundary conditions 
    3838   USE bdydyn         ! ocean open boundary conditions 
    … …  
    4848   USE prtctl         ! Print control 
    4949   USE timing         ! Timing 
     50   USE zdfdrg ,  ONLY : ln_drgice_imp, rCdU_top 
    5051#if defined key_agrif 
    5152   USE agrif_oce_interp 
    … …  
    99100      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
    100101      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
     102      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau 
    101103      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
    102104      !!---------------------------------------------------------------------- 
    … …  
    354356      ENDIF 
    355357      ! 
     358      IF ( iom_use("utau") ) THEN 
     359         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 
     360            ALLOCATE(zutau(jpi,jpj))  
     361            DO jj = 2, jpjm1 
     362               DO ji = 2, jpim1 
     363                  jk = miku(ji,jj)  
     364                  zutau(ji,jj) = utau(ji,jj) &  
     365                  &  + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk)  
     366               END DO 
     367            END DO 
     368            CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.) 
     369            CALL iom_put(  "utau", zutau(:,:) ) 
     370            DEALLOCATE(zutau) 
     371         ELSE 
     372            CALL iom_put(  "utau", utau(:,:) ) 
     373         ENDIF 
     374      ENDIF 
     375      ! 
     376      IF ( iom_use("vtau") ) THEN 
     377         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 
     378            ALLOCATE(zvtau(jpi,jpj)) 
     379            DO jj = 2, jpjm1 
     380               DO ji = 2, jpim1 
     381                  jk = mikv(ji,jj) 
     382                  zvtau(ji,jj) = vtau(ji,jj) & 
     383                  &  + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk) 
     384               END DO 
     385            END DO 
     386            CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.) 
     387            CALL iom_put(  "vtau", zvtau(:,:) ) 
     388            DEALLOCATE(zvtau) 
     389         ELSE 
     390            CALL iom_put(  "vtau", vtau(:,:) ) 
     391         ENDIF 
     392      ENDIF 
     393      ! 
    356394      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    357395         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynspg_ts.F90

    r12737 r13284  
    14651465      !                    !==  Set the barotropic drag coef.  ==! 
    14661466      ! 
    1467       IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities) 
     1467      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities) 
    14681468          
    14691469         DO jj = 2, jpjm1 
    … …  
    15281528      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case) 
    15291529      ! 
    1530       IF( ln_isfcav ) THEN 
     1530      IF( ln_isfcav.OR.ln_drgice_imp ) THEN 
    15311531         ! 
    15321532         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynzdf.F90

    r12292 r13284  
    141141            END DO 
    142142         END DO 
    143          IF( ln_isfcav ) THEN    ! Ocean cavities (ISF) 
     143         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF) 
    144144            DO jj = 2, jpjm1         
    145145               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    258258            END DO 
    259259         END DO 
    260          IF ( ln_isfcav ) THEN   ! top friction (always implicit) 
     260         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit) 
    261261            DO jj = 2, jpjm1 
    262262               DO ji = 2, jpim1 
    … …  
    423423            END DO 
    424424         END DO 
    425          IF ( ln_isfcav ) THEN 
     425         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 
    426426            DO jj = 2, jpjm1 
    427427               DO ji = 2, jpim1 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/LBC/lbc_lnk_multi_generic.h90

    r11536 r13284  
    1515#endif 
    1616 
    17    SUBROUTINE ROUTINE_MULTI( cdname                                                                             & 
    18       &                    , pt1, cdna1, psgn1, pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4, cdna4, psgn4   & 
    19       &                    , pt5, cdna5, psgn5, pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8, cdna8, psgn8   & 
    20       &                    , pt9, cdna9, psgn9, pt10, cdna10, psgn10, pt11, cdna11, psgn11                      & 
     17   SUBROUTINE ROUTINE_MULTI( cdname                                                                               & 
     18      &                    , pt1 , cdna1 , psgn1 , pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4 , cdna4 , psgn4   & 
     19      &                    , pt5 , cdna5 , psgn5 , pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8 , cdna8 , psgn8   & 
     20      &                    , pt9 , cdna9 , psgn9 , pt10, cdna10, psgn10, pt11, cdna11, psgn11, pt12, cdna12, psgn12  & 
     21      &                    , pt13, cdna13, psgn13, pt14, cdna14, psgn14, pt15, cdna15, psgn15, pt16, cdna16, psgn16  & 
    2122      &                    , kfillmode, pfillval, lsend, lrecv, ihlcom ) 
    2223      !!--------------------------------------------------------------------- 
    23       CHARACTER(len=*)   ,                   INTENT(in   ) :: cdname  ! name of the calling subroutine 
    24       ARRAY_TYPE(:,:,:,:)          , TARGET, INTENT(inout) :: pt1     ! arrays on which the lbc is applied 
    25       ARRAY_TYPE(:,:,:,:), OPTIONAL, TARGET, INTENT(inout) :: pt2  , pt3  , pt4  , pt5  , pt6  , pt7  , pt8  , pt9  , pt10  , pt11 
    26       CHARACTER(len=1)                     , INTENT(in   ) :: cdna1   ! nature of pt2D. array grid-points 
    27       CHARACTER(len=1)   , OPTIONAL        , INTENT(in   ) :: cdna2, cdna3, cdna4, cdna5, cdna6, cdna7, cdna8, cdna9, cdna10, cdna11 
    28       REAL(wp)                             , INTENT(in   ) :: psgn1   ! sign used across the north fold 
    29       REAL(wp)           , OPTIONAL        , INTENT(in   ) :: psgn2, psgn3, psgn4, psgn5, psgn6, psgn7, psgn8, psgn9, psgn10, psgn11 
    30       INTEGER            , OPTIONAL        , INTENT(in   ) :: kfillmode   ! filling method for halo over land (default = constant) 
    31       REAL(wp)           , OPTIONAL        , INTENT(in   ) :: pfillval    ! background value (used at closed boundaries) 
    32       LOGICAL, DIMENSION(4), OPTIONAL      , INTENT(in   ) :: lsend, lrecv   ! indicate how communications are to be carried out 
    33       INTEGER            , OPTIONAL        , INTENT(in   ) :: ihlcom         ! number of ranks and rows to be communicated 
     24      CHARACTER(len=*)     ,                   INTENT(in   ) ::   cdname  ! name of the calling subroutine 
     25      ARRAY_TYPE(:,:,:,:)            , TARGET, INTENT(inout) ::   pt1     ! arrays on which the lbc is applied 
     26      ARRAY_TYPE(:,:,:,:)  , OPTIONAL, TARGET, INTENT(inout) ::   pt2   , pt3   , pt4   , pt5   , pt6   , pt7   , pt8   , pt9  , & 
     27         &                                                        pt10  , pt11  , pt12  , pt13  , pt14  , pt15  , pt16 
     28      CHARACTER(len=1)                       , INTENT(in   ) ::   cdna1   ! nature of pt2D. array grid-points 
     29      CHARACTER(len=1)     , OPTIONAL        , INTENT(in   ) ::   cdna2 , cdna3 , cdna4 , cdna5 , cdna6 , cdna7 , cdna8 , cdna9, & 
     30         &                                                        cdna10, cdna11, cdna12, cdna13, cdna14, cdna15, cdna16 
     31      REAL(wp)                               , INTENT(in   ) ::   psgn1   ! sign used across the north fold 
     32      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   psgn2 , psgn3 , psgn4 , psgn5 , psgn6 , psgn7 , psgn8 , psgn9, & 
     33         &                                                        psgn10, psgn11, psgn12, psgn13, psgn14, psgn15, psgn16 
     34      INTEGER              , OPTIONAL        , INTENT(in   ) ::   kfillmode   ! filling method for halo over land (default = constant) 
     35      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries) 
     36      LOGICAL, DIMENSION(4), OPTIONAL        , INTENT(in   ) ::   lsend, lrecv   ! indicate how communications are to be carried out 
     37      INTEGER              , OPTIONAL        , INTENT(in   ) ::   ihlcom         ! number of ranks and rows to be communicated 
    3438      !! 
    3539      INTEGER                          ::   kfld        ! number of elements that will be attributed 
    36       PTR_TYPE         , DIMENSION(11) ::   ptab_ptr    ! pointer array 
    37       CHARACTER(len=1) , DIMENSION(11) ::   cdna_ptr    ! nature of ptab_ptr grid-points 
    38       REAL(wp)         , DIMENSION(11) ::   psgn_ptr    ! sign used across the north fold boundary 
     40      PTR_TYPE         , DIMENSION(16) ::   ptab_ptr    ! pointer array 
     41      CHARACTER(len=1) , DIMENSION(16) ::   cdna_ptr    ! nature of ptab_ptr grid-points 
     42      REAL(wp)         , DIMENSION(16) ::   psgn_ptr    ! sign used across the north fold boundary 
    3943      !!--------------------------------------------------------------------- 
    4044      ! 
    … …  
    5559      IF( PRESENT(psgn10) )   CALL ROUTINE_LOAD( pt10, cdna10, psgn10, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
    5660      IF( PRESENT(psgn11) )   CALL ROUTINE_LOAD( pt11, cdna11, psgn11, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     61      IF( PRESENT(psgn12) )   CALL ROUTINE_LOAD( pt12, cdna12, psgn12, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     62      IF( PRESENT(psgn13) )   CALL ROUTINE_LOAD( pt13, cdna13, psgn13, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     63      IF( PRESENT(psgn14) )   CALL ROUTINE_LOAD( pt14, cdna14, psgn14, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     64      IF( PRESENT(psgn15) )   CALL ROUTINE_LOAD( pt15, cdna15, psgn15, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     65      IF( PRESENT(psgn16) )   CALL ROUTINE_LOAD( pt16, cdna16, psgn16, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
    5766      ! 
    58       CALL lbc_lnk_ptr    ( cdname,              ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ihlcom ) 
     67      CALL lbc_lnk_ptr( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv, ihlcom ) 
    5968      ! 
    6069   END SUBROUTINE ROUTINE_MULTI 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbc_ice.F90

    r12395 r13284  
    6969   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s] 
    7070   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s] 
    71    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: wind speed module at T-point                 [m/s] 
     71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: sea surface freezing temperature            [degC] 
     72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rCdU_ice       !: ice-ocean drag at T-point (<0)               [m/s] 
    7273#endif 
    7374 
    … …  
    8990   ! variables used in the coupled interface 
    9091   INTEGER , PUBLIC, PARAMETER ::   jpl = ncat 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj 
     92   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice  
    9293    
    9394   ! already defined in ice.F90 for SI3 
    9495   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i 
    9596   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  h_i, h_s 
     97   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i_last_couple   !: Sea ice fraction on categories at the last coupling point 
    9698 
    9799   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K] 
    98100#endif 
    99101 
    100    REAL(wp), PUBLIC, SAVE ::   cldf_ice = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] 
     102   REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] 
    101103 
    102104   !! arrays relating to embedding ice in the ocean 
    … …  
    131133         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce   (jpi,jpj)     ,   & 
    132134         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce    (jpi,jpj)     ,   & 
    133          &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , STAT= ierr(2) ) 
     135         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , rCdU_ice   (jpi,jpj)     , STAT= ierr(2) ) 
    134136#endif 
    135137 
    … …  
    167169   LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .FALSE.  !: no SI3 ice model 
    168170   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE ice model 
    169    REAL(wp)        , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-] 
     171   REAL(wp)        , PUBLIC, PARAMETER ::   pp_cldf    = 0.81     !: cloud fraction over sea ice, summer CLIO value   [-] 
    170172   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1  
    171173   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbc_oce.F90

    r12132 r13284  
    136136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    137137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
     138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra         !: cloud cover (fraction of cloud in a gridcell) [-] 
    138139 
    139140   !!---------------------------------------------------------------------- 
    … …  
    178179         &      fwficb  (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) ) 
    179180         ! 
    180       ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
    181          &      atm_co2(jpi,jpj) ,                                        & 
    182          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    183          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
     181      ALLOCATE( tprecip(jpi,jpj) , sprecip  (jpi,jpj) , fr_i(jpi,jpj) ,   & 
     182         &      atm_co2(jpi,jpj) , cloud_fra(jpi,jpj) ,                   & 
     183         &      ssu_m  (jpi,jpj) , sst_m    (jpi,jpj) , frq_m(jpi,jpj) ,  & 
     184         &      ssv_m  (jpi,jpj) , sss_m    (jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
    184185         ! 
    185186      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90

    r12926 r13284  
    4646   USE lib_fortran    ! to use key_nosignedzero 
    4747#if defined key_si3 
    48    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif 
    49    USE icethd_dh      ! for CALL ice_thd_snwblow 
     48   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice 
     49   USE icevar         ! for CALL ice_var_snwblow 
    5050#endif 
    5151   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
    … …  
    8080   REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    8181 
    82    INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read 
     82   INTEGER , PARAMETER ::   jpfld   =11           ! maximum number of files to read 
    8383   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    8484   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    … …  
    9090   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s) 
    9191   INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa) 
    92    INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
     92   INTEGER , PARAMETER ::   jp_cc   =10           ! index of cloud cover                     (-)      range:0-1 
     93   INTEGER , PARAMETER ::   jp_tdif =11           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
    9394 
    9495   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
    … …  
    161162      !! 
    162163      !!---------------------------------------------------------------------- 
    163       INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument 
     164      INTEGER  ::   jfpr, jfld            ! dummy loop indice and argument 
    164165      INTEGER  ::   ios, ierror, ioptio   ! Local integer 
    165166      !! 
    … …  
    168169      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    169170      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
    170       TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        " 
     171      TYPE(FLD_N) ::   sn_slp , sn_tdif, sn_cc                 !       "                        " 
    171172      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    172          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
     173         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, sn_cc,         & 
    173174         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    174175         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    … …  
    214215      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    215216      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    216       slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif 
     217      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_cc)   = sn_cc 
     218      slf_i(jp_tdif) = sn_tdif 
    217219      ! 
    218220      lhftau = ln_taudif                     !- add an extra field if HF stress is used 
    … …  
    222224      ALLOCATE( sf(jfld), STAT=ierror ) 
    223225      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 
    224       DO ifpr= 1, jfld 
    225          ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
    226          IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
    227          IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   & 
    228             &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
    229             &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
    230  
    231       END DO 
     226 
    232227      !                                      !- fill the bulk structure with namelist informations 
    233228      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
    234229      ! 
     230      DO jfpr = 1, jfld 
     231         ! 
     232         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero) 
     233            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
     234            sf(jfpr)%fnow(:,:,1) = 0._wp 
     235         ELSE                                                  !-- used field --! 
     236            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
     237            IF( slf_i(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 
     238            IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )                      & 
     239               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 
     240               &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
     241         ENDIF 
     242      ENDDO 
     243      ! fill cloud cover array with constant value if "not used" 
     244      IF( TRIM(sf(jp_cc)%clrootname) == 'NOT USED' )   sf(jp_cc)%fnow(:,:,1) = pp_cldf 
     245          
    235246      IF ( ln_wave ) THEN 
    236247      !Activated wave module but neither drag nor stokes drift activated 
    … …  
    384395      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    385396 
     397      ! --- cloud cover --- ! 
     398      cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1) 
     399 
    386400      ! ----------------------------------------------------------------------------- ! 
    387401      !      0   Wind components and module at T-point relative to the moving ocean   ! 
    … …  
    797811      REAL(wp) ::   zst3                     ! local variable 
    798812      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    799       REAL(wp) ::   zztmp, z1_rLsub           !   -      - 
    800       REAL(wp) ::   zfr1, zfr2               ! local variables 
     813      REAL(wp) ::   zztmp, z1_rLsub          !   -      - 
    801814      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
    802815      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
    … …  
    807820      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
    808821      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2 
     822      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri 
    809823      !!--------------------------------------------------------------------- 
    810824      ! 
    … …  
    881895      ! --- evaporation minus precipitation --- ! 
    882896      zsnw(:,:) = 0._wp 
    883       CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
     897      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
    884898      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    885899      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
    … …  
    908922      END DO 
    909923 
    910       ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! 
    911       zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm 
    912       zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    913       ! 
    914       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
    915          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    916       ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    917          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    918       ELSEWHERE                                                         ! zero when hs>0 
    919          qtr_ice_top(:,:,:) = 0._wp  
    920       END WHERE 
     924      ! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- ! 
     925      IF( nn_qtrice == 0 ) THEN 
     926         ! formulation derived from Grenfell and Maykut (1977), where transmission rate 
     927         !    1) depends on cloudiness 
     928         !    2) is 0 when there is any snow 
     929         !    3) tends to 1 for thin ice 
     930         ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm 
     931         DO jl = 1, jpl 
     932            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm   
     933               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 
     934            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm 
     935               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 
     936            ELSEWHERE                                                         ! zero when hs>0 
     937               qtr_ice_top(:,:,jl) = 0._wp  
     938            END WHERE 
     939         ENDDO 
     940      ELSEIF( nn_qtrice == 1 ) THEN 
     941         ! formulation is derived from the thesis of M. Lebrun (2019). 
     942         !    It represents the best fit using several sets of observations 
     943         !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 
     944         qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:) 
     945      ENDIF 
    921946      ! 
    922947 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk_algo_ncar.F90

    r10190 r13284  
    1111   !! 
    1212   !!       Routine turb_ncar maintained and developed in AeroBulk 
    13    !!                     (http://aerobulk.sourceforge.net/) 
     13   !!                     (https://github.com/brodeau/aerobulk/) 
    1414   !! 
    15    !!                         L. Brodeau, 2015 
     15   !!                         L. Brodeau, 2020 
    1616   !!===================================================================== 
    17    !! History :  3.6  !  2016-02  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90 
     17   !! History :  4.0  !  2020-06  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90 
    1818   !!---------------------------------------------------------------------- 
    1919 
    … …  
    4242   PRIVATE 
    4343 
    44    PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90 
     44   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90 
    4545 
    4646   !                              ! NCAR own values for given constants: 
    4747   REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature... 
    48     
     48 
     49   INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
     50 
    4951   !!---------------------------------------------------------------------- 
    5052CONTAINS 
    5153 
    5254   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    53       &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    54       &                  Cdn, Chn, Cen                       ) 
    55       !!---------------------------------------------------------------------------------- 
     55      &                  Cd, Ch, Ce, t_zu, q_zu, Ub,         & 
     56      &                  CdN, ChN, CeN                       ) 
     57      !!---------------------------------------------------------------------- 
    5658      !!                      ***  ROUTINE  turb_ncar  *** 
    5759      !! 
    … …  
    5961      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008) 
    6062      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    61       !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
    62       !! 
    63       !! ** Method : Monin Obukhov Similarity Theory 
    64       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    65       !! 
    66       !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
    67       !! 
    68       !! ** Last update: Laurent Brodeau, June 2014: 
    69       !!    - handles both cases zt=zu and zt/=zu 
    70       !!    - optimized: less 2D arrays allocated and less operations 
    71       !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    72       !!       rather than temperature difference only... 
    73       !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    74       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    75       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
    76       !!      => 'vkarmn' and 'grav' 
    77       !! 
    78       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     63      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas 
    7964      !! 
    8065      !! INPUT : 
    8166      !! ------- 
    8267      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    83       !!    *  zu   : height for wind speed (generally 10m)                   [m] 
    84       !!    *  U_zu : scalar wind speed at 10m                                [m/s] 
    85       !!    *  sst  : SST                                                     [K] 
     68      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
     69      !!    *  sst  : bulk SST                                                [K] 
    8670      !!    *  t_zt : potential air temperature at zt                         [K] 
    8771      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    8872      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
    89       !! 
     73      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    9074      !! 
    9175      !! OUTPUT : 
    … …  
    9680      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    9781      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    98       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
     82      !!    *  Ub  : bulk wind speed at zu                                 [m/s] 
     83      !! 
     84      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    9985      !!---------------------------------------------------------------------------------- 
    10086      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    … …  
    10389      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin] 
    10490      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg] 
    105       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                   [kg/kg] 
     91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    10692      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
    10793      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    … …  
    11096      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    11197      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    112       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
    113       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114       ! 
    115       INTEGER ::   j_itt 
    116       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    117       INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
     98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ub    ! bulk wind speed at zu                     [m/s] 
     99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   CdN, ChN, CeN ! neutral transfer coefficients 
     100      ! 
     101      INTEGER :: j_itt 
     102      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    118103      ! 
    119104      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient 
    120       REAL(wp), DIMENSION(jpi,jpj) ::   sqrt_Cd_n10   ! root square of Cd_n10 
     105      REAL(wp), DIMENSION(jpi,jpj) ::   sqrtCdN10   ! root square of Cd_n10 
    121106      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu 
    122107      REAL(wp), DIMENSION(jpi,jpj) ::   zpsi_h_u 
    123108      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    124       REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer 
    125       !!---------------------------------------------------------------------------------- 
    126       ! 
    127       l_zt_equal_zu = .FALSE. 
    128       IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    129  
    130       U_blk = MAX( 0.5 , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    131  
    132       !! First guess of stability: 
    133       ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt 
    134       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
    135  
    136       !! Neutral coefficients at 10m: 
     109      REAL(wp), DIMENSION(jpi,jpj) ::   sqrtCd 
     110      !!---------------------------------------------------------------------------------- 
     111 
     112      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) 
     113 
     114      Ub = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
     115 
     116      !! Neutral drag coefficient at zu: 
    137117      IF( ln_cdgw ) THEN      ! wave drag case 
    138          cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
    139          ztmp0   (:,:) = cdn_wave(:,:) 
     118         CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp ) 
    140119      ELSE 
    141          ztmp0 = cd_neutral_10m( U_blk ) 
     120         CdN = CD_N10_NCAR( Ub ) 
    142121      ENDIF 
    143  
    144       sqrt_Cd_n10 = SQRT( ztmp0 ) 
     122      sqrtCdN10 = SQRT( CdN ) 
    145123 
    146124      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    147       Cd = ztmp0 
    148       Ce = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    149       Ch = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    150       stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd) 
    151   
    152       IF( ln_cdgw )   Cen = Ce  ; Chn = Ch 
     125      Cd = CdN 
     126      Ce = CE_N10_NCAR( sqrtCdN10 ) 
     127      ztmp0 = 0.5_wp + SIGN(0.5_wp, virt_temp(t_zt, q_zt) - virt_temp(sst, ssq)) ! we guess stability based on delta of virt. pot. temp. 
     128      Ch = CH_N10_NCAR( sqrtCdN10 , ztmp0 ) 
     129      sqrtCd = sqrtCdN10 
    153130 
    154131      !! Initializing values at z_u with z_t values: 
    155       t_zu = t_zt   ;   q_zu = q_zt 
    156  
    157       !!  * Now starting iteration loop 
    158       DO j_itt=1, nb_itt 
     132      t_zu = t_zt 
     133      q_zu = q_zt 
     134 
     135      !! ITERATION BLOCK 
     136      DO j_itt = 1, nb_itt 
    159137         ! 
    160138         ztmp1 = t_zu - sst   ! Updating air/sea differences 
    161139         ztmp2 = q_zu - ssq 
    162140 
    163          ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    164          ztmp1  = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
    165          ztmp2  = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    166  
    167          ztmp0 = 1. + rctv0*q_zu      ! multiply this with t and you have the virtual temperature 
    168  
    169          ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    170          ztmp0 =  (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 
    171          !                                                      ( Cd*U_blk*U_blk is U*^2 at zu ) 
     141         ! Updating turbulent scales :   (L&Y 2004 Eq. (7)) 
     142         ztmp0 = sqrtCd*Ub          ! u* 
     143         ztmp1 = Ch/sqrtCd*ztmp1    ! theta* 
     144         ztmp2 = Ce/sqrtCd*ztmp2    ! q* 
     145 
     146         ! Estimate the inverse of Obukov length (1/L) at height zu: 
     147         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 
    172148 
    173149         !! Stability parameters : 
    174          zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
    175          zpsi_h_u = psi_h( zeta_u ) 
    176  
    177          !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
     150         zeta_u   = zu*ztmp0 
     151         zeta_u   = sign( min(abs(zeta_u),10._wp), zeta_u ) 
     152 
     153         !! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c)) 
    178154         IF( .NOT. l_zt_equal_zu ) THEN 
    179             !! Array 'stab' is free for the moment so using it to store 'zeta_t' 
    180             stab = zt*ztmp0 ;  stab = SIGN( MIN(ABS(stab),10.0), stab )  ! Temporaty array stab == zeta_t !!! 
    181             stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again! 
    182             t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b) 
    183             q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c) 
    184             q_zu = max(0., q_zu) 
     155            ztmp0 = zt*ztmp0 ! zeta_t ! 
     156            ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 )  ! Temporaty array ztmp0 == zeta_t !!! 
     157            ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u) - psi_h_ncar(ztmp0)                   ! ztmp0 just used as temp array again! 
     158            t_zu = t_zt - ztmp1/vkarmn*ztmp0    ! ztmp1 is still theta*  L&Y 2004 Eq. (9b) 
     159            !! 
     160            q_zu = q_zt - ztmp2/vkarmn*ztmp0    ! ztmp2 is still q*      L&Y 2004 Eq. (9c) 
     161            q_zu = MAX(0._wp, q_zu) 
    185162         END IF 
    186163 
    187          ztmp2 = psi_m(zeta_u) 
    188          IF( ln_cdgw ) THEN      ! surface wave case 
    189             stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd)) 
    190             Cd   = stab * stab 
    191             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    192             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    193             ztmp1 = 1. + Chn * ztmp0      
    194             Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    195             ztmp1 = 1. + Cen * ztmp0 
    196             Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    197  
     164         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)... 
     165         !   In very rare low-wind conditions, the old way of estimating the 
     166         !   neutral wind speed at 10m leads to a negative value that causes the code 
     167         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     168         ztmp2 = psi_m_ncar(zeta_u) 
     169         ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ub, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u)) 
     170 
     171         IF( ln_cdgw ) THEN      ! wave drag case 
     172            CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp ) 
    198173         ELSE 
    199             ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
    200             !   In very rare low-wind conditions, the old way of estimating the 
    201             !   neutral wind speed at 10m leads to a negative value that causes the code 
    202             !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
    203             ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
    204             ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
    205             Cdn(:,:) = ztmp0 
    206             sqrt_Cd_n10 = sqrt(ztmp0) 
    207  
    208             stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    209             Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
    210             Chn(:,:) = Cx_n10 
    211  
    212             !! Update of transfer coefficients: 
    213             ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
    214             Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    215             stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
    216  
    217             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    218             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    219             ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
    220             Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
    221  
    222             Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
    223             Cen(:,:) = Cx_n10 
    224             ztmp1 = 1. + Cx_n10*ztmp0 
    225             Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    226             ENDIF 
    227          ! 
    228       END DO 
    229       ! 
     174            CdN   = CD_N10_NCAR(ztmp0)                                       ! Cd_n10 
     175         END IF 
     176         sqrtCdN10 = SQRT(CdN) 
     177 
     178         !! Update of transfer coefficients: 
     179         ztmp1  = 1._wp + sqrtCdN10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u)) 
     180         Cd     = MAX( CdN / ( ztmp1*ztmp1 ) , 0.1E-3_wp ) 
     181         sqrtCd = SQRT( Cd ) 
     182 
     183         ztmp0  = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / sqrtCdN10 
     184         ztmp2  = sqrtCd / sqrtCdN10 
     185 
     186         ztmp1  = 0.5_wp + sign(0.5_wp,zeta_u) ! stability flag 
     187         ChN    = CH_N10_NCAR( sqrtCdN10 , ztmp1 ) 
     188         ztmp1  = 1._wp + ChN*ztmp0 
     189         Ch     = MAX( ChN*ztmp2 / ztmp1 , 0.1E-3_wp )   ! L&Y 2004 Eq. (10b) 
     190 
     191         CeN = CE_N10_NCAR( sqrtCdN10 ) 
     192         ztmp1  = 1._wp + CeN*ztmp0 
     193         Ce     = MAX( CeN*ztmp2 / ztmp1 , 0.1E-3_wp )  ! L&Y 2004 Eq. (10c) 
     194 
     195      END DO !DO j_itt = 1, nb_itt 
     196 
    230197   END SUBROUTINE turb_ncar 
    231198 
    232199 
    233    FUNCTION cd_neutral_10m( pw10 ) 
    234       !!----------------------------------------------------------------------------------       
     200   FUNCTION CD_N10_NCAR( pw10 ) 
     201      !!---------------------------------------------------------------------------------- 
    235202      !! Estimate of the neutral drag coefficient at 10m as a function 
    236203      !! of neutral wind  speed at 10m 
    237204      !! 
    238       !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
    239       !! 
    240       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     205      !! Origin: Large & Yeager 2008, Eq. (11) 
     206      !! 
     207      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    241208      !!---------------------------------------------------------------------------------- 
    242209      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s) 
    243       REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m 
     210      REAL(wp), DIMENSION(jpi,jpj)             :: CD_N10_NCAR 
    244211      ! 
    245212      INTEGER  ::     ji, jj     ! dummy loop indices 
    … …  
    255222            ! 
    256223            ! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    257             zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) )   ! If pw10 < 33. => 0, else => 1 
    258             ! 
    259             cd_neutral_10m(ji,jj) = 1.e-3 * ( & 
    260                &       (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind <  33 m/s 
    261                &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s 
    262             ! 
    263             cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6) 
     224            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1 
     225            ! 
     226            CD_N10_NCAR(ji,jj) = 1.e-3_wp * ( & 
     227               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s 
     228               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s 
     229            ! 
     230            CD_N10_NCAR(ji,jj) = MAX( CD_N10_NCAR(ji,jj), 0.1E-3_wp ) 
    264231            ! 
    265232         END DO 
    266233      END DO 
    267234      ! 
    268    END FUNCTION cd_neutral_10m 
    269  
    270  
    271    FUNCTION psi_m( pzeta ) 
     235   END FUNCTION CD_N10_NCAR 
     236 
     237 
     238 
     239   FUNCTION CH_N10_NCAR( psqrtcdn10 , pstab ) 
     240      !!---------------------------------------------------------------------------------- 
     241      !! Estimate of the neutral heat transfer coefficient at 10m      !! 
     242      !! Origin: Large & Yeager 2008, Eq. (9) and (12) 
     243 
     244      !!---------------------------------------------------------------------------------- 
     245      REAL(wp), DIMENSION(jpi,jpj)             :: ch_n10_ncar 
     246      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 
     247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab      ! stable ABL => 1 / unstable ABL => 0 
     248      !!---------------------------------------------------------------------------------- 
     249      ! 
     250      ch_n10_ncar = MAX( 1.e-3_wp * psqrtcdn10*( 18._wp*pstab + 32.7_wp*(1._wp - pstab) )  , 0.1E-3_wp )   ! Eq. (9) & (12) Large & Yeager, 2008 
     251      ! 
     252   END FUNCTION CH_N10_NCAR 
     253 
     254   FUNCTION CE_N10_NCAR( psqrtcdn10 ) 
     255      !!---------------------------------------------------------------------------------- 
     256      !! Estimate of the neutral heat transfer coefficient at 10m      !! 
     257      !! Origin: Large & Yeager 2008, Eq. (9) and (13) 
     258      !!---------------------------------------------------------------------------------- 
     259      REAL(wp), DIMENSION(jpi,jpj)             :: ce_n10_ncar 
     260      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 
     261      !!---------------------------------------------------------------------------------- 
     262      ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , 0.1E-3_wp ) 
     263      ! 
     264   END FUNCTION CE_N10_NCAR 
     265 
     266 
     267   FUNCTION psi_m_ncar( pzeta ) 
    272268      !!---------------------------------------------------------------------------------- 
    273269      !! Universal profile stability function for momentum 
    274       !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    275       !!      
    276       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     270      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 
     271      !! 
     272      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    277273      !!         and L is M-O length 
    278274      !! 
    279       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    280       !!---------------------------------------------------------------------------------- 
    281       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pzeta 
    282       REAL(wp), DIMENSION(jpi,jpj)             ::   psi_m 
    283       ! 
    284       INTEGER  ::   ji, jj         ! dummy loop indices 
    285       REAL(wp) :: zx2, zx, zstab   ! local scalars 
    286       !!---------------------------------------------------------------------------------- 
    287       ! 
     275      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     276      !!---------------------------------------------------------------------------------- 
     277      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar 
     278      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
     279      ! 
     280      INTEGER  ::   ji, jj    ! dummy loop indices 
     281      REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab,  zstab   ! local scalars 
     282      !!---------------------------------------------------------------------------------- 
    288283      DO jj = 1, jpj 
    289284         DO ji = 1, jpi 
    290             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    291             zx2 = MAX ( zx2 , 1. ) 
    292             zx  = SQRT( zx2 ) 
    293             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    294             ! 
    295             psi_m(ji,jj) =        zstab  * (-5.*pzeta(ji,jj))       &          ! Stable 
    296                &          + (1. - zstab) * (2.*LOG((1. + zx)*0.5)   &          ! Unstable 
    297                &               + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5)  !    " 
     285 
     286            zzeta = pzeta(ji,jj) 
     287            ! 
     288            zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 - 16z)^0.5 
     289            zx2 = MAX( zx2 , 1._wp ) 
     290            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25 
     291            zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp )   & 
     292               &            + LOG( (1._wp + zx2)*0.5_wp )   & 
     293               &          - 2._wp*ATAN(zx) + rpi*0.5_wp 
     294            ! 
     295            zpsi_stab = -5._wp*zzeta 
     296            ! 
     297            zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1 
     298            ! 
     299            psi_m_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable 
     300               &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable 
    298301            ! 
    299302         END DO 
    300303      END DO 
    301       ! 
    302    END FUNCTION psi_m 
    303  
    304  
    305    FUNCTION psi_h( pzeta ) 
     304   END FUNCTION psi_m_ncar 
     305 
     306 
     307   FUNCTION psi_h_ncar( pzeta ) 
    306308      !!---------------------------------------------------------------------------------- 
    307309      !! Universal profile stability function for temperature and humidity 
    308       !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    309       !! 
    310       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     310      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 
     311      !! 
     312      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    311313      !!         and L is M-O length 
    312314      !! 
    313       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    314       !!---------------------------------------------------------------------------------- 
     315      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     316      !!---------------------------------------------------------------------------------- 
     317      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar 
    315318      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
    316       REAL(wp), DIMENSION(jpi,jpj)             :: psi_h 
    317       ! 
    318       INTEGER  ::   ji, jj    ! dummy loop indices 
    319       REAL(wp) :: zx2, zstab  ! local scalars 
     319      ! 
     320      INTEGER  ::   ji, jj     ! dummy loop indices 
     321      REAL(wp) :: zzeta, zx2, zpsi_unst, zpsi_stab, zstab  ! local scalars 
    320322      !!---------------------------------------------------------------------------------- 
    321323      ! 
    322324      DO jj = 1, jpj 
    323325         DO ji = 1, jpi 
    324             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    325             zx2 = MAX ( zx2 , 1. ) 
    326             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    327             ! 
    328             psi_h(ji,jj) =         zstab  * (-5.*pzeta(ji,jj))        &  ! Stable 
    329                &           + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5 ))   ! Unstable 
     326            ! 
     327            zzeta = pzeta(ji,jj) 
     328            ! 
     329            zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 -16z)^0.5 
     330            zx2 = MAX( zx2 , 1._wp ) 
     331            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 
     332            ! 
     333            zpsi_stab = -5._wp*zzeta 
     334            ! 
     335            zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1 
     336            ! 
     337            psi_h_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable 
     338               &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable 
    330339            ! 
    331340         END DO 
    332341      END DO 
    333       ! 
    334    END FUNCTION psi_h 
     342   END FUNCTION psi_h_ncar 
     343 
     344 
     345 
     346 
     347   FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi ) 
     348      !!---------------------------------------------------------------------------------- 
     349      !!  Provides the neutral-stability wind speed at 10 m 
     350      !!---------------------------------------------------------------------------------- 
     351      REAL(wp), DIMENSION(jpi,jpj)             :: UN10_from_CD  !: [m/s] 
     352      REAL(wp),                     INTENT(in) :: pzu  !: measurement heigh of bulk wind speed 
     353      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb  !: bulk wind speed at height pzu m   [m/s] 
     354      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd  !: drag coefficient 
     355      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum [] 
     356      !!---------------------------------------------------------------------------------- 
     357      !! Reminder: UN10 = u*/vkarmn * log(10/z0) 
     358      !!     and: u* = sqrt(Cd) * Ub 
     359      !!                                  u*/vkarmn * log(   10   /       z0    ) 
     360      UN10_from_CD(:,:) = SQRT(pCd(:,:))*pUb/vkarmn * LOG( 10._wp / z0_from_Cd( pzu, pCd(:,:), ppsi=ppsi(:,:) ) ) 
     361      !! 
     362   END FUNCTION UN10_from_CD 
     363 
     364 
     365   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 
     366      !!------------------------------------------------------------------------ 
     367      !! 
     368      !! Evaluates the 1./(Obukhov length) from air temperature, 
     369      !! air specific humidity, and frictional scales u*, t* and q* 
     370      !! 
     371      !! Author: L. Brodeau, June 2019 / AeroBulk 
     372      !!         (https://github.com/brodeau/aerobulk/) 
     373      !!------------------------------------------------------------------------ 
     374      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L     !: 1./(Obukhov length) [m^-1] 
     375      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha         !: reference potential temperature of air [K] 
     376      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa          !: reference specific humidity of air   [kg/kg] 
     377      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus          !: u*: friction velocity [m/s] 
     378      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs     !: \theta* and q* friction aka turb. scales for temp. and spec. hum. 
     379      ! 
     380      INTEGER  ::   ji, jj         ! dummy loop indices 
     381      REAL(wp) ::     zqa          ! local scalar 
     382      !!------------------------------------------------------------------- 
     383      ! 
     384      DO jj = 1, jpj 
     385         DO ji = 1, jpi 
     386            ! 
     387            zqa = (1._wp + rctv0*pqa(ji,jj)) 
     388            ! 
     389            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 
     390            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 
     391            !                      or 
     392            !  b/  -u* [ theta*              + 0.61 theta q* ] 
     393            ! 
     394            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 
     395               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
     396            ! 
     397         END DO 
     398      END DO 
     399      ! 
     400      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) 
     401      ! 
     402   END FUNCTION One_on_L 
     403 
     404 
     405   FUNCTION z0_from_Cd( pzu, pCd,  ppsi ) 
     406      REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd        !: roughness length [m] 
     407      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m] 
     408      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd   !: (neutral or non-neutral) drag coefficient [] 
     409      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum [] 
     410      !! 
     411      !! If pCd is the NEUTRAL-STABILITY drag coefficient then ppsi must be 0 or not given 
     412      !! If pCd is the drag coefficient (in stable or unstable conditions) then pssi must be provided 
     413      !!---------------------------------------------------------------------------------- 
     414      IF ( PRESENT(ppsi) ) THEN 
     415         !! Cd provided is the actual Cd (not the neutral-stability CdN) : 
     416         z0_from_Cd = pzu * EXP( - ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, double-checked! 
     417      ELSE 
     418         !! Cd provided is the neutral-stability Cd, aka CdN : 
     419         z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) )            !LB: ok, double-checked! 
     420      END IF 
     421   END FUNCTION z0_from_Cd 
     422 
     423   FUNCTION virt_temp( pta, pqa ) 
     424      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp !: virtual temperature [K] 
     425      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K] 
     426      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air   [kg/kg] 
     427      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 
     428   END FUNCTION virt_temp 
    335429 
    336430   !!====================================================================== 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbccpl.F90

    r13066 r13284  
    4141#endif 
    4242#if defined key_si3 
    43    USE icethd_dh      ! for CALL ice_thd_snwblow 
     43   USE icevar         ! for CALL ice_var_snwblow 
    4444#endif 
    4545   ! 
    … …  
    4848   USE lib_mpp        ! distribued memory computing library 
    4949   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     50 
     51#if defined key_oasis3  
     52   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut  
     53#endif  
    5054 
    5155   IMPLICIT NONE 
    … …  
    152156   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level  
    153157   INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes) 
    154    INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area 
     158   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area fraction 
    155159   INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness 
    156160   INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity 
    … …  
    159163 
    160164   INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent  
     165 
     166#if ! defined key_oasis3  
     167   ! Dummy variables to enable compilation when oasis3 is not being used  
     168   INTEGER                    ::   OASIS_Sent        = -1  
     169   INTEGER                    ::   OASIS_SentOut     = -1  
     170   INTEGER                    ::   OASIS_ToRest      = -1  
     171   INTEGER                    ::   OASIS_ToRestOut   = -1  
     172#endif  
    161173 
    162174   !                                  !!** namelist namsbc_cpl ** 
    … …  
    184196   LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models 
    185197                                          !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     198   LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)  
     199 
    186200   TYPE ::   DYNARR      
    187201      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3    
    … …  
    248262      REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos 
    249263      !! 
    250       NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   &  
     264      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             & 
     265         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  &  
    251266         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  &  
    252          &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   &  
    253          &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   &  
     267         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  &  
     268         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  &  
    254269         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  & 
    255          &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   & 
    256          &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   & 
    257          &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   & 
     270         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  & 
     271         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                & 
     272         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 & 
    258273         &                  sn_rcv_ts_ice 
    259  
    260274      !!--------------------------------------------------------------------- 
    261275      ! 
    … …  
    279293      ENDIF 
    280294      IF( lwp .AND. ln_cpl ) THEN                        ! control print 
     295         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
     296         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     297         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux 
     298         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl 
    281299         WRITE(numout,*)'  received fields (mutiple ice categogies)' 
    282300         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')' 
    … …  
    327345         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor  
    328346         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd  
    329          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    330          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
    331          WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl 
    332347      ENDIF 
    333348 
    … …  
    366381      !  
    367382      ! Vectors: change of sign at north fold ONLY if on the local grid 
    368       IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled 
     383      IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  & 
     384           .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled 
     385      ! 
    369386      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 
    370387       
    … …  
    821838      END SELECT 
    822839 
     840      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office) 
     841#if defined key_si3 || defined key_cice 
     842       a_i_last_couple(:,:,:) = 0._wp 
     843#endif 
    823844      !                                                      ! ------------------------- !  
    824845      !                                                      !      Ice Meltponds        !  
    … …  
    11081129      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient 
    11091130      REAL(wp) ::   zzx, zzy               ! temporary variables 
    1110       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr 
     1131      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra 
    11111132      !!---------------------------------------------------------------------- 
    11121133      ! 
    … …  
    12261247         ENDIF 
    12271248      ENDIF 
    1228  
     1249!!$      !                                                      ! ========================= ! 
     1250!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      ! 
     1251!!$      !                                                      ! ========================= ! 
     1252!!$      cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1) 
     1253!!$      END SELECT 
     1254!!$ 
     1255      zcloud_fra(:,:) = pp_cldf   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm. 
     1256      IF( ln_mixcpl ) THEN 
     1257         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:) 
     1258      ELSE 
     1259         cloud_fra(:,:) = zcloud_fra(:,:) 
     1260      ENDIF 
     1261      !                                                      ! ========================= ! 
    12291262      ! u(v)tau and taum will be modified by ice model 
    12301263      ! -> need to be reset before each call of the ice/fsbc       
    … …  
    16231656      ! 
    16241657      INTEGER  ::   ji, jj, jl   ! dummy loop index 
    1625       REAL(wp) ::   ztri         ! local scalar 
    16261658      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
    16271659      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice 
    16281660      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
     1661      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total 
    16291662      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu 
     1663      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri 
    16301664      !!---------------------------------------------------------------------- 
    16311665      ! 
    … …  
    16471681         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    16481682         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
    1649          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:) 
    16501683      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp 
    16511684         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1) 
    … …  
    16591692 
    16601693#if defined key_si3 
     1694 
     1695      ! --- evaporation over ice (kg/m2/s) --- ! 
     1696      IF (ln_scale_ice_flux) THEN ! typically met-office requirements 
     1697         IF (sn_rcv_emp%clcat == 'yes') THEN 
     1698            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     1699            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp 
     1700            END WHERE 
     1701            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:) 
     1702            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp 
     1703            END WHERE 
     1704         ELSE 
     1705            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:) 
     1706            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp 
     1707            END WHERE 
     1708            zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1709            DO jl = 2, jpl 
     1710               zevap_ice(:,:,jl) = zevap_ice(:,:,1) 
     1711            ENDDO 
     1712         ENDIF 
     1713      ELSE 
     1714         IF (sn_rcv_emp%clcat == 'yes') THEN 
     1715            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl) 
     1716            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:) 
     1717            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp 
     1718            END WHERE 
     1719         ELSE 
     1720            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) 
     1721            zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1722            DO jl = 2, jpl 
     1723               zevap_ice(:,:,jl) = zevap_ice(:,:,1) 
     1724            ENDDO 
     1725         ENDIF 
     1726      ENDIF 
     1727 
     1728      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN 
     1729         ! For conservative case zemp_ice has not been defined yet. Do it now. 
     1730         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:) 
     1731      ENDIF 
     1732 
    16611733      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing) 
    1662       zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw ) 
     1734      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw ) 
    16631735       
    16641736      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- ! 
    … …  
    16671739 
    16681740      ! --- evaporation over ocean (used later for qemp) --- ! 
    1669       zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) 
    1670  
    1671       ! --- evaporation over ice (kg/m2/s) --- ! 
    1672       DO jl=1,jpl 
    1673          IF (sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 
    1674          ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF 
    1675       ENDDO 
     1741      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) 
    16761742 
    16771743      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
    … …  
    17511817!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff 
    17521818!!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf 
    1753       IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving 
    1754       IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs 
    1755       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
    1756       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
    1757       IF ( iom_use('rain') ) CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
    1758       IF ( iom_use('snow_ao_cea') ) CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
    1759       IF ( iom_use('snow_ai_cea') ) CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
    1760       IF ( iom_use('rain_ao_cea') ) CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average) 
    1761       IF ( iom_use('subl_ai_cea') ) CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average) 
    1762       IF ( iom_use('evap_ao_cea') ) CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
    1763          &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 
     1819      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving 
     1820      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs 
     1821      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
     1822      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
     1823      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
     1824      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
     1825      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
     1826      IF( iom_use('rain_ao_cea') )  CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average) 
     1827      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )     ! Sublimation over sea-ice (cell average) 
     1828      IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
     1829         &                                                         - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 
    17641830      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf 
    17651831      ! 
    … …  
    17691835      CASE( 'oce only' )         ! the required field is directly provided 
    17701836         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) 
     1837         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero 
     1838         ! here so the only flux is the ocean only one. 
     1839         zqns_ice(:,:,:) = 0._wp  
    17711840      CASE( 'conservative' )     ! the required fields are directly provided 
    17721841         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
    … …  
    17981867               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    & 
    17991868                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
    1800                   &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1869                  &                                             + pist(:,:,jl) * picefr(:,:) ) ) 
    18011870            END DO 
    18021871         ELSE 
    … …  
    18041873               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    & 
    18051874                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
    1806                   &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1875                  &                                             + pist(:,:,jl) * picefr(:,:) ) ) 
    18071876            END DO 
    18081877         ENDIF 
    … …  
    19081977      CASE( 'oce only' ) 
    19091978         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) ) 
     1979         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero 
     1980         ! here so the only flux is the ocean only one. 
     1981         zqsr_ice(:,:,:) = 0._wp 
    19101982      CASE( 'conservative' ) 
    19111983         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
    … …  
    19932065            ENDDO 
    19942066         ENDIF 
     2067      CASE( 'none' )  
     2068         zdqns_ice(:,:,:) = 0._wp 
    19952069      END SELECT 
    19962070       
    … …  
    20082082      !                                                      ! ========================= ! 
    20092083      CASE ('coupled') 
    2010          IF( ln_mixcpl ) THEN 
    2011             DO jl=1,jpl 
    2012                qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:) 
    2013                qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:) 
    2014             ENDDO 
     2084         IF (ln_scale_ice_flux) THEN 
     2085            WHERE( a_i(:,:,:) > 1.e-10_wp ) 
     2086               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2087               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2088            ELSEWHERE 
     2089               qml_ice(:,:,:) = 0.0_wp 
     2090               qcn_ice(:,:,:) = 0.0_wp 
     2091            END WHERE 
    20152092         ELSE 
    20162093            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
    … …  
    20232100      IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==! 
    20242101         ! 
    2025          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    2026          ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77) 
    2027          ! 
    2028          WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
    2029             zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    2030          ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
    2031             zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri 
    2032          ELSEWHERE                                                         ! zero when hs>0 
    2033             zqtr_ice_top(:,:,:) = 0._wp 
    2034          END WHERE 
     2102         IF( nn_qtrice == 0 ) THEN 
     2103            ! formulation derived from Grenfell and Maykut (1977), where transmission rate 
     2104            !    1) depends on cloudiness 
     2105            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     2106            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm. 
     2107            !    2) is 0 when there is any snow 
     2108            !    3) tends to 1 for thin ice 
     2109            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm 
     2110            DO jl = 1, jpl 
     2111               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     2112                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 
     2113               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
     2114                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:) 
     2115               ELSEWHERE                                                           ! zero when hs>0 
     2116                  zqtr_ice_top(:,:,jl) = 0._wp  
     2117               END WHERE 
     2118            ENDDO 
     2119         ELSEIF( nn_qtrice == 1 ) THEN 
     2120            ! formulation is derived from the thesis of M. Lebrun (2019). 
     2121            !    It represents the best fit using several sets of observations 
     2122            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 
     2123            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:) 
     2124         ENDIF 
    20352125         !      
    20362126      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==! 
    20372127         ! 
    2038          !                    ! ===> here we must receive the qtr_ice_top array from the coupler 
    2039          !                           for now just assume zero (fully opaque ice) 
     2128         !          ! ===> here we must receive the qtr_ice_top array from the coupler 
     2129         !                 for now just assume zero (fully opaque ice) 
    20402130         zqtr_ice_top(:,:,:) = 0._wp 
    20412131         ! 
    … …  
    22312321      ENDIF 
    22322322 
     2323#if defined key_si3 || defined key_cice 
     2324      ! If this coupling was successful then save ice fraction for use between coupling points.  
     2325      ! This is needed for some calculations where the ice fraction at the last coupling point  
     2326      ! is needed.  
     2327      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. &  
     2328         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN  
     2329         IF ( sn_snd_thick%clcat == 'yes' ) THEN  
     2330           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl) 
     2331         ENDIF 
     2332      ENDIF 
     2333#endif 
     2334 
    22332335      IF( ssnd(jps_fice1)%laction ) THEN 
    22342336         SELECT CASE( sn_snd_thick1%clcat ) 
    … …  
    22942396            SELECT CASE( sn_snd_mpnd%clcat )   
    22952397            CASE( 'yes' )   
    2296                ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) 
     2398               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl) 
    22972399               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)   
    22982400            CASE( 'no' )   
    … …  
    23002402               ztmp4(:,:,:) = 0.0   
    23012403               DO jl=1,jpl   
    2302                  ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)   
    2303                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)  
     2404                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl) 
     2405                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
    23042406               ENDDO   
    23052407            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )   

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcmod.F90

    r12276 r13284  
    564564      ENDIF 
    565565      ! 
    566       CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice) 
    567       CALL iom_put( "vtau", vtau )   ! j-wind stress 
    568       ! 
    569566      IF(ln_ctl) THEN         ! print mean trends (used for debugging) 
    570567         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask ) 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfdrg.F90

    r13268 r13284  
    3232   USE lib_mpp        ! distributed memory computing 
    3333   USE prtctl         ! Print control 
     34   USE sbc_oce , ONLY : nn_ice  
    3435 
    3536   IMPLICIT NONE 
    … …  
    4647   LOGICAL          ::   ln_loglayer  ! logarithmic drag: Cd = vkarmn/log(z/z0) 
    4748   LOGICAL , PUBLIC ::   ln_drgimp    ! implicit top/bottom friction flag 
    48  
     49   LOGICAL , PUBLIC ::   ln_drgice_imp ! implicit ice-ocean drag  
    4950   !                                 !!* Namelist namdrg_top & _bot: TOP or BOTTOM coefficient namelist * 
    5051   REAL(wp)         ::   rn_Cd0       !: drag coefficient                                           [ - ] 
    … …  
    231232      INTEGER   ::   ios, ioptio   ! local integers 
    232233      !! 
    233       NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp 
     234      NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp, ln_drgice_imp 
    234235      !!---------------------------------------------------------------------- 
    235236      ! 
    … …  
    254255         WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer 
    255256         WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp 
     257         WRITE(numout,*) '      implicit ice-ocean drag                   ln_drgice_imp  =', ln_drgice_imp 
    256258      ENDIF 
    257259      ! 
    … …  
    264266      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' ) 
    265267      ! 
     268      IF ( ln_drgice_imp.AND.(.NOT.ln_drgimp) ) &  
     269         &                CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires ln_drgimp=T' ) 
     270      ! 
     271      IF ( ln_drgice_imp.AND.( nn_ice /=2 ) ) & 
     272         &  CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires si3' ) 
    266273      ! 
    267274      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor) 
    … …  
    274281      !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities) 
    275282      ! 
    276       IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting 
    277          ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) ) 
     283      IF( ln_isfcav.OR.ln_drgice_imp ) THEN              ! Ocean cavities: top friction setting 
     284         ALLOCATE( rCdU_top(jpi,jpj) ) 
     285      ENDIF 
     286      ! 
     287      IF( ln_isfcav ) THEN 
     288         ALLOCATE( rCd0_top(jpi,jpj)) 
    278289         CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in 
    279290            &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfgls.F90

    r13268 r13284  
    5454   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1) 
    5555   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation 
     56   INTEGER  ::   nn_z0_ice         ! Roughness accounting for sea ice 
    5657   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2) 
    5758   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen 
    … …  
    6263   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
    6364   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
     65   REAL(wp) ::   rn_hsri           ! Ice ocean roughness 
    6466   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)  
    6567 
    … …  
    151153      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves  
    152154      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves) 
     155      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice 
    153156      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before 
    154157      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before 
    … …  
    166169      ustar2_bot (:,:) = 0._wp 
    167170 
     171      SELECT CASE ( nn_z0_ice ) 
     172      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     173      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     174      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     175      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     176      END SELECT 
     177       
    168178      ! Compute surface, top and bottom friction at T-points 
    169179      DO jj = 2, jpjm1              !==  surface ocean friction 
    … …  
    211221      END SELECT 
    212222      ! 
     223      ! adapt roughness where there is sea ice 
     224      zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro 
     225      ! 
    213226      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==! 
    214227         DO jj = 1, jpjm1 
    … …  
    305318      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)  
    306319      ! First level 
    307       en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  ) 
     320      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  ) 
    308321      zd_lw(:,:,1) = en(:,:,1) 
    309322      zd_up(:,:,1) = 0._wp 
    … …  
    311324      !  
    312325      ! One level below 
    313       en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))  & 
    314          &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   ) 
     326      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 
     327         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   ) 
    315328      zd_lw(:,:,2) = 0._wp  
    316329      zd_up(:,:,2) = 0._wp 
    … …  
    321334      ! 
    322335      ! Dirichlet conditions at k=1 
    323       en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  ) 
     336      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  ) 
    324337      zd_lw(:,:,1) = en(:,:,1) 
    325338      zd_up(:,:,1) = 0._wp 
    … …  
    331344      zd_lw(:,:,2) = 0._wp 
    332345      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 
    333       zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
     346      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    334347          &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
    335348!!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf) 
    … …  
    582595         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    583596         zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    584          zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     597         zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) & 
     598            &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    585599         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 
    586600            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 
    … …  
    855869      REAL(wp)::   zcr   ! local scalar 
    856870      !! 
    857       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    858          &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    859          &            rn_crban, rn_charn, rn_frac_hs,        & 
    860          &            nn_bc_surf, nn_bc_bot, nn_z0_met,     & 
     871      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       & 
     872         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   & 
     873         &            rn_crban, rn_charn, rn_frac_hs,              & 
     874         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, & 
    861875         &            nn_stab_func, nn_clos 
    862876      !!---------------------------------------------------------- 
    … …  
    886900         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    887901         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     902         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice 
     903         SELECT CASE( nn_z0_ice ) 
     904         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking' 
     905         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     906         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)' 
     907         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     908         CASE DEFAULT 
     909            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3') 
     910         END SELECT 
    888911         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs 
    889912         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func 
    890913         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    891914         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
     915         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri 
    892916         WRITE(numout,*) 
    893917         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdfphy.F90

    r11536 r13284  
    2828   USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test) 
    2929   USE sbcrnf         ! surface boundary condition: runoff variables 
     30   USE sbc_ice        ! sea ice drag 
    3031#if defined key_agrif 
    3132   USE agrif_oce_interp   ! interpavm 
    … …  
    252253      ENDIF 
    253254      ! 
     255#if defined key_si3 
     256      IF ( ln_drgice_imp) THEN 
     257         IF ( ln_isfcav ) THEN 
     258            rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:) 
     259         ELSE 
     260            rCdU_top(:,:) = rCdU_ice(:,:) 
     261         ENDIF 
     262      ENDIF 
     263#endif 
     264      !  
    254265      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k) 
    255266      ! 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/ZDF/zdftke.F90

    r13268 r13284  
    4646   USE zdfmxl         ! vertical physics: mixed layer 
    4747   ! 
     48#if defined key_si3 
     49   USE ice, ONLY: hm_i, h_i 
     50#endif 
     51#if defined key_cice 
     52   USE sbc_ice, ONLY: h_i 
     53#endif 
    4854   USE in_out_manager ! I/O manager 
    4955   USE iom            ! I/O manager library 
    … …  
    6268   !                      !!** Namelist  namzdf_tke  ** 
    6369   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not 
     70   INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 
     71   REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice 
    6472   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3) 
    6573   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m] 
    … …  
    7482   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
    7583   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
    76    REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4    
    7784   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    7885   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     86   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)    
    7987 
    8088   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
    … …  
    190198      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    191199      ! 
    192       INTEGER ::   ji, jj, jk              ! dummy loop arguments 
     200      INTEGER ::   ji, jj, jk                  ! dummy loop arguments 
    193201      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars 
    194202      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3 
    195203      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient 
    196       REAL(wp) ::   zbbrau, zri                ! local scalars 
    197       REAL(wp) ::   zfact1, zfact2, zfact3     !   -         - 
    198       REAL(wp) ::   ztx2  , zty2  , zcof       !   -         - 
    199       REAL(wp) ::   ztau  , zdif               !   -         - 
    200       REAL(wp) ::   zus   , zwlc  , zind       !   -         - 
    201       REAL(wp) ::   zzd_up, zzd_lw             !   -         - 
     204      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars 
     205      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      - 
     206      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      - 
     207      REAL(wp) ::   ztau  , zdif               !   -      - 
     208      REAL(wp) ::   zus   , zwlc  , zind       !   -      - 
     209      REAL(wp) ::   zzd_up, zzd_lw             !   -      - 
    202210      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
    203       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i 
     211      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3 
    204212      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    205213      !!-------------------------------------------------------------------- 
    206214      ! 
    207       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
    208       zfact1 = -.5_wp * rdt  
    209       zfact2 = 1.5_wp * rdt * rn_ediss 
    210       zfact3 = 0.5_wp       * rn_ediss 
     215      zbbrau  =  rn_ebb / rau0       ! Local constant initialisation 
     216      zbbirau =  3.75_wp / rau0 
     217      zfact1  = -0.5_wp * rdt  
     218      zfact2  =  1.5_wp * rdt * rn_ediss 
     219      zfact3  =  0.5_wp       * rn_ediss 
     220      ! 
     221      ! ice fraction considered for attenuation of langmuir & wave breaking 
     222      SELECT CASE ( nn_eice ) 
     223      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     224      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     225      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     226      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     227      END SELECT 
    211228      ! 
    212229      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    213230      !                     !  Surface/top/bottom boundary condition on tke 
    214231      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    215        
     232      ! 
    216233      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    217234         DO ji = fs_2, fs_jpim1   ! vector opt. 
    218             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     235            en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 
     236               &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 
    219237         END DO 
    220238      END DO 
    … …  
    248266                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    249267                     &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
    250                   en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) & 
     268                  en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &      
    251269                     &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 
    252270               END DO 
    … …  
    286304            DO ji = fs_2, fs_jpim1   ! vector opt. 
    287305               zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    288                zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    289                IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     306               zus3(ji,jj) = ( 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    290307            END DO 
    291308         END DO          
    … …  
    293310            DO jj = 2, jpjm1 
    294311               DO ji = fs_2, fs_jpim1   ! vector opt. 
    295                   IF ( zfr_i(ji,jj) /= 0. ) THEN                
     312                  IF ( zus3(ji,jj) /= 0. ) THEN                
    296313                     ! vertical velocity due to LC    
    297314                     IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    298315                        !                                           ! vertical velocity due to LC 
    299                         zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     316                        zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    300317                        !                                           ! TKE Langmuir circulation source term 
    301                         en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     318                        en(ji,jj,jk) = en(ji,jj,jk) + rdt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    302319                     ENDIF 
    303320                  ENDIF 
    … …  
    399416       
    400417      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    401          DO jk = 2, jpkm1                       ! rn_eice =0 ON below sea-ice, =4 OFF when ice fraction > 0.25 
     418         DO jk = 2, jpkm1                       ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF 
    402419            DO jj = 2, jpjm1 
    403420               DO ji = fs_2, fs_jpim1   ! vector opt. 
    404421                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    405                      &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     422                     &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    406423               END DO 
    407424            END DO 
    … …  
    412429               jk = nmln(ji,jj) 
    413430               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    414                   &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     431                  &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    415432            END DO 
    416433         END DO 
    … …  
    425442                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    426443                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    427                      &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     444                     &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    428445               END DO 
    429446            END DO 
    … …  
    477494      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars 
    478495      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      - 
    479       REAL(wp) ::   zemxl, zemlm, zemlp       !   -      - 
     496      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      - 
    480497      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace 
    481498      !!-------------------------------------------------------------------- 
    … …  
    490507      zmxlm(:,:,:)  = rmxl_min     
    491508      zmxld(:,:,:)  = rmxl_min 
    492       ! 
     509      !  
    493510      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
     511         ! 
    494512         zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    495          DO jj = 2, jpjm1 
     513#if ! defined key_si3 && ! defined key_cice 
     514         DO jj = 2, jpjm1                     ! No sea-ice 
    496515            DO ji = fs_2, fs_jpim1 
    497                zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    498             END DO 
    499          END DO 
    500       ELSE  
     516               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1) 
     517            END DO 
     518         END DO 
     519#else 
     520 
     521         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     522         ! 
     523         CASE( 0 )                      ! No scaling under sea-ice 
     524            DO jj = 2, jpjm1 
     525               DO ji = fs_2, fs_jpim1 
     526                  zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 
     527               END DO 
     528            END DO 
     529            ! 
     530         CASE( 1 )                      ! scaling with constant sea-ice thickness 
     531            DO jj = 2, jpjm1 
     532               DO ji = fs_2, fs_jpim1 
     533                  zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     534                     &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1) 
     535               END DO 
     536            END DO 
     537            ! 
     538         CASE( 2 )                      ! scaling with mean sea-ice thickness 
     539            DO jj = 2, jpjm1 
     540               DO ji = fs_2, fs_jpim1 
     541#if defined key_si3 
     542                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     543                     &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 
     544#elif defined key_cice 
     545                  zmaxice = MAXVAL( h_i(ji,jj,:) ) 
     546                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     547                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
     548#endif 
     549               END DO 
     550            END DO 
     551            ! 
     552         CASE( 3 )                      ! scaling with max sea-ice thickness 
     553            DO jj = 2, jpjm1 
     554               DO ji = fs_2, fs_jpim1 
     555                  zmaxice = MAXVAL( h_i(ji,jj,:) ) 
     556                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     557                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
     558               END DO 
     559            END DO 
     560            ! 
     561         END SELECT 
     562#endif 
     563         ! 
     564         DO jj = 2, jpjm1 
     565            DO ji = fs_2, fs_jpim1 
     566               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 
     567            END DO 
     568         END DO 
     569         ! 
     570      ELSE 
    501571         zmxlm(:,:,1) = rn_mxl0 
    502572      ENDIF 
    … …  
    643713      INTEGER ::   ios 
    644714      !! 
    645       NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,  & 
    646          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,  & 
    647          &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc,      & 
    648          &                 nn_etau , nn_htau  , rn_efr , rn_eice   
     715      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  & 
     716         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  & 
     717         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             & 
     718         &                 nn_pdl  , ln_lc    , rn_lc,                 & 
     719         &                 nn_etau , nn_htau  , rn_efr   , nn_eice   
    649720      !!---------------------------------------------------------------------- 
    650721      ! 
    … …  
    675746         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0 
    676747         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0 
     748         IF( ln_mxl0 ) THEN 
     749            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice 
     750            IF( nn_mxlice == 1 ) & 
     751            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice 
     752            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     753            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice' 
     754            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness' 
     755            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness' 
     756            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness' 
     757            CASE DEFAULT 
     758               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 
     759            END SELECT 
     760         ENDIF 
    677761         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc 
    678762         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc 
    … …  
    680764         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau 
    681765         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr 
    682          WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice 
    683          WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   ' 
     766         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice 
     767         SELECT CASE( nn_eice )  
     768         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking' 
     769         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     770         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)' 
     771         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     772         CASE DEFAULT 
     773            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 
     774         END SELECT       
    684775         IF( .NOT.ln_drg_OFF ) THEN 
    685776            WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.