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 13472 for NEMO/trunk/src/OCE – NEMO

Changeset 13472 for NEMO/trunk/src/OCE


Ignore:
Timestamp:
2020-09-16T15:05:19+02:00 (4 years ago)
Author:
smasson
Message:

trunk: commit changes from r4.0-HEAD from 13284 to 13449, see #2523

Location:
NEMO/trunk/src/OCE
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/BDY/bdy_oce.F90

    r12377 r13472  
    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/trunk/src/OCE/BDY/bdydta.F90

    r13237 r13472  
    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 
     
    187188                        dta_bdy(jbdy)%aip(ib,jl) =  a_ip(ii,ij,jl) * tmask(ii,ij,1)  
    188189                        dta_bdy(jbdy)%hip(ib,jl) =  h_ip(ii,ij,jl) * tmask(ii,ij,1)  
     190                        dta_bdy(jbdy)%hil(ib,jl) =  h_il(ii,ij,jl) * tmask(ii,ij,1)  
    189191                     END DO 
    190192                  END DO 
     
    289291            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' )   bf_alias(jp_bdytsu)%fnow(:,1,:) = rice_tem (jbdy) 
    290292            IF( TRIM(bf_alias(jp_bdys_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdys_i)%fnow(:,1,:) = rice_sal (jbdy) 
    291             IF( TRIM(bf_alias(jp_bdyaip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyaip)%fnow(:,1,:) = rice_apnd(jbdy) * & ! rice_apnd is the pond fraction 
    292                &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i ) 
     293            IF( TRIM(bf_alias(jp_bdyaip)%clrootname) == 'NOT USED' )   &              ! rice_apnd is the pond fraction 
     294               &   bf_alias(jp_bdyaip)%fnow(:,1,:) = rice_apnd(jbdy) * bf_alias(jp_bdya_i)%fnow(:,1,:)   ! ( a_ip = rice_apnd*a_i ) 
    293295            IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy) 
    294              
     296            IF( TRIM(bf_alias(jp_bdyhil)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhil)%fnow(:,1,:) = rice_hlid(jbdy) 
     297 
    295298            ! if T_i is read and not T_su, set T_su = T_i 
    296299            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) & 
     
    316319               bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp 
    317320               bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp 
     321               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp 
     322            ENDIF 
     323            IF ( .NOT.ln_pnd_lids ) THEN 
     324               bf_alias(jp_bdyhil)%fnow(:,1,:) = 0._wp 
    318325            ENDIF 
    319326             
     
    321328            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)             
    322329            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output) 
    323                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,:), & 
    324                   &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & 
    325                   &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), & 
    326                   &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), & 
    327                   &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), & 
    328                   &              dta_alias%t_i                  , dta_alias%t_s                  , & 
    329                   &              dta_alias%tsu                  , dta_alias%s_i                  , & 
    330                   &              dta_alias%aip                  , dta_alias%hip ) 
     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,:), & ! in 
     331                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & ! out 
     332                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), &                                  ! in (optional) 
     333                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), &                                  ! in     - 
     334                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), bf_alias(jp_bdyhil)%fnow(:,1,:), & ! in     - 
     335                  &              dta_alias%t_i                  , dta_alias%t_s                  , &                                  ! out    - 
     336                  &              dta_alias%tsu                  , dta_alias%s_i                  , &                                  ! out    - 
     337                  &              dta_alias%aip                  , dta_alias%hip                  , dta_alias%hil )                    ! out    - 
    331338            ENDIF 
    332339         ENDIF 
     
    374381      !                                                         ! =F => baroclinic velocities in 3D boundary data 
    375382      LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta 
    376       REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd  
     383      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid 
    377384      INTEGER                                ::   ipk,ipl       ! 
    378385      INTEGER                                ::   idvar         ! variable ID 
     
    387394      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill 
    388395      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
    389       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        
     396      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        
    390397      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill 
    391398      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias 
    392399      ! 
    393       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    394       NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip 
    395       NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
    396       NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp 
     400      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d,                 & 
     401                         & 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, & 
     402                         & rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd, rn_ice_hlid,      & 
     403                         & ln_full_vel, ln_zinterp 
    397404      !!--------------------------------------------------------------------------- 
    398405      ! 
     
    464471#if defined key_si3 
    465472         IF( .NOT.ln_pnd ) THEN 
    466             rn_ice_apnd = 0. ; rn_ice_hpnd = 0. 
    467             CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' ) 
     473            rn_ice_apnd = 0. ; rn_ice_hpnd = 0. ; rn_ice_hlid = 0. 
     474            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 & rn_ice_hlid = 0 when no ponds' ) 
     475         ENDIF 
     476         IF( .NOT.ln_pnd_lids ) THEN 
     477            rn_ice_hlid = 0. 
    468478         ENDIF 
    469479#endif 
     
    475485         rice_apnd(jbdy) = rn_ice_apnd 
    476486         rice_hpnd(jbdy) = rn_ice_hpnd 
    477           
     487         rice_hlid(jbdy) = rn_ice_hlid 
     488 
    478489          
    479490         DO jfld = 1, jpbdyfld 
     
    576587            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. & 
    577588               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. & 
    578                & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip     ) THEN 
     589               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip .OR. jfld == jp_bdyhil ) THEN 
    579590               igrd = 1                                                    ! T point 
    580591               ipk = ipl                                                   ! jpl-cat data 
     
    627638               bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy 
    628639               bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta  
     640            ENDIF 
     641            IF( jfld == jp_bdyhil ) THEN 
     642               cl3 = 'hil' 
     643               bf_alias => bf(jp_bdyhil,jbdy:jbdy)                         ! alias for hil structure of bdy number jbdy 
     644               bn_alias => bn_hil                                          ! alias for hil structure of nambdy_dta  
    629645            ENDIF 
    630646 
     
    696712                  ENDIF 
    697713               ENDIF 
     714               IF( jfld == jp_bdyhil ) THEN 
     715                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hil => bf_alias(1)%fnow(:,1,:) 
     716                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hil(iszdim,jpl) ) 
     717                  ENDIF 
     718               ENDIF 
    698719            ENDIF 
    699720 
  • NEMO/trunk/src/OCE/BDY/bdyice.F90

    r13226 r13472  
    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.0_wp, h_i , 'T', 1.0_wp, h_s , 'T', 1.0_wp, oa_i, 'T', 1.0_wp & 
    97                  &                      , a_ip, 'T', 1.0_wp, v_ip, 'T', 1.0_wp, s_i , 'T', 1.0_wp, t_su, 'T', 1.0_wp & 
    98                  &                      , v_i , 'T', 1.0_wp, v_s , 'T', 1.0_wp, sv_i, 'T', 1.0_wp                & 
    99                  &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
     96            CALL lbc_lnk_multi('bdyice', a_i , 'T', 1._wp, h_i , 'T', 1._wp, h_s , 'T', 1._wp, oa_i, 'T', 1._wp                  & 
     97               &                       , s_i , 'T', 1._wp, t_su, 'T', 1._wp, v_i , 'T', 1._wp, v_s , 'T', 1._wp, sv_i, 'T', 1._wp & 
     98               &                       , a_ip, 'T', 1._wp, v_ip, 'T', 1._wp, v_il, 'T', 1._wp                                     & 
     99               &                       , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    100100            ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk 
    101             CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1.0_wp, e_s , 'T', 1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    102             CALL lbc_lnk_multi( 'bdyice', t_i , 'T', 1.0_wp, e_i , 'T', 1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
     101            CALL lbc_lnk_multi('bdyice', t_s , 'T', 1._wp, e_s , 'T', 1._wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
     102            CALL lbc_lnk_multi('bdyice', t_i , 'T', 1._wp, e_i , 'T', 1._wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    103103         END IF 
    104104      END DO   ! ir 
     
    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) 
     
    265272               ! 
    266273               ! melt ponds 
    267                IF( a_i(ji,jj,jl) > epsi10 ) THEN 
    268                   a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl) 
    269                ELSE 
    270                   a_ip_frac(ji,jj,jl) = 0._wp 
    271                ENDIF 
    272274               v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl) 
     275               v_il(ji,jj,jl) = h_il(ji,jj,jl) * a_ip(ji,jj,jl) 
    273276               ! 
    274277            ELSE   ! no ice at the boundary 
     
    278281               h_s (ji,jj,  jl) = 0._wp 
    279282               oa_i(ji,jj,  jl) = 0._wp 
    280                a_ip(ji,jj,  jl) = 0._wp 
    281                v_ip(ji,jj,  jl) = 0._wp 
    282283               t_su(ji,jj,  jl) = rt0 
    283284               t_s (ji,jj,:,jl) = rt0 
    284285               t_i (ji,jj,:,jl) = rt0  
    285286 
    286                a_ip_frac(ji,jj,jl) = 0._wp 
    287                h_ip     (ji,jj,jl) = 0._wp 
    288                a_ip     (ji,jj,jl) = 0._wp 
    289                v_ip     (ji,jj,jl) = 0._wp 
     287               a_ip(ji,jj,jl) = 0._wp 
     288               h_ip(ji,jj,jl) = 0._wp 
     289               h_il(ji,jj,jl) = 0._wp 
    290290                
    291291               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     
    303303               e_s (ji,jj,:,jl) = 0._wp 
    304304               e_i (ji,jj,:,jl) = 0._wp 
     305               v_ip(ji,jj,  jl) = 0._wp 
     306               v_il(ji,jj,  jl) = 0._wp 
    305307 
    306308            ENDIF 
  • NEMO/trunk/src/OCE/CRS/crsfld.F90

    r13295 r13472  
    146146      CALL iom_put( "voces" , zs_crs )   ! vS 
    147147 
    148       IF( iom_use( "eken") ) THEN     !      kinetic energy 
     148      IF( iom_use( "ke") ) THEN     !      kinetic energy 
    149149         z3d(:,:,jk) = 0._wp  
    150150         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     
    159159         ! 
    160160         CALL crs_dom_ope( z3d, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0_wp ) 
    161          CALL iom_put( "eken", zt_crs ) 
     161         CALL iom_put( "ke", zt_crs ) 
    162162      ENDIF 
    163163      !  Horizontal divergence ( following OCE/DYN/divhor.F90 )  
  • NEMO/trunk/src/OCE/DIA/diawri.F90

    r13295 r13472  
    118118      INTEGER ::   ji, jj, jk       ! dummy loop indices 
    119119      INTEGER ::   ikbot            ! local integer 
     120      REAL(wp)::   ze3 
    120121      REAL(wp)::   zztmp , zztmpx   ! local scalar 
    121122      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     
    175176      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    176177      IF ( iom_use("sbt") ) THEN 
    177          DO_2D( 1, 1, 1, 1 ) 
     178         DO_2D( 0, 0, 0, 0 ) 
    178179            ikbot = mbkt(ji,jj) 
    179180            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
     
    185186      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    186187      IF ( iom_use("sbs") ) THEN 
    187          DO_2D( 1, 1, 1, 1 ) 
     188         DO_2D( 0, 0, 0, 0 ) 
    188189            ikbot = mbkt(ji,jj) 
    189190            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
     
    207208            ! 
    208209         END_2D 
    209          CALL lbc_lnk( 'diawri', z2d, 'T', 1.0_wp ) 
    210210         CALL iom_put( "taubot", z2d )            
    211211      ENDIF 
     
    214214      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    215215      IF ( iom_use("sbu") ) THEN 
    216          DO_2D( 1, 1, 1, 1 ) 
     216         DO_2D( 0, 0, 0, 0 ) 
    217217            ikbot = mbku(ji,jj) 
    218218            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     
    224224      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    225225      IF ( iom_use("sbv") ) THEN 
    226          DO_2D( 1, 1, 1, 1 ) 
     226         DO_2D( 0, 0, 0, 0 ) 
    227227            ikbot = mbkv(ji,jj) 
    228228            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
     
    253253      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    254254 
     255      IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN 
     256         z3d(:,:,jpk) = 0. 
     257         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     258            zztmp  = ts(ji,jj,jk,jp_sal,Kmm) 
     259            zztmpx = (ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm)) * r1_e1u(ji-1,jj) 
     260            zztmpy = (ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm)) * r1_e2v(ji,jj-1) 
     261            z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
     262               &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk) 
     263         END_3D 
     264         CALL iom_put( "socegrad2",  z3d )          ! square of module of sal gradient 
     265         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     266            z3d(ji,jj,jk) = SQRT( z3d(ji,jj,jk) ) 
     267         END_3D 
     268         CALL iom_put( "socegrad" ,  z3d )          ! module of sal gradient 
     269      ENDIF 
     270          
    255271      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    256272         DO_2D( 0, 0, 0, 0 ) 
     
    261277               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    262278         END_2D 
    263          CALL lbc_lnk( 'diawri', z2d, 'T', 1.0_wp ) 
    264279         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
    265          z2d(:,:) = SQRT( z2d(:,:) ) 
     280         DO_2D( 0, 0, 0, 0 ) 
     281            z2d(ji,jj) = SQRT( z2d(ji,jj) ) 
     282         END_2D 
    266283         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient 
    267284      ENDIF 
     
    270287      IF( iom_use("heatc") ) THEN 
    271288         z2d(:,:)  = 0._wp  
    272          DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
     289         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    273290            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
    274291         END_3D 
     
    278295      IF( iom_use("saltc") ) THEN 
    279296         z2d(:,:)  = 0._wp  
    280          DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
     297         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    281298            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    282299         END_3D 
     
    284301      ENDIF 
    285302      ! 
    286       IF ( iom_use("eken") ) THEN 
     303      IF( iom_use("salt2c") ) THEN 
     304         z2d(:,:)  = 0._wp  
     305         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     306            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
     307         END_3D 
     308         CALL iom_put( "salt2c", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
     309      ENDIF 
     310      ! 
     311      IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN 
    287312         z3d(:,:,jpk) = 0._wp  
    288313         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    289             zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    290             z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
    291                &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
    292                &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
    293                &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
    294          END_3D 
    295          CALL lbc_lnk( 'diawri', z3d, 'T', 1.0_wp ) 
    296          CALL iom_put( "eken", z3d )                 ! kinetic energy 
     314            zztmpx = 0.5 * ( uu(ji-1,jj  ,jk,Kmm) + uu(ji,jj,jk,Kmm) ) 
     315            zztmpy = 0.5 * ( vv(ji  ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm) ) 
     316            z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy ) 
     317         END_3D 
     318         CALL iom_put( "ke", z3d )                 ! kinetic energy 
     319 
     320         z2d(:,:)  = 0._wp  
     321         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     322            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk) 
     323         END_3D 
     324         CALL iom_put( "ke_int", z2d )   ! vertically integrated kinetic energy 
    297325      ENDIF 
    298326      ! 
    299327      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
     328 
     329      IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN 
     330          
     331         z3d(:,:,jpk) = 0._wp  
     332         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     333            z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
     334               &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
     335         END_3D 
     336         CALL iom_put( "relvor", z3d )                  ! relative vorticity 
     337 
     338         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     339            z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)  
     340         END_3D 
     341         CALL iom_put( "absvor", z3d )                  ! absolute vorticity 
     342 
     343         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     344            ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     345               &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     346            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     347            ELSE                      ;   ze3 = 0._wp 
     348            ENDIF 
     349            z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)  
     350         END_3D 
     351         CALL iom_put( "potvor", z3d )                  ! potential vorticity 
     352 
     353      ENDIF 
    300354      ! 
    301355      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
     
    315369            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 
    316370         END_3D 
    317          CALL lbc_lnk( 'diawri', z2d, 'U', -1.0_wp ) 
    318371         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
    319372      ENDIF 
     
    324377            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 
    325378         END_3D 
    326          CALL lbc_lnk( 'diawri', z2d, 'U', -1.0_wp ) 
    327379         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
    328380      ENDIF 
     
    342394            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 
    343395         END_3D 
    344          CALL lbc_lnk( 'diawri', z2d, 'V', -1.0_wp ) 
    345396         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
    346397      ENDIF 
     
    351402            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 
    352403         END_3D 
    353          CALL lbc_lnk( 'diawri', z2d, 'V', -1.0_wp ) 
    354404         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
    355405      ENDIF 
     
    360410            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
    361411         END_3D 
    362          CALL lbc_lnk( 'diawri', z2d, 'T', -1.0_wp ) 
    363412         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
    364413      ENDIF 
     
    368417            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
    369418         END_3D 
    370          CALL lbc_lnk( 'diawri', z2d, 'T', -1.0_wp ) 
    371419         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
    372420      ENDIF 
  • NEMO/trunk/src/OCE/DOM/domain.F90

    r13458 r13472  
    120120         WRITE(numout,*)     '         cn_cfg = ', TRIM( cn_cfg ), '   nn_cfg = ', nn_cfg 
    121121      ENDIF 
    122       lwxios = .FALSE. 
     122      nn_wxios = 0 
    123123      ln_xios_read = .FALSE. 
    124124      ! 
  • NEMO/trunk/src/OCE/DYN/dynatf.F90

    r13295 r13472  
    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 
     
    5050   USE prtctl         ! Print control 
    5151   USE timing         ! Timing 
     52   USE zdfdrg ,  ONLY : ln_drgice_imp, rCdU_top 
    5253#if defined key_agrif 
    5354   USE agrif_oce_interp 
     
    120121      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
    121122      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld 
     123      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau 
    122124      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
    123125      !!---------------------------------------------------------------------- 
     
    321323      ENDIF 
    322324      ! 
     325      IF ( iom_use("utau") ) THEN 
     326         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 
     327            ALLOCATE(zutau(jpi,jpj))  
     328            DO_2D( 0, 0, 0, 0 ) 
     329               jk = miku(ji,jj)  
     330               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa) 
     331            END_2D 
     332            CALL iom_put(  "utau", zutau(:,:) ) 
     333            DEALLOCATE(zutau) 
     334         ELSE 
     335            CALL iom_put(  "utau", utau(:,:) ) 
     336         ENDIF 
     337      ENDIF 
     338      ! 
     339      IF ( iom_use("vtau") ) THEN 
     340         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 
     341            ALLOCATE(zvtau(jpi,jpj)) 
     342            DO_2D( 0, 0, 0, 0 ) 
     343               jk = mikv(ji,jj) 
     344               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa) 
     345            END_2D 
     346            CALL iom_put(  "vtau", zvtau(:,:) ) 
     347            DEALLOCATE(zvtau) 
     348         ELSE 
     349            CALL iom_put(  "vtau", vtau(:,:) ) 
     350         ENDIF 
     351      ENDIF 
     352      ! 
    323353      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   & 
    324354         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask ) 
  • NEMO/trunk/src/OCE/DYN/dynspg_ts.F90

    r13295 r13472  
    264264         IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg 
    265265            CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy 
    266             DO_2D( 0, 0, 0, 0 ) 
     266            DO_2D( 0, 0, 0, 0 )                                ! SPG with the application of W/D gravity filters 
    267267               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj ,Kmm) )   & 
    268268                  &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
     
    14051405      !                    !==  Set the barotropic drag coef.  ==! 
    14061406      ! 
    1407       IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities) 
     1407      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities) 
    14081408          
    14091409         DO_2D( 0, 0, 0, 0 ) 
     
    14561456      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case) 
    14571457      ! 
    1458       IF( ln_isfcav ) THEN 
     1458      IF( ln_isfcav.OR.ln_drgice_imp ) THEN 
    14591459         ! 
    14601460         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity 
  • NEMO/trunk/src/OCE/DYN/dynzdf.F90

    r13295 r13472  
    141141            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
    142142         END_2D 
    143          IF( ln_isfcav ) THEN    ! Ocean cavities (ISF) 
     143         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF) 
    144144            DO_2D( 0, 0, 0, 0 ) 
    145145               iku = miku(ji,jj)         ! top ocean level at u- and v-points  
     
    247247            zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
    248248         END_2D 
    249          IF ( ln_isfcav ) THEN   ! top friction (always implicit) 
     249         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit) 
    250250            DO_2D( 0, 0, 0, 0 ) 
    251251               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
     
    273273      !----------------------------------------------------------------------- 
    274274      ! 
    275       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     275      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    276276         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    277277      END_3D 
    278278      ! 
    279       DO_2D( 0, 0, 0, 0 ) 
     279      DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    280280         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    & 
    281281            &             + r_vvl   * e3u(ji,jj,1,Kaa)  
     
    287287      END_3D 
    288288      ! 
    289       DO_2D( 0, 0, 0, 0 ) 
     289      DO_2D( 0, 0, 0, 0 )             !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==! 
    290290         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
    291291      END_2D 
     
    329329            END_3D 
    330330         END SELECT 
    331          DO_2D( 0, 0, 0, 0 ) 
     331         DO_2D( 0, 0, 0, 0 )   !* Surface boundary conditions 
    332332            zwi(ji,jj,1) = 0._wp 
    333333            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    & 
     
    385385            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
    386386         END_2D 
    387          IF ( ln_isfcav ) THEN 
     387         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 
    388388            DO_2D( 0, 0, 0, 0 ) 
    389389               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     
    410410      !----------------------------------------------------------------------- 
    411411      ! 
    412       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     412      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    413413         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    414414      END_3D 
    415415      ! 
    416       DO_2D( 0, 0, 0, 0 ) 
     416      DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
    417417         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    & 
    418418            &             + r_vvl   * e3v(ji,jj,1,Kaa)  
     
    424424      END_3D 
    425425      ! 
    426       DO_2D( 0, 0, 0, 0 ) 
     426      DO_2D( 0, 0, 0, 0 )             !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==! 
    427427         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 
    428428      END_2D 
  • NEMO/trunk/src/OCE/IOM/iom.F90

    r13295 r13472  
    350350           rst_file = TRIM(clpath)//TRIM(cn_ocerst_in) 
    351351        ELSE 
    352            rst_file = TRIM(clpath)//'1_'//TRIM(cn_ocerst_in) 
     352           rst_file = TRIM(clpath)//TRIM(Agrif_CFixed())//'_'//TRIM(cn_ocerst_in) 
    353353        ENDIF 
    354354!set name of the restart file and enable available fields 
  • NEMO/trunk/src/OCE/ISF/isfcavmlt.F90

    r13295 r13472  
    136136      !! ** Method     : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux. 
    137137      !!                 From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being : 
    138       !!                   qfwf  = Gammat * Rau0 * Cp * ( Tw - Tfrz ) / Lf  
     138      !!                   qfwf  = Gammat * rho0 * Cp * ( Tw - Tfrz ) / Lf  
    139139      !!                   qhoce = qlat 
    140140      !!                   qhc   = qfwf * Cp * Tfrz 
  • NEMO/trunk/src/OCE/LBC/lbc_lnk_multi_generic.h90

    r13286 r13472  
    3535#endif 
    3636 
    37    SUBROUTINE ROUTINE_MULTI( cdname                                                                             & 
    38       &                    , pt1, cdna1, psgn1, pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4, cdna4, psgn4   & 
    39       &                    , pt5, cdna5, psgn5, pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8, cdna8, psgn8   & 
    40       &                    , pt9, cdna9, psgn9, pt10, cdna10, psgn10, pt11, cdna11, psgn11                      & 
     37   SUBROUTINE ROUTINE_MULTI( cdname                                                                               & 
     38      &                    , pt1 , cdna1 , psgn1 , pt2 , cdna2 , psgn2 , pt3 , cdna3 , psgn3 , pt4 , cdna4 , psgn4   & 
     39      &                    , pt5 , cdna5 , psgn5 , pt6 , cdna6 , psgn6 , pt7 , cdna7 , psgn7 , pt8 , cdna8 , psgn8   & 
     40      &                    , pt9 , cdna9 , psgn9 , pt10, cdna10, psgn10, pt11, cdna11, psgn11, pt12, cdna12, psgn12  & 
     41      &                    , pt13, cdna13, psgn13, pt14, cdna14, psgn14, pt15, cdna15, psgn15, pt16, cdna16, psgn16  & 
    4142      &                    , kfillmode, pfillval, lsend, lrecv ) 
    4243      !!--------------------------------------------------------------------- 
    43       CHARACTER(len=*)   ,                   INTENT(in   ) :: cdname  ! name of the calling subroutine 
    44       ARRAY_TYPE(:,:,:,:)          , TARGET, INTENT(inout) :: pt1     ! arrays on which the lbc is applied 
    45       ARRAY_TYPE(:,:,:,:), OPTIONAL, TARGET, INTENT(inout) :: pt2  , pt3  , pt4  , pt5  , pt6  , pt7  , pt8  , pt9  , pt10  , pt11 
    46       CHARACTER(len=1)                     , INTENT(in   ) :: cdna1   ! nature of pt2D. array grid-points 
    47       CHARACTER(len=1)   , OPTIONAL        , INTENT(in   ) :: cdna2, cdna3, cdna4, cdna5, cdna6, cdna7, cdna8, cdna9, cdna10, cdna11 
    48       REAL(wp)                             , INTENT(in   ) :: psgn1   ! sign used across the north fold 
    49       REAL(wp)           , OPTIONAL        , INTENT(in   ) :: psgn2, psgn3, psgn4, psgn5, psgn6, psgn7, psgn8, psgn9, psgn10, psgn11 
    50       INTEGER            , OPTIONAL        , INTENT(in   ) :: kfillmode   ! filling method for halo over land (default = constant) 
    51       REAL(wp)           , OPTIONAL        , INTENT(in   ) :: pfillval    ! background value (used at closed boundaries) 
    52       LOGICAL, DIMENSION(4), OPTIONAL      , INTENT(in   ) :: lsend, lrecv   ! indicate how communications are to be carried out 
     44      CHARACTER(len=*)     ,                   INTENT(in   ) ::   cdname  ! name of the calling subroutine 
     45      ARRAY_TYPE(:,:,:,:)            , TARGET, INTENT(inout) ::   pt1     ! arrays on which the lbc is applied 
     46      ARRAY_TYPE(:,:,:,:)  , OPTIONAL, TARGET, INTENT(inout) ::   pt2   , pt3   , pt4   , pt5   , pt6   , pt7   , pt8   , pt9  , & 
     47         &                                                        pt10  , pt11  , pt12  , pt13  , pt14  , pt15  , pt16 
     48      CHARACTER(len=1)                       , INTENT(in   ) ::   cdna1   ! nature of pt2D. array grid-points 
     49      CHARACTER(len=1)     , OPTIONAL        , INTENT(in   ) ::   cdna2 , cdna3 , cdna4 , cdna5 , cdna6 , cdna7 , cdna8 , cdna9, & 
     50         &                                                        cdna10, cdna11, cdna12, cdna13, cdna14, cdna15, cdna16 
     51      REAL(wp)                               , INTENT(in   ) ::   psgn1   ! sign used across the north fold 
     52      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   psgn2 , psgn3 , psgn4 , psgn5 , psgn6 , psgn7 , psgn8 , psgn9, & 
     53         &                                                        psgn10, psgn11, psgn12, psgn13, psgn14, psgn15, psgn16 
     54      INTEGER              , OPTIONAL        , INTENT(in   ) ::   kfillmode   ! filling method for halo over land (default = constant) 
     55      REAL(wp)             , OPTIONAL        , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries) 
     56      LOGICAL, DIMENSION(4), OPTIONAL        , INTENT(in   ) ::   lsend, lrecv   ! indicate how communications are to be carried out 
    5357      !! 
    5458      INTEGER                          ::   kfld        ! number of elements that will be attributed 
    55       PTR_TYPE         , DIMENSION(11) ::   ptab_ptr    ! pointer array 
    56       CHARACTER(len=1) , DIMENSION(11) ::   cdna_ptr    ! nature of ptab_ptr grid-points 
    57       REAL(wp)         , DIMENSION(11) ::   psgn_ptr    ! sign used across the north fold boundary 
     59      PTR_TYPE         , DIMENSION(16) ::   ptab_ptr    ! pointer array 
     60      CHARACTER(len=1) , DIMENSION(16) ::   cdna_ptr    ! nature of ptab_ptr grid-points 
     61      REAL(wp)         , DIMENSION(16) ::   psgn_ptr    ! sign used across the north fold boundary 
    5862      !!--------------------------------------------------------------------- 
    5963      ! 
     
    7478      IF( PRESENT(psgn10) )   CALL ROUTINE_LOAD( pt10, cdna10, psgn10, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
    7579      IF( PRESENT(psgn11) )   CALL ROUTINE_LOAD( pt11, cdna11, psgn11, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     80      IF( PRESENT(psgn12) )   CALL ROUTINE_LOAD( pt12, cdna12, psgn12, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     81      IF( PRESENT(psgn13) )   CALL ROUTINE_LOAD( pt13, cdna13, psgn13, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     82      IF( PRESENT(psgn14) )   CALL ROUTINE_LOAD( pt14, cdna14, psgn14, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     83      IF( PRESENT(psgn15) )   CALL ROUTINE_LOAD( pt15, cdna15, psgn15, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
     84      IF( PRESENT(psgn16) )   CALL ROUTINE_LOAD( pt16, cdna16, psgn16, ptab_ptr, cdna_ptr, psgn_ptr, kfld ) 
    7685      ! 
    77       CALL lbc_lnk_ptr    ( cdname,              ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv ) 
     86      CALL lbc_lnk_ptr( cdname, ptab_ptr, cdna_ptr, psgn_ptr, kfld, kfillmode, pfillval, lsend, lrecv ) 
    7887      ! 
    7988   END SUBROUTINE ROUTINE_MULTI 
  • NEMO/trunk/src/OCE/SBC/sbc_ice.F90

    r12396 r13472  
    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 
     
    9899#endif 
    99100 
    100    REAL(wp), PUBLIC, SAVE ::   cldf_ice = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] 
     101   REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] 
    101102 
    102103   !! arrays relating to embedding ice in the ocean 
     
    131132         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce   (jpi,jpj)     ,   & 
    132133         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce    (jpi,jpj)     ,   & 
    133          &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , STAT= ierr(2) ) 
     134         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , rCdU_ice   (jpi,jpj)     , STAT= ierr(2) ) 
    134135#endif 
    135136 
     
    167168   LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .FALSE.  !: no SI3 ice model 
    168169   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   [-] 
     170   REAL(wp)        , PUBLIC, PARAMETER ::   pp_cldf    = 0.81     !: cloud fraction over sea ice, summer CLIO value   [-] 
    170171   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1  
    171172   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj 
  • NEMO/trunk/src/OCE/SBC/sbc_oce.F90

    r13295 r13472  
    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   !!--------------------------------------------------------------------- 
     
    188189      ! 
    189190      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
    190          &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) ,                       & 
     191         &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) , cloud_fra(jpi,jpj),   & 
    191192         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    192193         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
  • NEMO/trunk/src/OCE/SBC/sbcblk.F90

    r13305 r13472  
    4444   USE lib_fortran    ! to use key_nosignedzero 
    4545#if defined key_si3 
    46    USE ice     , ONLY :   jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif 
    47    USE icethd_dh      ! for CALL ice_thd_snwblow 
     46   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice 
     47   USE icevar         ! for CALL ice_var_snwblow 
    4848#endif 
    4949   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
     
    8787   INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component) 
    8888   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
    89    INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 12   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
    90    INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 13   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
    91    INTEGER , PUBLIC, PARAMETER ::   jpfld    = 13   ! maximum number of files to read 
     89   INTEGER , PUBLIC, PARAMETER ::   jp_cc    = 12   ! index of cloud cover                     (-)      range:0-1 
     90   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 13   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
     91   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 14   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
     92   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 14   ! maximum number of files to read 
    9293 
    9394   ! Warning: keep this structure allocatable for Agrif... 
     
    175176      TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        " 
    176177      TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        " 
    177       TYPE(FLD_N) ::   sn_hpgi, sn_hpgj                        !       "                        " 
     178      TYPE(FLD_N) ::   sn_cc, sn_hpgi, sn_hpgj                 !       "                        " 
    178179      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable 
    179180      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    180181         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     & 
    181          &                 sn_hpgi, sn_hpgj,                                          & 
     182         &                 sn_cc, sn_hpgi, sn_hpgj,                                   & 
    182183         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    183184         &                 cn_dir , rn_zqt, rn_zu,                                    & 
     
    260261      slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi 
    261262      slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow 
    262       slf_i(jp_slp  ) = sn_slp 
     263      slf_i(jp_slp  ) = sn_slp     ;   slf_i(jp_cc   ) = sn_cc 
    263264      slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm 
    264265      slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj 
     
    289290         ! 
    290291         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default) 
    291             IF(     jfpr == jp_slp  ) THEN 
     292            IF(     jfpr == jp_slp ) THEN 
    292293               sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa 
    293294            ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 
     
    295296            ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN 
    296297               DEALLOCATE( sf(jfpr)%fnow )              ! deallocate as not used in this case 
     298            ELSEIF( jfpr == jp_cc  ) THEN 
     299               sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf 
    297300            ELSE 
    298301               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr 
     
    303306            ! 
    304307            IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   & 
    305                &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
    306                &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 
     308         &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
     309         &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 
    307310         ENDIF 
    308311      END DO 
     
    559562      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 
    560563 
     564      ! --- cloud cover --- ! 
     565      cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1) 
     566 
    561567      ! ----------------------------------------------------------------------------- ! 
    562568      !      0   Wind components and module at T-point relative to the moving ocean   ! 
     
    10191025      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    10201026      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      - 
    1021       REAL(wp) ::   zfr1, zfr2               ! local variables 
    10221027      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
    10231028      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
     
    10281033      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB 
    10291034      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2 
     1035      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri 
    10301036      !!--------------------------------------------------------------------- 
    10311037      ! 
     
    11121118      ! --- evaporation minus precipitation --- ! 
    11131119      zsnw(:,:) = 0._wp 
    1114       CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
     1120      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
    11151121      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    11161122      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     
    11391145      END DO 
    11401146 
    1141       ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! 
    1142       zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm 
    1143       zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    1144       ! 
    1145       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
    1146          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    1147       ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    1148          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    1149       ELSEWHERE                                                         ! zero when hs>0 
    1150          qtr_ice_top(:,:,:) = 0._wp 
    1151       END WHERE 
    1152       ! 
    1153  
     1147      ! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- ! 
     1148      IF( nn_qtrice == 0 ) THEN 
     1149         ! formulation derived from Grenfell and Maykut (1977), where transmission rate 
     1150         !    1) depends on cloudiness 
     1151         !    2) is 0 when there is any snow 
     1152         !    3) tends to 1 for thin ice 
     1153         ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm 
     1154         DO jl = 1, jpl 
     1155            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm   
     1156               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 
     1157            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm 
     1158               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 
     1159            ELSEWHERE                                                         ! zero when hs>0 
     1160               qtr_ice_top(:,:,jl) = 0._wp  
     1161            END WHERE 
     1162         ENDDO 
     1163      ELSEIF( nn_qtrice == 1 ) THEN 
     1164         ! formulation is derived from the thesis of M. Lebrun (2019). 
     1165         !    It represents the best fit using several sets of observations 
     1166         !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 
     1167         qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:) 
     1168      ENDIF 
     1169      ! 
    11541170      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 
    11551171         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
  • NEMO/trunk/src/OCE/SBC/sbccpl.F90

    r13461 r13472  
    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    
     
    191205 
    192206   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
     207#if defined key_si3 || defined key_cice 
     208   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time 
     209#endif 
    193210 
    194211   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]  
     
    211228      !!             ***  FUNCTION sbc_cpl_alloc  *** 
    212229      !!---------------------------------------------------------------------- 
    213       INTEGER :: ierr(4) 
     230      INTEGER :: ierr(5) 
    214231      !!---------------------------------------------------------------------- 
    215232      ierr(:) = 0 
     
    221238#endif 
    222239      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
    223       ! 
    224       IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) )  
     240#if defined key_si3 || defined key_cice 
     241      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(4) ) 
     242#endif 
     243      ! 
     244      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(5) ) 
    225245 
    226246      sbc_cpl_alloc = MAXVAL( ierr ) 
     
    249269      REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos 
    250270      !! 
    251       NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   &  
     271      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             & 
     272         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  &  
    252273         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  &  
    253          &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   &  
    254          &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   &  
     274         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  &  
     275         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  &  
    255276         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  & 
    256          &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   & 
    257          &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   & 
    258          &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   & 
     277         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  & 
     278         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                & 
     279         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 & 
    259280         &                  sn_rcv_ts_ice 
    260  
    261281      !!--------------------------------------------------------------------- 
    262282      ! 
     
    278298      ENDIF 
    279299      IF( lwp .AND. ln_cpl ) THEN                        ! control print 
     300         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
     301         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     302         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux 
     303         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl 
    280304         WRITE(numout,*)'  received fields (mutiple ice categogies)' 
    281305         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')' 
     
    326350         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor  
    327351         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd  
    328          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    329          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
    330          WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl 
    331352      ENDIF 
    332353 
     
    367388      IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  & 
    368389           .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled 
    369  
     390      ! 
    370391      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 
    371392       
     
    822843      END SELECT 
    823844 
     845      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office) 
     846#if defined key_si3 || defined key_cice 
     847       a_i_last_couple(:,:,:) = 0._wp 
     848#endif 
    824849      !                                                      ! ------------------------- !  
    825850      !                                                      !      Ice Meltponds        !  
     
    11101135      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient 
    11111136      REAL(wp) ::   zzx, zzy               ! temporary variables 
    1112       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr 
     1137      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra 
    11131138      !!---------------------------------------------------------------------- 
    11141139      ! 
     
    12241249         ENDIF 
    12251250      ENDIF 
    1226  
     1251!!$      !                                                      ! ========================= ! 
     1252!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      ! 
     1253!!$      !                                                      ! ========================= ! 
     1254!!$      cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1) 
     1255!!$      END SELECT 
     1256!!$ 
     1257      zcloud_fra(:,:) = pp_cldf   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm. 
     1258      IF( ln_mixcpl ) THEN 
     1259         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:) 
     1260      ELSE 
     1261         cloud_fra(:,:) = zcloud_fra(:,:) 
     1262      ENDIF 
     1263      !                                                      ! ========================= ! 
    12271264      ! u(v)tau and taum will be modified by ice model 
    12281265      ! -> need to be reset before each call of the ice/fsbc       
     
    16231660      ! 
    16241661      INTEGER  ::   ji, jj, jl   ! dummy loop index 
    1625       REAL(wp) ::   ztri         ! local scalar 
    16261662      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
    16271663      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice 
    16281664      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
     1665      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total 
    16291666      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu 
     1667      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri 
    16301668      !!---------------------------------------------------------------------- 
    16311669      ! 
     
    16471685         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    16481686         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
    1649          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:) 
    16501687      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp 
    16511688         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1) 
     
    16591696 
    16601697#if defined key_si3 
     1698 
     1699      ! --- evaporation over ice (kg/m2/s) --- ! 
     1700      IF (ln_scale_ice_flux) THEN ! typically met-office requirements 
     1701         IF (sn_rcv_emp%clcat == 'yes') THEN 
     1702            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     1703            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp 
     1704            END WHERE 
     1705            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:) 
     1706            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp 
     1707            END WHERE 
     1708         ELSE 
     1709            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:) 
     1710            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp 
     1711            END WHERE 
     1712            zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1713            DO jl = 2, jpl 
     1714               zevap_ice(:,:,jl) = zevap_ice(:,:,1) 
     1715            ENDDO 
     1716         ENDIF 
     1717      ELSE 
     1718         IF (sn_rcv_emp%clcat == 'yes') THEN 
     1719            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl) 
     1720            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:) 
     1721            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp 
     1722            END WHERE 
     1723         ELSE 
     1724            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) 
     1725            zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1726            DO jl = 2, jpl 
     1727               zevap_ice(:,:,jl) = zevap_ice(:,:,1) 
     1728            ENDDO 
     1729         ENDIF 
     1730      ENDIF 
     1731 
     1732      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN 
     1733         ! For conservative case zemp_ice has not been defined yet. Do it now. 
     1734         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:) 
     1735      ENDIF 
     1736 
    16611737      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing) 
    1662       zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw ) 
     1738      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw ) 
    16631739       
    16641740      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- ! 
     
    16671743 
    16681744      ! --- 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 
     1745      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) 
    16761746 
    16771747      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
     
    17511821!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff 
    17521822!!      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       IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
    1756       IF( iom_use('precip') )       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) 
     1823      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving 
     1824      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs 
     1825      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
     1826      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
     1827      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
     1828      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
     1829      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
     1830      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average) 
     1831      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) 
     1832      IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
     1833         &                                                         - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 
    17641834      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf 
    17651835      ! 
     
    17691839      CASE( 'oce only' )         ! the required field is directly provided 
    17701840         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) 
     1841         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero 
     1842         ! here so the only flux is the ocean only one. 
     1843         zqns_ice(:,:,:) = 0._wp  
    17711844      CASE( 'conservative' )     ! the required fields are directly provided 
    17721845         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
     
    17981871               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    & 
    17991872                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
    1800                   &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1873                  &                                             + pist(:,:,jl) * picefr(:,:) ) ) 
    18011874            END DO 
    18021875         ELSE 
     
    18041877               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    & 
    18051878                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
    1806                   &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1879                  &                                             + pist(:,:,jl) * picefr(:,:) ) ) 
    18071880            END DO 
    18081881         ENDIF 
     
    19101983      CASE( 'oce only' ) 
    19111984         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) ) 
     1985         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero 
     1986         ! here so the only flux is the ocean only one. 
     1987         zqsr_ice(:,:,:) = 0._wp 
    19121988      CASE( 'conservative' ) 
    19131989         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
     
    19952071            ENDDO 
    19962072         ENDIF 
     2073      CASE( 'none' )  
     2074         zdqns_ice(:,:,:) = 0._wp 
    19972075      END SELECT 
    19982076       
     
    20102088      !                                                      ! ========================= ! 
    20112089      CASE ('coupled') 
    2012          IF( ln_mixcpl ) THEN 
    2013             DO jl=1,jpl 
    2014                qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:) 
    2015                qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:) 
    2016             ENDDO 
     2090         IF (ln_scale_ice_flux) THEN 
     2091            WHERE( a_i(:,:,:) > 1.e-10_wp ) 
     2092               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2093               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2094            ELSEWHERE 
     2095               qml_ice(:,:,:) = 0.0_wp 
     2096               qcn_ice(:,:,:) = 0.0_wp 
     2097            END WHERE 
    20172098         ELSE 
    20182099            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
     
    20252106      IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==! 
    20262107         ! 
    2027          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    2028          ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77) 
    2029          ! 
    2030          WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
    2031             zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    2032          ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
    2033             zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri 
    2034          ELSEWHERE                                                         ! zero when hs>0 
    2035             zqtr_ice_top(:,:,:) = 0._wp 
    2036          END WHERE 
     2108         IF( nn_qtrice == 0 ) THEN 
     2109            ! formulation derived from Grenfell and Maykut (1977), where transmission rate 
     2110            !    1) depends on cloudiness 
     2111            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     2112            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm. 
     2113            !    2) is 0 when there is any snow 
     2114            !    3) tends to 1 for thin ice 
     2115            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm 
     2116            DO jl = 1, jpl 
     2117               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     2118                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 
     2119               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
     2120                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:) 
     2121               ELSEWHERE                                                           ! zero when hs>0 
     2122                  zqtr_ice_top(:,:,jl) = 0._wp  
     2123               END WHERE 
     2124            ENDDO 
     2125         ELSEIF( nn_qtrice == 1 ) THEN 
     2126            ! formulation is derived from the thesis of M. Lebrun (2019). 
     2127            !    It represents the best fit using several sets of observations 
     2128            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 
     2129            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:) 
     2130         ENDIF 
    20372131         !      
    20382132      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==! 
    20392133         ! 
    2040          !                    ! ===> here we must receive the qtr_ice_top array from the coupler 
    2041          !                           for now just assume zero (fully opaque ice) 
     2134         !          ! ===> here we must receive the qtr_ice_top array from the coupler 
     2135         !                 for now just assume zero (fully opaque ice) 
    20422136         zqtr_ice_top(:,:,:) = 0._wp 
    20432137         ! 
     
    20962190      ! 
    20972191      isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges 
     2192      info = OASIS_idle 
    20982193 
    20992194      zfr_l(:,:) = 1.- fr_i(:,:) 
     
    22342329      ENDIF 
    22352330 
     2331#if defined key_si3 || defined key_cice 
     2332      ! If this coupling was successful then save ice fraction for use between coupling points.  
     2333      ! This is needed for some calculations where the ice fraction at the last coupling point  
     2334      ! is needed.  
     2335      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. &  
     2336         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN  
     2337         IF ( sn_snd_thick%clcat == 'yes' ) THEN  
     2338           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl) 
     2339         ENDIF 
     2340      ENDIF 
     2341#endif 
     2342 
    22362343      IF( ssnd(jps_fice1)%laction ) THEN 
    22372344         SELECT CASE( sn_snd_thick1%clcat ) 
     
    22972404            SELECT CASE( sn_snd_mpnd%clcat )   
    22982405            CASE( 'yes' )   
    2299                ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) 
     2406               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl) 
    23002407               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)   
    23012408            CASE( 'no' )   
     
    23032410               ztmp4(:,:,:) = 0.0   
    23042411               DO jl=1,jpl   
    2305                  ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)   
    2306                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)  
     2412                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl) 
     2413                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
    23072414               ENDDO   
    23082415            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )   
  • NEMO/trunk/src/OCE/SBC/sbcmod.F90

    r13286 r13472  
    563563      ENDIF 
    564564      ! 
    565       CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice) 
    566       CALL iom_put( "vtau", vtau )   ! j-wind stress 
    567       ! 
    568565      IF(sn_cfctl%l_prtctl) THEN     ! print mean trends (used for debugging) 
    569566         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask ) 
  • NEMO/trunk/src/OCE/ZDF/zdfdrg.F90

    r13461 r13472  
    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                                           [ - ] 
     
    226227      INTEGER   ::   ios, ioptio   ! local integers 
    227228      !! 
    228       NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp 
     229      NAMELIST/namdrg/ ln_drg_OFF, ln_lin, ln_non_lin, ln_loglayer, ln_drgimp, ln_drgice_imp 
    229230      !!---------------------------------------------------------------------- 
    230231      ! 
     
    247248         WRITE(numout,*) '      logarithmic drag: Cd = vkarmn/log(z/z0)   ln_loglayer = ', ln_loglayer 
    248249         WRITE(numout,*) '      implicit friction                         ln_drgimp   = ', ln_drgimp 
     250         WRITE(numout,*) '      implicit ice-ocean drag                   ln_drgice_imp  =', ln_drgice_imp 
    249251      ENDIF 
    250252      ! 
     
    257259      IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' ) 
    258260      ! 
     261      IF ( ln_drgice_imp.AND.(.NOT.ln_drgimp) ) &  
     262         &                CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires ln_drgimp=T' ) 
     263      ! 
     264      IF ( ln_drgice_imp.AND.( nn_ice /=2 ) ) & 
     265         &  CALL ctl_stop( 'zdf_drg_init: ln_drgice_imp=T requires si3' ) 
    259266      ! 
    260267      !                     !==  BOTTOM drag setting  ==!   (applied at seafloor) 
     
    267274      !                     !==  TOP drag setting  ==!   (applied at the top of ocean cavities) 
    268275      ! 
    269       IF( ln_isfcav ) THEN              ! Ocean cavities: top friction setting 
    270          ALLOCATE( rCd0_top(jpi,jpj), rCdU_top(jpi,jpj) ) 
     276      IF( ln_isfcav.OR.ln_drgice_imp ) THEN              ! Ocean cavities: top friction setting 
     277         ALLOCATE( rCdU_top(jpi,jpj) ) 
     278      ENDIF 
     279      ! 
     280      IF( ln_isfcav ) THEN 
     281         ALLOCATE( rCd0_top(jpi,jpj)) 
    271282         CALL drg_init( 'TOP   '   , mikt       ,                                         &   ! <== in 
    272283            &           r_Cdmin_top, r_Cdmax_top, r_z0_top, r_ke0_top, rCd0_top, rCdU_top )   ! ==> out 
  • NEMO/trunk/src/OCE/ZDF/zdfgls.F90

    r13461 r13472  
    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 
     
    153155      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves  
    154156      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves) 
     157      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice 
    155158      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before 
    156159      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before 
     
    168171      ustar2_bot (:,:) = 0._wp 
    169172 
     173      SELECT CASE ( nn_z0_ice ) 
     174      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     175      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     176      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     177      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     178      END SELECT 
     179       
    170180      ! Compute surface, top and bottom friction at T-points 
    171181      DO_2D( 0, 0, 0, 0 ) 
     
    207217      END SELECT 
    208218      ! 
     219      ! adapt roughness where there is sea ice 
     220      zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro 
     221      ! 
    209222      DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 
    210223         eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 
     
    291304      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)  
    292305      ! First level 
    293       en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  ) 
     306      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  ) 
    294307      zd_lw(:,:,1) = en(:,:,1) 
    295308      zd_up(:,:,1) = 0._wp 
     
    297310      !  
    298311      ! One level below 
    299       en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm))  & 
    300          &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   ) 
     312      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 
     313         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   ) 
    301314      zd_lw(:,:,2) = 0._wp  
    302315      zd_up(:,:,2) = 0._wp 
     
    307320      ! 
    308321      ! Dirichlet conditions at k=1 
    309       en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  ) 
     322      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  ) 
    310323      zd_lw(:,:,1) = en(:,:,1) 
    311324      zd_up(:,:,1) = 0._wp 
     
    317330      zd_lw(:,:,2) = 0._wp 
    318331      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 
    319       zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
     332      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 
    320333          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf) 
    321334!!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 
     
    530543         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 
    531544         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 
    532          zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     545         zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) & 
     546            &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    533547         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 
    534548            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 
     
    753767      REAL(wp)::   zcr   ! local scalar 
    754768      !! 
    755       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    756          &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    757          &            rn_crban, rn_charn, rn_frac_hs,        & 
    758          &            nn_bc_surf, nn_bc_bot, nn_z0_met,     & 
     769      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim,       & 
     770         &            rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri,   & 
     771         &            rn_crban, rn_charn, rn_frac_hs,              & 
     772         &            nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, & 
    759773         &            nn_stab_func, nn_clos 
    760774      !!---------------------------------------------------------- 
     
    782796         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    783797         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     798         WRITE(numout,*) '      surface wave breaking under ice               nn_z0_ice      = ', nn_z0_ice 
     799         SELECT CASE( nn_z0_ice ) 
     800         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on surface wave breaking' 
     801         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     802         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-fr_i(:,:)' 
     803         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     804         CASE DEFAULT 
     805            CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3') 
     806         END SELECT 
    784807         WRITE(numout,*) '      Wave height frac. (used if nn_z0_met=2)       rn_frac_hs     = ', rn_frac_hs 
    785808         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func 
    786809         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    787810         WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro 
     811         WRITE(numout,*) '      Ice-ocean roughness (used if nn_z0_ice/=0)    rn_hsri        = ', rn_hsri 
    788812         WRITE(numout,*) 
    789813         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:' 
  • NEMO/trunk/src/OCE/ZDF/zdfphy.F90

    r13226 r13472  
    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 
     
    253254      ENDIF 
    254255      ! 
     256#if defined key_si3 
     257      IF ( ln_drgice_imp) THEN 
     258         IF ( ln_isfcav ) THEN 
     259            rCdU_top(:,:) = rCdU_top(:,:) + ssmask(:,:) * tmask(:,:,1) * rCdU_ice(:,:) 
     260         ELSE 
     261            rCdU_top(:,:) = rCdU_ice(:,:) 
     262         ENDIF 
     263      ENDIF 
     264#endif 
     265      !  
    255266      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k) 
    256267      ! 
  • NEMO/trunk/src/OCE/ZDF/zdftke.F90

    r13461 r13472  
    5252#endif 
    5353   ! 
     54#if defined key_si3 
     55   USE ice, ONLY: hm_i, h_i 
     56#endif 
     57#if defined key_cice 
     58   USE sbc_ice, ONLY: h_i 
     59#endif 
    5460   USE in_out_manager ! I/O manager 
    5561   USE iom            ! I/O manager library 
     
    6874   !                      !!** Namelist  namzdf_tke  ** 
    6975   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not 
     76   INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 
     77   REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice 
    7078   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3) 
    7179   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m] 
    72    INTEGER  ::      nn_mxlice ! type of scaling under sea-ice 
    73    REAL(wp) ::      rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1) 
    7480   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1) 
    7581   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
     
    8288   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
    8389   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
    84    REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4    
    8590   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8691   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     92   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)    
    8793 
    8894   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
     
    199205      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    200206      ! 
    201       INTEGER ::   ji, jj, jk              ! dummy loop arguments 
     207      INTEGER ::   ji, jj, jk                  ! dummy loop arguments 
    202208      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars 
    203209      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3 
    204210      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient 
    205       REAL(wp) ::   zbbrau, zri                ! local scalars 
    206       REAL(wp) ::   zfact1, zfact2, zfact3     !   -         - 
    207       REAL(wp) ::   ztx2  , zty2  , zcof       !   -         - 
    208       REAL(wp) ::   ztau  , zdif               !   -         - 
    209       REAL(wp) ::   zus   , zwlc  , zind       !   -         - 
    210       REAL(wp) ::   zzd_up, zzd_lw             !   -         - 
     211      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars 
     212      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      - 
     213      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      - 
     214      REAL(wp) ::   ztau  , zdif               !   -      - 
     215      REAL(wp) ::   zus   , zwlc  , zind       !   -      - 
     216      REAL(wp) ::   zzd_up, zzd_lw             !   -      - 
    211217      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
    212       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i 
     218      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3 
    213219      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    214220      !!-------------------------------------------------------------------- 
    215221      ! 
    216       zbbrau = rn_ebb / rho0       ! Local constant initialisation 
    217       zfact1 = -.5_wp * rn_Dt  
    218       zfact2 = 1.5_wp * rn_Dt * rn_ediss 
    219       zfact3 = 0.5_wp       * rn_ediss 
     222      zbbrau  = rn_ebb / rho0       ! Local constant initialisation 
     223      zbbirau = 3.75_wp / rho0 
     224      zfact1  = -.5_wp * rn_Dt  
     225      zfact2  = 1.5_wp * rn_Dt * rn_ediss 
     226      zfact3  = 0.5_wp         * rn_ediss 
     227      ! 
     228      ! ice fraction considered for attenuation of langmuir & wave breaking 
     229      SELECT CASE ( nn_eice ) 
     230      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     231      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     232      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     233      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     234      END SELECT 
    220235      ! 
    221236      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    222237      !                     !  Surface/top/bottom boundary condition on tke 
    223238      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    224       !  
     239      ! 
    225240      DO_2D( 0, 0, 0, 0 ) 
     241!! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 
     242!!       one way around would be to increase zbbirau  
     243!!          en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 
     244!!             &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 
    226245         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    227246      END_2D 
     
    273292         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    274293         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    275          DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 
    276             zus  = zcof * taum(ji,jj) 
     294         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! Last w-level at which zpelc>=0.5*us*us  
     295            zus = zcof * taum(ji,jj)          !      with us=0.016*wind(starting from jpk-1) 
    277296            IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk 
    278297         END_3D 
     
    284303         DO_2D( 0, 0, 0, 0 ) 
    285304            zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    286             zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    287             IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     305            zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    288306         END_2D 
    289          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    290             IF ( zfr_i(ji,jj) /= 0. ) THEN                
    291                ! vertical velocity due to LC    
     307         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en 
     308            IF ( zus3(ji,jj) /= 0._wp ) THEN                
    292309               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    293310                  !                                           ! vertical velocity due to LC 
    294                   zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     311                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) 
    295312                  !                                           ! TKE Langmuir circulation source term 
    296                   en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     313                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    297314               ENDIF 
    298315            ENDIF 
     
    326343         !                                   ! eddy coefficient (ensure numerical stability) 
    327344         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal 
    328             &          /    (  e3t(ji,jj,jk  ,Kmm)   & 
    329             &                * e3w(ji,jj,jk  ,Kmm)  ) 
     345            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    330346         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal 
    331             &          /    (  e3t(ji,jj,jk-1,Kmm)   & 
    332             &                * e3w(ji,jj,jk  ,Kmm)  ) 
     347            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  ) 
    333348         ! 
    334349         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     
    370385       
    371386      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    372          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     387         DO_3D( 0, 0, 0, 0, 2, jpkm1 )     ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF  
    373388            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    374                &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     389               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    375390         END_3D 
    376391      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction) 
     
    378393            jk = nmln(ji,jj) 
    379394            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    380                &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     395               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    381396         END_2D 
    382397      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability) 
     
    388403            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    389404            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   & 
    390                &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     405               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    391406         END_3D 
    392407      ENDIF 
     
    450465      zmxlm(:,:,:)  = rmxl_min     
    451466      zmxld(:,:,:)  = rmxl_min 
    452       ! 
     467      !  
    453468     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
    454469         ! 
     
    468483         CASE( 1 )                           ! scaling with constant sea-ice thickness 
    469484            DO_2D( 0, 0, 0, 0 ) 
    470                zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 
     485               zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     486                  &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1) 
    471487            END_2D 
    472488            ! 
     
    474490            DO_2D( 0, 0, 0, 0 ) 
    475491#if defined key_si3 
    476                zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1) 
     492               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     493                  &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 
    477494#elif defined key_cice 
    478495               zmaxice = MAXVAL( h_i(ji,jj,:) ) 
    479                zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     496               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     497                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    480498#endif 
    481499            END_2D 
     
    484502            DO_2D( 0, 0, 0, 0 ) 
    485503               zmaxice = MAXVAL( h_i(ji,jj,:) ) 
    486                zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     504               zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     505                  &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    487506            END_2D 
    488507            ! 
     
    610629         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             & 
    611630         &                 nn_pdl  , ln_lc    , rn_lc    ,             & 
    612          &                 nn_etau , nn_htau  , rn_efr   , rn_eice   
     631         &                 nn_etau , nn_htau  , rn_efr   , nn_eice   
    613632      !!---------------------------------------------------------------------- 
    614633      ! 
     
    642661         ENDIF          
    643662         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0 
     663         IF( ln_mxl0 ) THEN 
     664            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice 
     665            IF( nn_mxlice == 1 ) & 
     666            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice 
     667            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     668            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice' 
     669            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness' 
     670            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness' 
     671            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness' 
     672            CASE DEFAULT 
     673               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 
     674            END SELECT 
     675         ENDIF 
    644676         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc 
    645677         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc 
     
    647679         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau 
    648680         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr 
    649          WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice 
    650          WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   ' 
     681         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice 
     682         SELECT CASE( nn_eice )  
     683         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking' 
     684         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     685         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)' 
     686         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     687         CASE DEFAULT 
     688            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 
     689         END SELECT       
    651690         IF( .NOT.ln_drg_OFF ) THEN 
    652691            WRITE(numout,*) 
  • NEMO/trunk/src/OCE/module_example

    r11536 r13472  
    9393      INTEGER  ::   ji, jj, jk       ! dummy loop arguments  (DOCTOR : start with j, but not jp) 
    9494      INTEGER  ::   itoto, itata     ! temporary integers    (DOCTOR : start with i 
    95       REAL(wp) ::   zmlmin, zbbrau   ! temporary scalars     (DOCTOR : start with z) 
     95      REAL(wp) ::   zmlmin, zbbrho   ! temporary scalars     (DOCTOR : start with z) 
    9696      REAL(wp) ::   zfact1, zfact2   ! do not use continuation lines in declaration 
    9797      REAL(wp), DIMENSION(jpi,jpj) ::   zwrk_2d   ! 2D workspace 
     
    101101 
    102102      zmlmin = 1.e-8                             ! Local constant initialization 
    103       zbbrau =  .5 * ebb / rau0 
     103      zbbrho =  .5 * ebb / rho0 
    104104      zfact1 = -.5 * rdt * efave 
    105105      zfact2 = 1.5 * rdt * ediss 
Note: See TracChangeset for help on using the changeset viewer.