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 15440 – NEMO

Changeset 15440


Ignore:
Timestamp:
2021-10-23T12:18:24+02:00 (2 years ago)
Author:
cetlod
Message:

dev_PISCO : merge with trunk@15439

Location:
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO
Files:
38 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg

    r15127 r15440  
    325325   sn_sal      = 'dyna_grid_T'           ,       120.        , 'vosaline'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    326326   sn_mld      = 'dyna_grid_T'           ,       120.        , 'somixhgt'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    327    sn_emp      = 'dyna_grid_T'           ,       120.        , 'sowaflup'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
     327   sn_emp      = 'dyna_grid_T'           ,       120.        , 'sowaflcd'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    328328   sn_fmf      = 'dyna_grid_T'           ,       120.        , 'iowaflup'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    329329   sn_ice      = 'dyna_grid_T'           ,       120.        , 'soicecov'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/ice.F90

    r15349 r15440  
    196196   !                                      !   = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow) 
    197197   !                                      !   = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities) 
     198 
     199   !                                     !!** namelist (namthd) ** 
     200   LOGICAL , PUBLIC ::   ln_icedH         ! activate ice thickness change from growing/melting (T) or not (F) 
     201   LOGICAL , PUBLIC ::   ln_icedA         ! activate lateral melting param. (T) or not (F) 
     202   LOGICAL , PUBLIC ::   ln_icedO         ! activate ice growth in open-water (T) or not (F) 
     203   LOGICAL , PUBLIC ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F) 
     204   LOGICAL , PUBLIC ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean 
     205   ! 
     206   !                                     !!** namelist (namthd_do) ** 
     207   REAL(wp), PUBLIC ::   rn_hinew         ! thickness for new ice formation (m) 
     208   LOGICAL , PUBLIC ::   ln_frazil        ! use of frazil ice collection as function of wind (T) or not (F) 
     209   REAL(wp), PUBLIC ::   rn_maxfraz       ! maximum portion of frazil ice collecting at the ice bottom 
     210   REAL(wp), PUBLIC ::   rn_vfraz         ! threshold drift speed for collection of bottom frazil ice 
     211   REAL(wp), PUBLIC ::   rn_Cfraz         ! squeezing coefficient for collection of bottom frazil ice 
    198212   ! 
    199213   !                                     !!** ice-vertical diffusion namelist (namthd_zdf) ** 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icectl.F90

    r15127 r15440  
    8484      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
    8585      !! 
    86       REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
    87          &          zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 
    88          &          zdiag_eimin, zdiag_esmin, zdiag_simin 
    89       REAL(wp) ::   zvtrp, zetrp 
    90       REAL(wp) ::   zarea 
    91       !!------------------------------------------------------------------- 
    92       ! 
     86      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat 
     87      REAL(wp), DIMENSION(jpi,jpj,10)     ::   ztmp3 
     88      REAL(wp), DIMENSION(jpi,jpj,jpl,8)  ::   ztmp4 
     89      REAL(wp), DIMENSION(10)             ::   zchk3          
     90      REAL(wp), DIMENSION(8)              ::   zchk4          
     91      !!------------------------------------------------------------------- 
     92      ! 
     93      ! -- quantities -- ! 
     94      ztmp3(:,:,1) = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t        ! volume 
     95      ztmp3(:,:,2) = SUM( sv_i * rhoi, dim=3 ) * e1e2t                                             ! salt 
     96      ztmp3(:,:,3) = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ! heat 
     97      ! 
     98      ! -- fluxes -- ! 
     99      ztmp3(:,:,4) = ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd &  ! mass 
     100         &          + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t 
     101      ztmp3(:,:,5) = ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw &                                ! salt 
     102         &          + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t 
     103      ztmp3(:,:,6) = ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &                                ! heat 
     104         &          - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t 
     105      ! 
     106      ! -- global sum -- ! 
     107      zchk3(1:6) = glob_sum_vec( 'icectl', ztmp3(:,:,1:6) ) 
     108 
    93109      IF( icount == 0 ) THEN 
    94  
    95          pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 
    96          pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) 
    97          pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 
    98  
    99          ! mass flux 
    100          pdiag_fv = glob_sum( 'icectl',  & 
    101             &                         ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
    102             &                           wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t ) 
    103          ! salt flux 
    104          pdiag_fs = glob_sum( 'icectl',  & 
    105             &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 
    106             &                           sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 
    107          ! heat flux 
    108          pdiag_ft = glob_sum( 'icectl',  & 
    109             &                         (   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  & 
    110             &                           - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) 
    111  
     110         ! 
     111         pdiag_v  = zchk3(1) 
     112         pdiag_s  = zchk3(2) 
     113         pdiag_t  = zchk3(3) 
     114         pdiag_fv = zchk3(4) 
     115         pdiag_fs = zchk3(5) 
     116         pdiag_ft = zchk3(6) 
     117         ! 
    112118      ELSEIF( icount == 1 ) THEN 
    113  
    114          ! -- mass diag -- ! 
    115          zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t )      & 
    116             &            - pdiag_v ) * r1_Dt_ice                                                                          & 
    117             &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       & 
    118             &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 
    119             &                                 wfx_ice_sub + wfx_spr ) * e1e2t )                                           & 
    120             &         - pdiag_fv 
    121119         ! 
    122          ! -- salt diag -- ! 
    123          zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_Dt_ice  & 
    124             &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni +           & 
    125             &                                 sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 
    126             &         - pdiag_fs 
    127          ! 
    128          ! -- heat diag -- ! 
    129          zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t & 
    130             &         ) * r1_Dt_ice                                                                                           & 
    131             &         + glob_sum( 'icectl', (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                      & 
    132             &                                - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )                    & 
    133             &         - pdiag_ft 
    134  
    135          ! -- min/max diag -- ! 
    136          zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    137          zdiag_vimin = glob_min( 'icectl', v_i  ) 
    138          zdiag_vsmin = glob_min( 'icectl', v_s  ) 
    139          zdiag_vpmin = glob_min( 'icectl', v_ip ) 
    140          zdiag_vlmin = glob_min( 'icectl', v_il ) 
    141          zdiag_aimin = glob_min( 'icectl', a_i  ) 
    142          zdiag_simin = glob_min( 'icectl', sv_i ) 
    143          zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
    144          zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
     120         ! -- mass, salt and heat diags -- ! 
     121         zdiag_mass = ( zchk3(1) - pdiag_v ) * r1_Dt_ice + ( zchk3(4) - pdiag_fv ) 
     122         zdiag_salt = ( zchk3(2) - pdiag_s ) * r1_Dt_ice + ( zchk3(5) - pdiag_fs ) 
     123         zdiag_heat = ( zchk3(3) - pdiag_t ) * r1_Dt_ice + ( zchk3(6) - pdiag_ft ) 
     124 
     125         ! -- max concentration diag -- ! 
     126         ztmp3(:,:,7) = SUM( a_i, dim=3 ) 
     127         zchk3(7)     = glob_max( 'icectl', ztmp3(:,:,7) ) 
    145128 
    146129         ! -- advection scheme is conservative? -- ! 
    147          zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 
    148          zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
    149  
    150          ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
    151          zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     130         ztmp3(:,:,8 ) = diag_adv_mass * e1e2t  
     131         ztmp3(:,:,9 ) = diag_adv_heat * e1e2t  
     132         ztmp3(:,:,10) = SUM( a_i + epsi10, dim=3 ) * e1e2t ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
     133         zchk3(8:10)   = glob_sum_vec( 'icectl', ztmp3(:,:,8:10) ) 
     134          
     135         ! -- min diags -- ! 
     136         ztmp4(:,:,:,1) = v_i 
     137         ztmp4(:,:,:,2) = v_s 
     138         ztmp4(:,:,:,3) = v_ip 
     139         ztmp4(:,:,:,4) = v_il 
     140         ztmp4(:,:,:,5) = a_i 
     141         ztmp4(:,:,:,6) = sv_i 
     142         ztmp4(:,:,:,7) = SUM( e_i, dim=3 ) 
     143         ztmp4(:,:,:,8) = SUM( e_s, dim=3 ) 
     144         zchk4(1:8)     = glob_min_vec( 'icectl', ztmp4(:,:,:,1:8) ) 
    152145 
    153146         IF( lwp ) THEN 
    154147            ! check conservation issues 
    155             IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 
     148            IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zchk3(10) ) & 
    156149               &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rDt_ice 
    157             IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
     150            IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zchk3(10) ) & 
    158151               &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rDt_ice 
    159             IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 
     152            IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zchk3(10) ) & 
    160153               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
    161154            ! check negative values 
    162             IF( zdiag_vimin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i  < 0        = ',zdiag_vimin 
    163             IF( zdiag_vsmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_s  < 0        = ',zdiag_vsmin 
    164             IF( zdiag_vpmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_ip < 0        = ',zdiag_vpmin 
    165             IF( zdiag_vlmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_il < 0        = ',zdiag_vlmin 
    166             IF( zdiag_aimin < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i  < 0        = ',zdiag_aimin 
    167             IF( zdiag_simin < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i  < 0        = ',zdiag_simin 
    168             IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i  < 0        = ',zdiag_eimin 
    169             IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s  < 0        = ',zdiag_esmin 
     155            IF( zchk4(1) < 0. )   WRITE(numout,*)   cd_routine,' : violation v_i  < 0        = ',zchk4(1) 
     156            IF( zchk4(2) < 0. )   WRITE(numout,*)   cd_routine,' : violation v_s  < 0        = ',zchk4(2) 
     157            IF( zchk4(3) < 0. )   WRITE(numout,*)   cd_routine,' : violation v_ip < 0        = ',zchk4(3) 
     158            IF( zchk4(4) < 0. )   WRITE(numout,*)   cd_routine,' : violation v_il < 0        = ',zchk4(4) 
     159            IF( zchk4(5) < 0. )   WRITE(numout,*)   cd_routine,' : violation a_i  < 0        = ',zchk4(5) 
     160            IF( zchk4(6) < 0. )   WRITE(numout,*)   cd_routine,' : violation s_i  < 0        = ',zchk4(6) 
     161            IF( zchk4(7) < 0. )   WRITE(numout,*)   cd_routine,' : violation e_i  < 0        = ',zchk4(7) 
     162            IF( zchk4(8) < 0. )   WRITE(numout,*)   cd_routine,' : violation e_s  < 0        = ',zchk4(8) 
    170163            ! check maximum ice concentration 
    171             IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
    172                &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_aimax 
     164            IF( zchk3(7)>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
     165               &                  WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zchk3(7) 
    173166            ! check if advection scheme is conservative 
    174             IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    175                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 
    176             IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    177                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rDt_ice 
     167            IF( ABS(zchk3(8)) > zchk_m * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) & 
     168               &                  WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zchk3(8) * rDt_ice 
     169            IF( ABS(zchk3(9)) > zchk_t * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) & 
     170               &                  WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zchk3(9) * rDt_ice 
    178171         ENDIF 
    179172         ! 
    … …  
    195188      !!------------------------------------------------------------------- 
    196189      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
    197       REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat 
    198       REAL(wp) ::   zarea 
    199       !!------------------------------------------------------------------- 
    200  
    201       ! water flux 
    202       ! -- mass diag -- ! 
    203       zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr   + wfx_sub + wfx_pnd & 
    204          &                              + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 
    205  
    206       ! -- salt diag -- ! 
    207       zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 
    208  
    209       ! -- heat diag -- ! 
    210       zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     190      !! 
     191      REAL(wp), DIMENSION(jpi,jpj,4)     ::   ztmp 
     192      REAL(wp), DIMENSION(4)             ::   zchk          
     193      !!------------------------------------------------------------------- 
     194 
     195      ztmp(:,:,1) = ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ! mass diag 
     196      ztmp(:,:,2) = ( sfx + diag_sice - diag_adv_salt ) * e1e2t                                                                     ! salt 
     197      ztmp(:,:,3) = ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t                                                   ! heat 
    211198      ! equivalent to this: 
    212       !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 
    213       !!   &                                          - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 
    214       !!   &                                          ) * e1e2t ) 
    215  
    216       ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
    217       zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
    218  
     199      !! ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 
     200      !!   &                                        - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) 
     201      ztmp(:,:,4) =  SUM( a_i + epsi10, dim=3 ) * e1e2t      ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
     202 
     203      ! global sums 
     204      zchk(1:4)   = glob_sum_vec( 'icectl', ztmp(:,:,1:4) ) 
     205       
    219206      IF( lwp ) THEN 
    220          IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 
    221             &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rDt_ice 
    222          IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
    223             &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rDt_ice 
    224          IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 
    225             &                   WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
     207         IF( ABS(zchk(1)) > zchk_m * rn_icechk_glo * zchk(4) ) & 
     208            &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zchk(1) * rDt_ice 
     209         IF( ABS(zchk(2)) > zchk_s * rn_icechk_glo * zchk(4) ) & 
     210            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zchk(2) * rDt_ice 
     211         IF( ABS(zchk(3)) > zchk_t * rn_icechk_glo * zchk(4) ) & 
     212            &                   WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zchk(3) * rDt_ice 
    226213      ENDIF 
    227214      ! 
    … …  
    762749      INTEGER, INTENT(in) ::   kt   ! ice time-step index 
    763750      ! 
    764       INTEGER  ::   ji, jj 
    765       REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 
    766       REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 
     751      REAL(wp), DIMENSION(jpi,jpj,6) ::   ztmp 
     752      REAL(wp), DIMENSION(6)         ::   zchk 
    767753      !!------------------------------------------------------------------- 
    768754      ! 
    … …  
    773759      ENDIF 
    774760      ! 
    775       ! 2D budgets (must be close to 0) 
    776       IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 
    777          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    778             zdiag_mass2D(ji,jj) = wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj)   + wfx_sub(ji,jj) + wfx_pnd(ji,jj) & 
    779                &                + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj) 
    780             zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 
    781             zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 
    782          END_2D 
    783          ! 
    784          ! write outputs 
    785          CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 
    786          CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 
    787          CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 
    788       ENDIF 
    789  
    790       ! -- mass diag -- ! 
    791       zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr   + wfx_sub + wfx_pnd & 
    792          &                                  + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) * rDt_ice 
    793       zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 
    794  
    795       ! -- salt diag -- ! 
    796       zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rDt_ice * 1.e-3 
    797       zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 
    798  
    799       ! -- heat diag -- ! 
    800       zdiag_heat     = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
    801       zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
    802  
     761      ! -- 2D budgets (must be close to 0) -- ! 
     762      ztmp(:,:,1) =  wfx_ice  (:,:) + wfx_snw  (:,:) + wfx_spr  (:,:) + wfx_sub(:,:) + wfx_pnd(:,:) & 
     763         &         + diag_vice(:,:) + diag_vsnw(:,:) + diag_vpnd(:,:) - diag_adv_mass(:,:) 
     764      ztmp(:,:,2) = sfx(:,:) + diag_sice(:,:) - diag_adv_salt(:,:) 
     765      ztmp(:,:,3) = qt_oce_ai(:,:) - qt_atm_oi(:,:) + diag_heat(:,:) - diag_adv_heat(:,:) 
     766 
     767      ! write outputs 
     768      CALL iom_put( 'icedrift_mass', ztmp(:,:,1) ) 
     769      CALL iom_put( 'icedrift_salt', ztmp(:,:,2) ) 
     770      CALL iom_put( 'icedrift_heat', ztmp(:,:,3) ) 
     771 
     772      ! -- 1D budgets -- ! 
     773      ztmp(:,:,1) = ztmp(:,:,1) * e1e2t * rDt_ice         ! mass 
     774      ztmp(:,:,2) = ztmp(:,:,2) * e1e2t * rDt_ice * 1.e-3 ! salt 
     775      ztmp(:,:,3) = ztmp(:,:,3) * e1e2t                   ! heat 
     776 
     777      ztmp(:,:,4) = diag_adv_mass * e1e2t * rDt_ice 
     778      ztmp(:,:,5) = diag_adv_salt * e1e2t * rDt_ice * 1.e-3 
     779      ztmp(:,:,6) = diag_adv_heat * e1e2t 
     780 
     781      ! global sums 
     782      zchk(1:6) = glob_sum_vec( 'icectl', ztmp(:,:,1:6) ) 
     783       
    803784      !                    ! write out to file 
    804785      IF( lwp ) THEN 
    805786         ! check global drift (must be close to 0) 
    806          WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift     [kg]', zdiag_mass 
    807          WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift     [kg]', zdiag_salt 
    808          WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift     [W] ', zdiag_heat 
     787         WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift     [kg]', zchk(1) 
     788         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift     [kg]', zchk(2) 
     789         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift     [W] ', zchk(3) 
    809790         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
    810          WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'mass drift adv [kg]', zdiag_adv_mass 
    811          WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift adv [kg]', zdiag_adv_salt 
    812          WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift adv [W] ', zdiag_adv_heat 
     791         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'mass drift adv [kg]', zchk(4) 
     792         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift adv [kg]', zchk(5) 
     793         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift adv [W] ', zchk(6) 
    813794      ENDIF 
    814795      !                    ! drifts 
    815       rdiag_icemass = rdiag_icemass + zdiag_mass 
    816       rdiag_icesalt = rdiag_icesalt + zdiag_salt 
    817       rdiag_iceheat = rdiag_iceheat + zdiag_heat 
    818       rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 
    819       rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 
    820       rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 
     796      rdiag_icemass = rdiag_icemass + zchk(1) 
     797      rdiag_icesalt = rdiag_icesalt + zchk(2) 
     798      rdiag_iceheat = rdiag_iceheat + zchk(3) 
     799      rdiag_adv_icemass = rdiag_adv_icemass + zchk(4) 
     800      rdiag_adv_icesalt = rdiag_adv_icesalt + zchk(5) 
     801      rdiag_adv_iceheat = rdiag_adv_iceheat + zchk(6) 
    821802      ! 
    822803      !                    ! output drifts and close ascii file 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icedyn_rhg_evp.F90

    r15349 r15440  
    134134      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    135135      ! 
    136       REAL(wp) ::   zintb, zintn                                        ! dummy argument 
    137136      REAL(wp) ::   zfac_x, zfac_y 
    138       REAL(wp) ::   zshear, zdum1, zdum2 
    139137      ! 
    140138      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points 
    … …  
    181179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s) 
    182180      !! -- advect fields at the rheology time step for the calculation of strength 
     181      !!    it seems that convergence is worse when ll_advups=true. So it not really a good idea 
    183182      LOGICAL  ::   ll_advups = .FALSE. 
    184183      REAL(wp) ::   zdt_ups 
    … …  
    704703            ENDIF 
    705704            ! 
    706             CALL rhg_upstream( zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i 
    707             CALL rhg_upstream( zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i 
     705            CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i 
     706            CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i 
    708707            ! 
    709708            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! strength 
    … …  
    967966      REAL(wp)          ::   zresm           ! local real 
    968967      CHARACTER(len=20) ::   clname 
     968      LOGICAL           ::   ll_maxcvg 
     969      REAL(wp), DIMENSION(jpi,jpj,2) ::   zres 
     970      REAL(wp), DIMENSION(2)         ::   ztmp 
    969971      !!---------------------------------------------------------------------- 
    970  
     972      ll_maxcvg = .FALSE. 
     973      ! 
    971974      ! create file 
    972975      IF( kt == nit000 .AND. kiter == 1 ) THEN 
    … …  
    983986            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
    984987            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
    985             istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid ) 
     988            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 
    986989            istatus = NF90_ENDDEF(ncvgid) 
    987990         ENDIF 
    … …  
    9971000      ELSE 
    9981001         zresm = 0._wp 
    999          DO_2D( 0, 0, 0, 0 ) 
    1000             zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
    1001                &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 
    1002          END_2D 
    1003          CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
     1002         IF( ll_maxcvg ) THEN   ! error max over the domain 
     1003            DO_2D( 0, 0, 0, 0 ) 
     1004               zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1005                  &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 
     1006            END_2D 
     1007            CALL mpp_max( 'icedyn_rhg_evp', zresm ) 
     1008         ELSE                   ! error averaged over the domain 
     1009            DO_2D( 0, 0, 0, 0 ) 
     1010               zres(ji,jj,1) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1011                  &                 ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) 
     1012               zres(ji,jj,2) = pmsk15(ji,jj) 
     1013            END_2D 
     1014            ztmp(:) = glob_sum_vec( 'icedyn_rhg_evp', zres ) 
     1015            IF( ztmp(2) /= 0._wp )   zresm = ztmp(1) / ztmp(2) 
     1016         ENDIF 
    10041017      ENDIF 
    10051018 
    … …  
    10691082   END SUBROUTINE rhg_evp_rst 
    10701083 
    1071    SUBROUTINE rhg_upstream( pdt, pu, pv, pt ) 
     1084   SUBROUTINE rhg_upstream( jter, pdt, pu, pv, pt ) 
    10721085      !!--------------------------------------------------------------------- 
    10731086      !!                    ***  ROUTINE rhg_upstream  *** 
    … …  
    10751088      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer 
    10761089      !!---------------------------------------------------------------------- 
     1090      INTEGER                    , INTENT(in   ) ::   jter 
    10771091      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step 
    10781092      REAL(wp), DIMENSION(:,:  ) , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
    … …  
    10811095      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    10821096      REAL(wp) ::   ztra          ! local scalar 
    1083       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups   ! upstream fluxes 
     1097      LOGICAL  ::   ll_upsxy = .TRUE. 
     1098      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zpt   ! upstream fluxes and tracer guess 
    10841099      !!---------------------------------------------------------------------- 
    10851100      DO jl = 1, jpl 
    1086          DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    1087             zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + & 
    1088                &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj  ,jl) 
    1089             zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + & 
    1090                &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj+1,jl) 
    1091          END_2D 
     1101         IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **! 
     1102            ! 
     1103            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
     1104               zfu_ups(ji,jj) = MAX(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji+1,jj,jl) 
     1105               zfv_ups(ji,jj) = MAX(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj+1,jl) 
     1106            END_2D 
     1107            ! 
     1108         ELSE                              !** alternate directions **! 
     1109            ! 
     1110            IF( MOD(jter,2) == 1 ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     1111               ! 
     1112               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )       !-- flux in x-direction 
     1113                  zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj,jl) + & 
     1114                     &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     1115               END_2D 
     1116               ! 
     1117               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls )     !-- first guess of tracer from u-flux 
     1118                  ztra       = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) 
     1119                  zpt(ji,jj) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1120               END_2D 
     1121               ! 
     1122               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )   !-- flux in y-direction 
     1123                  zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj  ) + & 
     1124                     &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     1125               END_2D 
     1126               ! 
     1127            ELSE                          !==  even ice time step:  adv_y then adv_x  ==! 
     1128               ! 
     1129               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )       !-- flux in y-direction 
     1130                  zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj  ,jl) + & 
     1131                     &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     1132               END_2D 
     1133               ! 
     1134               DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 )     !-- first guess of tracer from v-flux 
     1135                  ztra       = - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
     1136                  zpt(ji,jj) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1137               END_2D 
     1138               ! 
     1139               DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )   !-- flux in x-direction 
     1140                  zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji  ,jj) + & 
     1141                     &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     1142               END_2D 
     1143               ! 
     1144            ENDIF 
     1145            ! 
     1146         ENDIF 
     1147         ! 
    10921148         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    1093             ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
    1094             ! 
    1095             pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1149            ztra         = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
     1150            pt(ji,jj,jl) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    10961151         END_2D 
    10971152      END DO 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icesbc.F90

    r15127 r15440  
    109109      !!                dqns_ice                                 = non solar  heat sensistivity                  [W/m2] 
    110110      !!                qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2] 
     111      !!            + these fields 
     112      !!                qsb_ice_bot                              = sensible heat at the ice bottom               [W/m2] 
     113      !!                fhld, qlead                              = heat budget in the leads                      [W/m2] 
    111114      !!            + some fields that are not used outside this module: 
    112115      !!                qla_ice                                  = latent heat flux over ice                     [W/m2] 
    … …  
    117120      INTEGER, INTENT(in) ::   kt     ! ocean time step 
    118121      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled) 
    119       ! 
    120       INTEGER  ::   ji, jj, jl      ! dummy loop index 
    121       REAL(wp) ::   zmiss_val       ! missing value retrieved from xios 
    122       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zalb, zmsk00      ! 2D workspace 
    123122      !!-------------------------------------------------------------------- 
    124123      ! 
    … …  
    130129         WRITE(numout,*)'~~~~~~~~~~~~~~~' 
    131130      ENDIF 
    132  
    133       ! get missing value from xml 
    134       CALL iom_miss_val( "icetemp", zmiss_val ) 
    135  
    136       ! --- ice albedo --- ! 
     131      !                     !== ice albedo ==! 
    137132      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice ) 
    138  
    139133      ! 
    140134      SELECT CASE( ksbc )   !== fluxes over sea ice ==! 
    … …  
    142136      CASE( jp_usr )              !--- user defined formulation 
    143137                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i ) 
    144       CASE( jp_blk, jp_abl )  !--- bulk formulation & ABL formulation 
     138      CASE( jp_blk, jp_abl )      !--- bulk formulation & ABL formulation 
    145139                                  CALL blk_ice_2    ( t_su, h_s, h_i, alb_ice, & 
    146140            &                                         theta_air_zt(:,:), q_air_zt(:,:),    &   ! #LB: known from "sbc_oce" module... 
    … …  
    156150         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) 
    157151      END SELECT 
    158  
    159       !--- output ice albedo and surface albedo ---! 
    160       IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN 
    161  
    162          ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) ) 
    163  
    164          WHERE( at_i_b < 1.e-03 ) 
    165             zmsk00(:,:) = 0._wp 
    166             zalb  (:,:) = rn_alb_oce 
    167          ELSEWHERE 
    168             zmsk00(:,:) = 1._wp 
    169             zalb  (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
    170          END WHERE 
    171          ! ice albedo 
    172          CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) 
    173          ! ice+ocean albedo 
    174          zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 
    175          CALL iom_put( 'albedo' , zalb ) 
    176  
    177          DEALLOCATE( zalb, zmsk00 ) 
    178  
    179       ENDIF 
     152      !                     !== some fluxes at the ice-ocean interface and in the leads 
     153      CALL ice_flx_other 
    180154      ! 
    181155      IF( ln_timing )   CALL timing_stop('icesbc') 
    … …  
    270244 
    271245 
     246   SUBROUTINE ice_flx_other 
     247      !!----------------------------------------------------------------------- 
     248      !!                   ***  ROUTINE ice_flx_other *** 
     249      !! 
     250      !! ** Purpose :   prepare necessary fields for thermo calculations 
     251      !! 
     252      !! ** Inputs  :   u_ice, v_ice, ssu_m, ssv_m, utau, vtau 
     253      !!                frq_m, qsr_oce, qns_oce, qemp_oce, e3t_m, sst_m 
     254      !! ** Outputs :   qsb_ice_bot, fhld, qlead 
     255      !!----------------------------------------------------------------------- 
     256      INTEGER  ::   ji, jj             ! dummy loop indices 
     257      REAL(wp) ::   zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1 
     258      REAL(wp), PARAMETER ::   zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04) 
     259      REAL(wp), PARAMETER ::   zch        = 0.0057_wp   ! heat transfer coefficient 
     260      REAL(wp), DIMENSION(jpi,jpj) ::  zfric, zvel      ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
     261      !!----------------------------------------------------------------------- 
     262      ! 
     263      ! computation of friction velocity at T points 
     264      IF( ln_icedyn ) THEN 
     265         DO_2D( 0, 0, 0, 0 ) 
     266            zu_io   = u_ice(ji  ,jj  ) - ssu_m(ji  ,jj  ) 
     267            zu_iom1 = u_ice(ji-1,jj  ) - ssu_m(ji-1,jj  ) 
     268            zv_io   = v_ice(ji  ,jj  ) - ssv_m(ji  ,jj  ) 
     269            zv_iom1 = v_ice(ji  ,jj-1) - ssv_m(ji  ,jj-1) 
     270            ! 
     271            zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * tmask(ji,jj,1) 
     272            zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) + & 
     273               &                          ( v_ice(ji  ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) ) ) 
     274         END_2D 
     275      ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean 
     276         DO_2D( 0, 0, 0, 0 ) 
     277            zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  & 
     278               &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
     279               &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
     280            zvel(ji,jj) = 0._wp 
     281         END_2D 
     282      ENDIF 
     283      CALL lbc_lnk( 'icesbc', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp ) 
     284      ! 
     285      !--------------------------------------------------------------------! 
     286      ! Partial computation of forcing for the thermodynamic sea ice model 
     287      !--------------------------------------------------------------------! 
     288      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )   ! needed for qlead 
     289         rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
     290         ! 
     291         ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
     292         zqld =  tmask(ji,jj,1) * rDt_ice *  & 
     293            &    ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) +  & 
     294            &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
     295 
     296         ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 
     297         !     (mostly<0 but >0 if supercooling) 
     298         zqfr     = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 
     299         zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0 
     300         zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0 
     301 
     302         ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 
     303         !     (mostly>0 but <0 if supercooling) 
     304         zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
     305         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
     306 
     307         ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
     308         !                              the freezing point, so that we do not have SST < T_freeze 
     309         !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 
     310         !                              The following formulation is ok for both normal conditions and supercooling 
     311         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     312 
     313         ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 
     314         ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 
     315         IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 
     316            zqfr               = 0._wp 
     317            zqfr_pos           = 0._wp 
     318            qsb_ice_bot(ji,jj) = 0._wp 
     319         ENDIF 
     320         ! 
     321         ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
     322         !     qlead is the energy received from the atm. in the leads. 
     323         !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2) 
     324         !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2) 
     325         IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     326            ! upper bound for fhld: fhld should be equal to zqld 
     327            !                        but we have to make sure that this heat will not make the sst drop below the freezing point 
     328            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 
     329            !                        The following formulation is ok for both normal conditions and supercooling 
     330            fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
     331               &                                 - qsb_ice_bot(ji,jj) ) 
     332            qlead(ji,jj) = 0._wp 
     333         ELSE 
     334            fhld (ji,jj) = 0._wp 
     335            ! upper bound for qlead: qlead should be equal to zqld 
     336            !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 
     337            !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 
     338            !                        and freezing point is reached if zqfr = zqld - qsb*a/dt 
     339            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 
     340            !                        The following formulation is ok for both normal conditions and supercooling 
     341            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
     342         ENDIF 
     343         ! 
     344         ! If ice is landfast and ice concentration reaches its max 
     345         ! => stop ice formation in open water 
     346         IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp 
     347         ! 
     348         ! If the grid cell is almost fully covered by ice (no leads) 
     349         ! => stop ice formation in open water 
     350         IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp 
     351         ! 
     352         ! If ln_leadhfx is false 
     353         ! => do not use energy of the leads to melt sea-ice 
     354         IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp 
     355         ! 
     356      END_2D 
     357 
     358      ! In case we bypass open-water ice formation 
     359      IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp 
     360      ! In case we bypass growing/melting from top and bottom 
     361      IF( .NOT. ln_icedH ) THEN 
     362         qsb_ice_bot(:,:) = 0._wp 
     363         fhld       (:,:) = 0._wp 
     364      ENDIF 
     365       
     366   END SUBROUTINE ice_flx_other 
     367    
     368    
    272369   SUBROUTINE ice_sbc_init 
    273370      !!------------------------------------------------------------------- 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icethd.F90

    r15349 r15440  
    2020   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 
    2121   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    22       &                 qml_ice, qcn_ice, qtr_ice_top, utau_ice, vtau_ice 
     22      &                 qml_ice, qcn_ice, qtr_ice_top 
    2323   USE ice1D          ! sea-ice: thermodynamics variables 
    2424   USE icethd_zdf     ! sea-ice: vertical heat diffusion 
    … …  
    4848   PUBLIC   ice_thd_init    ! called by ice_init 
    4949 
    50    !!** namelist (namthd) ** 
    51    LOGICAL ::   ln_icedH         ! activate ice thickness change from growing/melting (T) or not (F) 
    52    LOGICAL ::   ln_icedA         ! activate lateral melting param. (T) or not (F) 
    53    LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F) 
    54    LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F) 
    55    LOGICAL ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean 
    56  
    5750   !! for convergence tests 
    5851   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztice_cvgerr, ztice_cvgstp 
    … …  
    9285      ! 
    9386      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    94       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 
    95       REAL(wp), PARAMETER :: zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04) 
    96       REAL(wp), PARAMETER :: zch        = 0.0057_wp   ! heat transfer coefficient 
    97       REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
    98       ! 
    99       ! for collection thickness 
    100       INTEGER  ::   iter             !   -       - 
    101       REAL(wp) ::   zvfrx, zvgx, ztaux, zf                     !   -      - 
    102       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp   !   -      - 
    103       REAL(wp), PARAMETER ::   zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used) 
    104       REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                                             ! frazil ice thickness 
    105       REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )                      ! 1/SQRT(airdensity*drag) 
    106       REAL(wp), PARAMETER ::   zgamafr = 0.03_wp 
    10787      !!------------------------------------------------------------------- 
    10888      ! controls 
    … …  
    122102         ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp 
    123103      ENDIF 
    124  
    125       !---------------------------------------------! 
    126       ! computation of friction velocity at T points 
    127       !---------------------------------------------! 
    128       IF( ln_icedyn ) THEN 
    129          zu_io(:,:) = u_ice(:,:) - ssu_m(:,:) 
    130          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
    131          DO_2D( 0, 0, 0, 0 ) 
    132             zfric(ji,jj) = rn_cio * ( 0.5_wp *  & 
    133                &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    134                &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 
    135             zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 
    136                &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 
    137          END_2D 
    138       ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean 
    139          DO_2D( 0, 0, 0, 0 ) 
    140             zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  & 
    141                &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    142                &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
    143             zvel(ji,jj) = 0._wp 
    144          END_2D 
    145       ENDIF 
    146       CALL lbc_lnk( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp ) 
    147       ! 
    148       !--------------------------------------------------------------------! 
    149       ! Partial computation of forcing for the thermodynamic sea ice model 
    150       !--------------------------------------------------------------------! 
    151       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )   ! needed for qlead 
    152          rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    153          ! 
    154          ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
    155          zqld =  tmask(ji,jj,1) * rDt_ice *  & 
    156             &    ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) +  & 
    157             &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    158  
    159          ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 
    160          !     (mostly<0 but >0 if supercooling) 
    161          zqfr     = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 
    162          zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0 
    163          zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0 
    164  
    165          ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 
    166          !     (mostly>0 but <0 if supercooling) 
    167          zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
    168          qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
    169  
    170          ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
    171          !                              the freezing point, so that we do not have SST < T_freeze 
    172          !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 
    173          !                              The following formulation is ok for both normal conditions and supercooling 
    174          qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
    175  
    176          ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 
    177          ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 
    178          IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 
    179             zqfr               = 0._wp 
    180             zqfr_pos           = 0._wp 
    181             qsb_ice_bot(ji,jj) = 0._wp 
    182          ENDIF 
    183          ! 
    184          ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
    185          !     qlead is the energy received from the atm. in the leads. 
    186          !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2) 
    187          !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2) 
    188          IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    189             ! upper bound for fhld: fhld should be equal to zqld 
    190             !                        but we have to make sure that this heat will not make the sst drop below the freezing point 
    191             !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 
    192             !                        The following formulation is ok for both normal conditions and supercooling 
    193             fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    194                &                                 - qsb_ice_bot(ji,jj) ) 
    195             qlead(ji,jj) = 0._wp 
    196          ELSE 
    197             fhld (ji,jj) = 0._wp 
    198             ! upper bound for qlead: qlead should be equal to zqld 
    199             !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 
    200             !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 
    201             !                        and freezing point is reached if zqfr = zqld - qsb*a/dt 
    202             !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 
    203             !                        The following formulation is ok for both normal conditions and supercooling 
    204             qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    205          ENDIF 
    206          ! 
    207          ! If ice is landfast and ice concentration reaches its max 
    208          ! => stop ice formation in open water 
    209          IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp 
    210          ! 
    211          ! If the grid cell is almost fully covered by ice (no leads) 
    212          ! => stop ice formation in open water 
    213          IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp 
    214          ! 
    215          ! If ln_leadhfx is false 
    216          ! => do not use energy of the leads to melt sea-ice 
    217          IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp 
    218          ! 
    219       END_2D 
    220  
    221       ! In case we bypass open-water ice formation 
    222       IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp 
    223       ! In case we bypass growing/melting from top and bottom 
    224       IF( .NOT. ln_icedH ) THEN 
    225          qsb_ice_bot(:,:) = 0._wp 
    226          fhld       (:,:) = 0._wp 
    227       ENDIF 
    228  
    229  
    230       !---------------------------------------------------------! 
    231       ! Collection thickness of ice formed in leads and polynyas 
    232       !---------------------------------------------------------!     
    233       ! ht_i_new is the thickness of new ice formed in open water 
    234       ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 
    235       ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge 
    236       ! where it forms aggregates of a specific thickness called collection thickness. 
    237       ! 
    238       fraz_frac(:,:) = 0._wp 
    239       ! 
    240       ! Default new ice thickness 
    241       WHERE( qlead(:,:) < 0._wp ) ! cooling 
    242          ht_i_new(:,:) = rn_hinew 
    243       ELSEWHERE 
    244          ht_i_new(:,:) = 0._wp 
    245       END WHERE 
    246  
    247       IF( ln_frazil ) THEN 
    248          ztwogp  = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) )  ! reduced grav 
    249          ! 
    250          DO_2D( 0, 0, 0, 0 ) 
    251             IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
    252                ! -- Wind stress -- ! 
    253                ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
    254                   &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp 
    255                ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   & 
    256                   &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp 
    257                ! Square root of wind stress 
    258                ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
    259  
    260                ! -- Frazil ice velocity -- ! 
    261                rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
    262                zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
    263                zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
    264  
    265                ! -- Pack ice velocity -- ! 
    266                zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
    267                zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
    268  
    269                ! -- Relative frazil/pack ice velocity -- ! 
    270                rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
    271                zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
    272                   &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15_wp * 0.15_wp ) * rswitch 
    273  
    274                !!clem 
    275                fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 
    276                !!clem 
    277                 
    278                ! -- new ice thickness (iterative loop) -- ! 
    279                ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1_wp )    & 
    280                   &                      / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
    281                iter = 1 
    282                DO WHILE ( iter < 20 )  
    283                   zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   & 
    284                      &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 
    285                   zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
    286  
    287                   ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 
    288                   iter = iter + 1 
    289                END DO 
    290                ! 
    291                ! bound ht_i_new (though I don't see why it should be necessary) 
    292                ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 
    293                ! 
    294             ELSE 
    295                ht_i_new(ji,jj) = 0._wp 
    296             ENDIF 
    297             ! 
    298          END_2D 
    299          !  
    300          CALL lbc_lnk( 'icethd', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  ) 
    301  
    302       ENDIF 
    303  
     104      ! 
     105      CALL ice_thd_frazil             !--- frazil ice: collection thickness (ht_i_new) & fraction of frazil (fraz_frac) 
     106      ! 
    304107      !-------------------------------------------------------------------------------------------! 
    305108      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
    … …  
    369172         ENDIF 
    370173      END_2D 
    371       CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
     174      CALL lbc_lnk( 'icethd', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
    372175      ! 
    373176      ! convergence tests 
    … …  
    447250      ! 
    448251   END SUBROUTINE ice_thd_mono 
    449  
    450252 
    451253   SUBROUTINE ice_thd_1d2d( kl, kn ) 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icethd_da.F90

    r12489 r15440  
    109109      !!--------------------------------------------------------------------- 
    110110      INTEGER  ::   ji     ! dummy loop indices 
    111       REAL(wp)            ::   zastar, zdfloe, zperi, zwlat, zda 
     111      REAL(wp)            ::   zastar, zdfloe, zperi, zwlat, zda, zda_tot 
    112112      REAL(wp), PARAMETER ::   zdmax = 300._wp 
    113113      REAL(wp), PARAMETER ::   zcs   = 0.66_wp 
    114114      REAL(wp), PARAMETER ::   zm1   = 3.e-6_wp 
    115115      REAL(wp), PARAMETER ::   zm2   = 1.36_wp 
    116       ! 
    117       REAL(wp), DIMENSION(jpij) ::   zda_tot 
    118116      !!--------------------------------------------------------------------- 
    119117      ! 
    … …  
    128126         zwlat  = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2  ! Melt speed rate [m/s] 
    129127         ! 
    130          zda_tot(ji) = MIN( zwlat * zperi * rDt_ice, at_i_1d(ji) )                 ! sea ice concentration decrease (>0) 
     128         zda_tot = MIN( zwlat * zperi * rDt_ice, at_i_1d(ji) )                     ! sea ice concentration decrease (>0) 
    131129       
    132130         ! --- Distribute reduction among ice categories and calculate associated ice-ocean fluxes --- ! 
    … …  
    134132            ! decrease of concentration for the category jl 
    135133            !    each category contributes to melting in proportion to its concentration 
    136             zda = MIN( a_i_1d(ji), zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji) ) 
     134            zda = MIN( a_i_1d(ji), zda_tot * a_i_1d(ji) / at_i_1d(ji) ) 
    137135             
    138136            ! Contribution to salt flux 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icethd_do.F90

    r15349 r15440  
    1717   USE phycst         ! physical constants 
    1818   USE sbc_oce , ONLY : sss_m 
     19   USE sbc_ice , ONLY : utau_ice, vtau_ice 
    1920   USE ice1D          ! sea-ice: thermodynamics variables 
    2021   USE ice            ! sea-ice: variables 
    … …  
    3435 
    3536   PUBLIC   ice_thd_do        ! called by ice_thd 
     37   PUBLIC   ice_thd_frazil    ! called by ice_thd 
    3638   PUBLIC   ice_thd_do_init   ! called by ice_stp 
    37  
    38    !                          !!** namelist (namthd_do) ** 
    39    REAL(wp), PUBLIC ::   rn_hinew      ! thickness for new ice formation (m) 
    40    LOGICAL , PUBLIC ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F) 
    41    REAL(wp), PUBLIC ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom 
    42    REAL(wp), PUBLIC ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice 
    43    REAL(wp), PUBLIC ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice 
    4439 
    4540   !! * Substitutions 
    … …  
    337332 
    338333 
     334   SUBROUTINE ice_thd_frazil 
     335      !!----------------------------------------------------------------------- 
     336      !!                   ***  ROUTINE ice_thd_frazil *** 
     337      !! 
     338      !! ** Purpose :   frazil ice collection thickness and fraction 
     339      !! 
     340      !! ** Inputs  :   u_ice, v_ice, utau_ice, vtau_ice 
     341      !! ** Ouputs  :   ht_i_new, fraz_frac 
     342      !!----------------------------------------------------------------------- 
     343      INTEGER  ::   ji, jj             ! dummy loop indices 
     344      INTEGER  ::   iter 
     345      REAL(wp) ::   zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp 
     346      REAL(wp), PARAMETER ::   zcai    = 1.4e-3_wp                       ! ice-air drag (clem: should be dependent on coupling/forcing used) 
     347      REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                         ! frazil ice thickness 
     348      REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )  ! 1/SQRT(airdensity*drag) 
     349      REAL(wp), PARAMETER ::   zgamafr = 0.03_wp 
     350      !!----------------------------------------------------------------------- 
     351      ! 
     352      !---------------------------------------------------------! 
     353      ! Collection thickness of ice formed in leads and polynyas 
     354      !---------------------------------------------------------!     
     355      ! ht_i_new is the thickness of new ice formed in open water 
     356      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 
     357      ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge 
     358      ! where it forms aggregates of a specific thickness called collection thickness. 
     359      ! 
     360      fraz_frac(:,:) = 0._wp 
     361      ! 
     362      ! Default new ice thickness 
     363      WHERE( qlead(:,:) < 0._wp ) ! cooling 
     364         ht_i_new(:,:) = rn_hinew 
     365      ELSEWHERE 
     366         ht_i_new(:,:) = 0._wp 
     367      END WHERE 
     368 
     369      IF( ln_frazil ) THEN 
     370         ztwogp  = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) )  ! reduced grav 
     371         ! 
     372         DO_2D( 0, 0, 0, 0 ) 
     373            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
     374               ! -- Wind stress -- ! 
     375               ztaux = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
     376               ztauy = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
     377               ! Square root of wind stress 
     378               ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
     379 
     380               ! -- Frazil ice velocity -- ! 
     381               rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
     382               zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
     383               zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
     384 
     385               ! -- Pack ice velocity -- ! 
     386               zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
     387               zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
     388 
     389               ! -- Relative frazil/pack ice velocity -- ! 
     390               rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
     391               zvrel2  = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch 
     392 
     393               ! -- fraction of frazil ice -- ! 
     394               fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 
     395                
     396               ! -- new ice thickness (iterative loop) -- ! 
     397               ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1_wp )    & 
     398                  &                      / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
     399               iter = 1 
     400               DO WHILE ( iter < 20 )  
     401                  zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   & 
     402                     &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 
     403                  zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
     404 
     405                  ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 
     406                  iter = iter + 1 
     407               END DO 
     408               ! 
     409               ! bound ht_i_new (though I don't see why it should be necessary) 
     410               ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 
     411               ! 
     412            ELSE 
     413               ht_i_new(ji,jj) = 0._wp 
     414            ENDIF 
     415            ! 
     416         END_2D 
     417         !  
     418         CALL lbc_lnk( 'icethd_frazil', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  ) 
     419 
     420      ENDIF 
     421   END SUBROUTINE ice_thd_frazil 
     422 
     423    
    339424   SUBROUTINE ice_thd_do_init 
    340425      !!----------------------------------------------------------------------- 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/iceupdate.F90

    r15127 r15440  
    9292      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    9393      REAL(wp) ::   zqsr             ! New solar flux received by the ocean 
    94       REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
     94      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z2d                  ! 2D workspace 
    9595      !!--------------------------------------------------------------------- 
    9696      IF( ln_timing )   CALL timing_start('iceupdate') 
    … …  
    239239 
    240240      IF ( iom_use( 'vfxthin' ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations 
     241         ALLOCATE( z2d(jpi,jpj) ) 
    241242         WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog 
    242243         ELSEWHERE                                     ; z2d = 0._wp 
    243244         END WHERE 
    244245         CALL iom_put( 'vfxthin', wfx_opw + z2d ) 
     246         DEALLOCATE( z2d ) 
    245247      ENDIF 
    246248 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icevar.F90

    r15127 r15440  
    343343      REAL(wp) ::   z1_dS 
    344344      REAL(wp) ::   ztmp1, ztmp2, zs0, zs 
    345       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only 
     345      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_slope_s, zalpha    ! case 2 only 
    346346      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
    347347      REAL(wp), PARAMETER :: zsi1 = 4.5_wp 
    … …  
    361361      CASE( 2 )       !  time varying salinity with linear profile  ! 
    362362         !            !---------------------------------------------! 
    363          ! 
    364          ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) 
     363         z1_dS = 1._wp / ( zsi1 - zsi0 ) 
     364         ! 
     365         ALLOCATE( z_slope_s(jpi,jpj) , zalpha(jpi,jpj) ) 
    365366         ! 
    366367         DO jl = 1, jpl 
    367             DO jk = 1, nlay_i 
    368                sz_i(:,:,jk,jl)  = s_i(:,:,jl) 
    369             END DO 
    370          END DO 
    371          !                                      ! Slope of the linear profile 
    372          WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:) 
    373          ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp 
    374          END WHERE 
    375          ! 
    376          z1_dS = 1._wp / ( zsi1 - zsi0 ) 
    377          DO jl = 1, jpl 
     368 
    378369            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    379                zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  ) 
     370               !                                      ! Slope of the linear profile 
     371               IF( h_i(ji,jj,jl) > epsi20 ) THEN 
     372                  z_slope_s(ji,jj) = 2._wp * s_i(ji,jj,jl) / h_i(ji,jj,jl) 
     373               ELSE 
     374                  z_slope_s(ji,jj) = 0._wp 
     375               ENDIF 
     376               ! 
     377               zalpha(ji,jj) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  ) 
    380378               !                             ! force a constant profile when SSS too low (Baltic Sea) 
    381                IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
     379               IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj) = 0._wp 
    382380            END_2D 
    383          END DO 
    384          ! 
    385          ! Computation of the profile 
    386          DO jl = 1, jpl 
     381            ! 
     382            ! Computation of the profile 
    387383            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i ) 
    388384               !                          ! linear profile with 0 surface value 
    389                zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 
    390                zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile 
     385               zs0 = z_slope_s(ji,jj) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 
     386               zs  = zalpha(ji,jj) * zs0 + ( 1._wp - zalpha(ji,jj) ) * s_i(ji,jj,jl)     ! weighting the profile 
    391387               sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 
    392388            END_3D 
    … …  
    448444      CASE( 2 )       !  time varying salinity with linear profile  ! 
    449445         !            !---------------------------------------------! 
     446         z1_dS = 1._wp / ( zsi1 - zsi0 ) 
    450447         ! 
    451448         ALLOCATE( z_slope_s(jpij), zalpha(jpij) ) 
    452449         ! 
    453          !                                      ! Slope of the linear profile 
    454          WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti) 
    455          ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp 
    456          END WHERE 
    457  
    458          z1_dS = 1._wp / ( zsi1 - zsi0 ) 
    459450         DO ji = 1, npti 
     451            !                                      ! Slope of the linear profile 
     452            IF( h_i_1d(ji) > epsi20 ) THEN 
     453               z_slope_s(ji) = 2._wp * s_i_1d(ji) / h_i_1d(ji) 
     454            ELSE 
     455               z_slope_s(ji) = 0._wp 
     456            ENDIF 
     457            ! 
    460458            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  ) 
    461459            !                             ! force a constant profile when SSS too low (Baltic Sea) 
    462460            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp 
     461            ! 
    463462         END DO 
    464463         ! 
    … …  
    715714      bv_i (:,:,:) = 0._wp 
    716715      DO jl = 1, jpl 
    717          DO jk = 1, nlay_i 
    718             WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 
    719                bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
    720             END WHERE 
    721          END DO 
     716         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i ) 
     717            IF( t_i(ji,jj,jk,jl) < rt0 - epsi10 ) THEN 
     718               bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rTmlt * sz_i(ji,jj,jk,jl) * r1_nlay_i / ( t_i(ji,jj,jk,jl) - rt0 ) 
     719            ENDIF 
     720         END_3D 
    722721      END DO 
    723722      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
    … …  
    782781      ! temporary 
    783782      REAL(wp) :: zintn, zintb                     ! time interpolation weights [] 
    784       REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m] 
    785783      ! 
    786784      ! compute ice load used to define the equivalent ssh in lead 
    … …  
    795793         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    796794         ! 
    797          zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rho0 
     795         ! compute equivalent ssh in lead 
     796         ice_var_sshdyn(:,:) = pssh(:,:) + ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rho0 
    798797         ! 
    799798      ELSE 
    800          zsnwiceload(:,:) = 0.0_wp 
     799         ! compute equivalent ssh in lead 
     800         ice_var_sshdyn(:,:) = pssh(:,:) 
    801801      ENDIF 
    802       ! compute equivalent ssh in lead 
    803       ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:) 
    804802      ! 
    805803   END FUNCTION ice_var_sshdyn 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icewri.F90

    r15127 r15440  
    2020   USE ice            ! sea-ice: variables 
    2121   USE icevar         ! sea-ice: operations 
     22   USE icealb , ONLY : rn_alb_oce 
    2223   ! 
    2324   USE ioipsl         ! 
    … …  
    5354      REAL(wp) ::   z2da, z2db, zrho1, zrho2 
    5455      REAL(wp) ::   zmiss_val       ! missing value retrieved from xios 
    55       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast                     ! 2D workspace 
     56      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                            ! 2D workspace 
    5657      REAL(wp), DIMENSION(jpi,jpj)     ::   zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask 
    5758      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zmsk00l, zmsksnl               ! cat masks 
     59      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zfast, zalb, zmskalb      ! 2D workspace 
    5860      ! 
    5961      ! Global ice diagnostics (SIMIP) 
    … …  
    131133      IF( iom_use('vice'    ) )   CALL iom_put( 'vice'   , v_ice    )                                                       ! ice velocity v 
    132134      ! 
    133       IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN                                                              ! module of ice velocity 
     135      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN                                                              ! module of ice velocity & fast ice 
     136         ALLOCATE( zfast(jpi,jpj) ) 
    134137         DO_2D( 0, 0, 0, 0 ) 
    135138            z2da  = u_ice(ji,jj) + u_ice(ji-1,jj) 
    … …  
    144147         END WHERE 
    145148         CALL iom_put( 'fasticepres', zfast ) 
    146       ENDIF 
    147  
     149         DEALLOCATE( zfast ) 
     150      ENDIF 
     151      ! 
     152      IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN                                                                   ! ice albedo and surface albedo 
     153         ALLOCATE( zalb(jpi,jpj), zmskalb(jpi,jpj) ) 
     154         ! ice albedo 
     155         WHERE( at_i_b < 1.e-03 ) 
     156            zmskalb(:,:) = 0._wp 
     157            zalb   (:,:) = rn_alb_oce 
     158         ELSEWHERE 
     159            zmskalb(:,:) = 1._wp 
     160            zalb   (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
     161         END WHERE 
     162         CALL iom_put( 'icealb' , zalb * zmskalb + zmiss_val * ( 1._wp - zmskalb ) ) 
     163         ! ice+ocean albedo 
     164         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 
     165         CALL iom_put( 'albedo' , zalb ) 
     166         DEALLOCATE( zalb, zmskalb ) 
     167      ENDIF 
     168      ! 
    148169      ! --- category-dependent fields --- ! 
    149170      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0% 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_oce_interp.F90

    r15127 r15440  
    15781578         DO jj = j1, j2 
    15791579            DO ji = i1, i2 
    1580                zh = 0._wp 
    1581                DO jk = 1, mbkt_parent(ji, jj)-1 
    1582                   zh = zh + e3t0_parent(ji,jj,jk) 
    1583                END DO 
    1584                e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh 
     1580               IF ( mbkt_parent(ji,jj) > 1 ) THEN 
     1581                  zh = 0._wp 
     1582                  DO jk = 1, mbkt_parent(ji, jj)-1 
     1583                     zh = zh + e3t0_parent(ji,jj,jk) 
     1584                  END DO 
     1585                  e3t0_parent(ji,jj,mbkt_parent(ji,jj)) = ht0_parent(ji, jj) - zh 
     1586               ENDIF   
    15851587            END DO 
    15861588         END DO 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_oce_sponge.F90

    r15349 r15440  
    386386      INTEGER  ::   iku, ikv 
    387387      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot 
     388      REAl(wp) :: zflag, zdmod, zdtot 
    388389      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv 
    389390      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
    … …  
    401402               DO jj=j1,j2 
    402403                  DO ji=i1,i2 
    403                      tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a) 
     404                     ! JC: masking is mandatory here: before tracer field seems  
     405                     !     to hold non zero values where tmask=0 
     406                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a) * tmask(ji,jj,jk) 
    404407                  END DO 
    405408               END DO 
    … …  
    536539                  DO ji = i1,i2-1 
    537540                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    538                      ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
     541                     zdtot =  tsbdiff(ji+1,jj,jk,jn) -  tsbdiff(ji,jj,jk,jn)  
     542                     zdmod =       ts(ji+1,jj,jk,jn,Kbb_a) - ts(ji,jj,jk,jn,Kbb_a) 
     543                     zflag = 0.5_wp + SIGN(0.5_wp, zdtot*zdmod) 
     544                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( zflag * zdtot + (1._wp - zflag) * zdmod )  
    539545                  END DO 
    540546               END DO 
    … …  
    544550                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
    545551                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     552                     zdtot =  tsbdiff(ji,jj+1,jk,jn) -  tsbdiff(ji,jj,jk,jn)  
     553                     zdmod =       ts(ji,jj+1,jk,jn,Kbb_a) - ts(ji,jj,jk,jn,Kbb_a) 
     554                     zflag = 0.5_wp + SIGN(0.5_wp, zdtot*zdmod) 
     555                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( zflag * zdtot + (1._wp - zflag) * zdmod )  
    546556                  END DO 
    547557               END DO 
    … …  
    612622            DO jj=j1,j2 
    613623               DO ji=i1,i2 
    614                   tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) 
     624                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) * umask(ji,jj,jk) 
    615625               END DO 
    616626            END DO 
    … …  
    797807            DO jj=j1,j2 
    798808               DO ji=i1,i2 
    799                   tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) 
     809                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) * vmask(ji,jj,jk) 
    800810               END DO 
    801811            END DO 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_top_sponge.F90

    r14170 r15440  
    9090               DO jj=j1,j2 
    9191                  DO ji=i1,i2 
    92                      tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a) 
     92                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a) * tmask(ji,jj,jk)  
    9393                  END DO 
    9494               END DO 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ASM/asmbkg.F90

    r12377 r15440  
    3131   USE eosbn2             ! Equation of state (eos_bn2 routine) 
    3232   USE zdfmxl             ! Mixed layer depth 
    33    USE dom_oce     , ONLY :   ndastp 
     33   USE dom_oce     , ONLY : ndastp, l_istiled 
    3434   USE in_out_manager     ! I/O manager 
    3535   USE iom                ! I/O module 
    … …  
    7474      !!----------------------------------------------------------------------- 
    7575 
    76       !                                !------------------------------------------- 
    77       IF( kt == nitbkg_r ) THEN        ! Write out background at time step nitbkg_r 
    78          !                             !-----------------------------------======== 
    79          ! 
    80          WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg ) 
    81          cl_asmbkg = TRIM( cl_asmbkg ) 
    82          INQUIRE( FILE = cl_asmbkg, EXIST = llok ) 
    83          ! 
    84          IF( .NOT. llok ) THEN 
    85             IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmbkg ) 
    86             ! 
    87             !                                      ! Define the output file         
    88             CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE. ) 
    89             ! 
    90             IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
    91                zdate = REAL( ndastp ) 
    92                IF( ln_zdftke ) THEN                   ! read turbulent kinetic energy ( en ) 
    93                   IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
    94                   CALL tke_rst( nit000, 'READ' ) 
    95                ENDIF 
    96             ELSE 
    97                zdate = REAL( ndastp ) 
    98             ENDIF 
    99             ! 
    100             !                                      ! Write the information 
    101             CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate                ) 
    102             CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , uu(:,:,:,Kmm)        ) 
    103             CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vv(:,:,:,Kmm)        ) 
    104             CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , ts(:,:,:,jp_tem,Kmm) ) 
    105             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , ts(:,:,:,jp_sal,Kmm) ) 
    106             CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , ssh(:,:,Kmm)         ) 
    107             IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en ) 
    108             ! 
    109             CALL iom_close( inum ) 
    110          ENDIF 
    111          ! 
    112       ENDIF 
    11376 
    114       !                                !------------------------------------------- 
    115       IF( kt == nitdin_r ) THEN        ! Write out background at time step nitdin_r 
    116          !                             !-----------------------------------======== 
    117          ! 
    118          WRITE(cl_asmdin, FMT='(A,".nc")' ) TRIM( c_asmdin ) 
    119          cl_asmdin = TRIM( cl_asmdin ) 
    120          INQUIRE( FILE = cl_asmdin, EXIST = llok ) 
    121          ! 
    122          IF( .NOT. llok ) THEN 
    123             IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmdin ) 
    124             ! 
    125             !                                      ! Define the output file         
    126             CALL iom_open( c_asmdin, inum, ldwrt = .TRUE. ) 
    127             ! 
    128             IF( nitdin_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
    129  
    130                zdate = REAL( ndastp ) 
    131             ELSE 
    132                zdate = REAL( ndastp ) 
    133             ENDIF 
    134             ! 
    135             !                                      ! Write the information 
    136             CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate                ) 
    137             CALL iom_rstput( kt, nitdin_r, inum, 'un'     , uu(:,:,:,Kmm)        ) 
    138             CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vv(:,:,:,Kmm)        ) 
    139             CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , ts(:,:,:,jp_tem,Kmm) ) 
    140             CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , ts(:,:,:,jp_sal,Kmm) ) 
    141             CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , ssh(:,:,Kmm)         ) 
     77      IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile 
     78          !                                !------------------------------------------- 
     79          IF( kt == nitbkg_r ) THEN        ! Write out background at time step nitbkg_r 
     80             !                             !-----------------------------------======== 
     81             ! 
     82             WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg ) 
     83             cl_asmbkg = TRIM( cl_asmbkg ) 
     84             INQUIRE( FILE = cl_asmbkg, EXIST = llok ) 
     85             ! 
     86             IF( .NOT. llok ) THEN 
     87                IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmbkg ) 
     88                ! 
     89                !                                      ! Define the output file         
     90                CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE. ) 
     91                ! 
     92                IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
     93                   zdate = REAL( ndastp ) 
     94                   IF( ln_zdftke ) THEN                   ! read turbulent kinetic energy ( en ) 
     95                      IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     96                      CALL tke_rst( nit000, 'READ' ) 
     97                   ENDIF 
     98                ELSE 
     99                   zdate = REAL( ndastp ) 
     100                ENDIF 
     101                ! 
     102                !                                      ! Write the information 
     103                CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate                ) 
     104                CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , uu(:,:,:,Kmm)        ) 
     105                CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vv(:,:,:,Kmm)        ) 
     106                CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , ts(:,:,:,jp_tem,Kmm) ) 
     107                CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , ts(:,:,:,jp_sal,Kmm) ) 
     108                CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , ssh(:,:,Kmm)         ) 
     109                IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en ) 
     110                ! 
     111                CALL iom_close( inum ) 
     112             ENDIF 
     113             ! 
     114          ENDIF 
     115     
     116          !                                !------------------------------------------- 
     117          IF( kt == nitdin_r ) THEN        ! Write out background at time step nitdin_r 
     118             !                             !-----------------------------------======== 
     119             ! 
     120             WRITE(cl_asmdin, FMT='(A,".nc")' ) TRIM( c_asmdin ) 
     121             cl_asmdin = TRIM( cl_asmdin ) 
     122             INQUIRE( FILE = cl_asmdin, EXIST = llok ) 
     123             ! 
     124             IF( .NOT. llok ) THEN 
     125                IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmdin ) 
     126                ! 
     127                !                                      ! Define the output file         
     128                CALL iom_open( c_asmdin, inum, ldwrt = .TRUE. ) 
     129                ! 
     130                IF( nitdin_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
     131     
     132                   zdate = REAL( ndastp ) 
     133                ELSE 
     134                   zdate = REAL( ndastp ) 
     135                ENDIF 
     136                ! 
     137                !                                      ! Write the information 
     138                CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate                ) 
     139                CALL iom_rstput( kt, nitdin_r, inum, 'un'     , uu(:,:,:,Kmm)        ) 
     140                CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vv(:,:,:,Kmm)        ) 
     141                CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , ts(:,:,:,jp_tem,Kmm) ) 
     142                CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , ts(:,:,:,jp_sal,Kmm) ) 
     143                CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , ssh(:,:,Kmm)         ) 
    142144#if defined key_si3 
    143             IF( nn_ice == 2 ) THEN 
    144                IF( ALLOCATED(at_i) ) THEN 
    145                   CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', at_i(:,:)   ) 
    146                ELSE 
    147                   CALL ctl_warn('asm_bkg_wri: Ice concentration not written to background ',   & 
    148                      &          'as ice variable at_i not allocated on this timestep') 
    149                ENDIF 
    150             ENDIF 
     145                IF( nn_ice == 2 ) THEN 
     146                  IF( ALLOCATED(at_i) ) THEN 
     147                      CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', at_i(:,:)   ) 
     148                   ELSE 
     149                     CALL ctl_warn('asm_bkg_wri: Ice concentration not written to background ',   & 
     150                        &          'as ice variable at_i not allocated on this timestep') 
     151                  ENDIF 
     152                ENDIF 
    151153#endif 
    152             ! 
    153             CALL iom_close( inum ) 
    154          ENDIF 
    155          ! 
    156       ENDIF 
     154                ! 
     155                CALL iom_close( inum ) 
     156             ENDIF 
     157             ! 
     158          ENDIF 
     159      ENDIF ! check for last tile 
    157160      !                     
    158161   END SUBROUTINE asm_bkg_wri 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdy_oce.F90

    r13472 r15440  
    137137   TYPE(OBC_DATA) , DIMENSION(jp_bdy), TARGET      ::   dta_bdy           !: bdy external data (local process) 
    138138!$AGRIF_END_DO_NOT_TREAT 
    139    LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lsend_bdy      !: mark needed communication for given boundary, grid and neighbour 
    140    LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lrecv_bdy      !:  when searching in any direction 
     139   LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lsend_bdyolr   !: mark needed communication for given boundary, grid and neighbour 
     140   LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lrecv_bdyolr   !:  when searching in any direction (only for orlansky) 
    141141   LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lsend_bdyint   !: mark needed communication for given boundary, grid and neighbour 
    142142   LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lrecv_bdyint   !:  when searching towards the interior of the computational domain 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdydta.F90

    r13472 r15440  
    243243         ! If full velocities in boundary data, then split it into barotropic and baroclinic component 
    244244         IF( bf_alias(jp_bdyu3d)%ltotvel ) THEN     ! if we read 3D total velocity (can be true only if u3d was read) 
    245             ! 
    246245            igrd = 2                       ! zonal velocity 
    247246            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    … …  
    258257               END DO 
    259258            END DO 
     259         ENDIF   ! ltotvel 
     260         IF( bf_alias(jp_bdyv3d)%ltotvel ) THEN     ! if we read 3D total velocity (can be true only if v3d was read) 
    260261            igrd = 3                       ! meridional velocity 
    261262            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdydyn2d.F90

    r15349 r15440  
    1818   USE bdylib          ! BDY library routines 
    1919   USE phycst          ! physical constants 
    20    USE lib_mpp, ONLY: jpfillnothing 
     20   USE lib_mpp 
    2121   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2222   USE wet_dry         ! Use wet dry to get reference ssh level 
    2323   USE in_out_manager  ! 
    24    USE lib_mpp, ONLY: ctl_stop 
    2524 
    2625   IMPLICIT NONE 
    … …  
    5049      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pssh 
    5150      !! 
    52       INTEGER  ::   ib_bdy, ir     ! BDY set index, rim index 
    53       LOGICAL  ::   llrim0         ! indicate if rim 0 is treated 
    54       LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
     51      INTEGER               ::   ib_bdy, ir     ! BDY set index, rim index 
     52      INTEGER, DIMENSION(3) ::   idir3 
     53      INTEGER, DIMENSION(6) ::   idir6 
     54      LOGICAL               ::   llrim0         ! indicate if rim 0 is treated 
     55      LOGICAL, DIMENSION(8) ::   llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
     56      !!---------------------------------------------------------------------- 
    5557       
    5658      llsend2(:) = .false.   ;   llrecv2(:) = .false. 
    … …  
    8789            SELECT CASE( cn_dyn2d(ib_bdy) ) 
    8890            CASE('flather') 
    89                llsend2(1:2) = llsend2(1:2) .OR. lsend_bdyint(ib_bdy,2,1:2,ir)   ! west/east, U points 
    90                llsend2(1)   = llsend2(1)   .OR. lsend_bdyext(ib_bdy,2,1,ir)     ! neighbour might search point towards its east bdy 
    91                llrecv2(1:2) = llrecv2(1:2) .OR. lrecv_bdyint(ib_bdy,2,1:2,ir)   ! west/east, U points 
    92                llrecv2(2)   = llrecv2(2)   .OR. lrecv_bdyext(ib_bdy,2,2,ir)     ! might search point towards bdy on the east 
    93                llsend3(3:4) = llsend3(3:4) .OR. lsend_bdyint(ib_bdy,3,3:4,ir)   ! north/south, V points 
    94                llsend3(3)   = llsend3(3)   .OR. lsend_bdyext(ib_bdy,3,3,ir)     ! neighbour might search point towards its north bdy 
    95                llrecv3(3:4) = llrecv3(3:4) .OR. lrecv_bdyint(ib_bdy,3,3:4,ir)   ! north/south, V points 
    96                llrecv3(4)   = llrecv3(4)   .OR. lrecv_bdyext(ib_bdy,3,4,ir)     ! might search point towards bdy on the north 
     91               idir6 = (/ jpwe, jpea, jpsw, jpse, jpnw, jpne /) 
     92               llsend2(idir6) = llsend2(idir6) .OR. lsend_bdyint(ib_bdy,2,idir6,ir)   ! west/east, U points 
     93               idir3 = (/ jpwe, jpsw, jpnw /) 
     94               llsend2(idir3) = llsend2(idir3) .OR. lsend_bdyext(ib_bdy,2,idir3,ir)   ! nei might search point towards its east bdy 
     95               llrecv2(idir6) = llrecv2(idir6) .OR. lrecv_bdyint(ib_bdy,2,idir6,ir)   ! west/east, U points 
     96               idir3 = (/ jpea, jpse, jpne /) 
     97               llrecv2(idir3) = llrecv2(idir3) .OR. lrecv_bdyext(ib_bdy,2,idir3,ir)   ! might search point towards bdy on the east 
     98               idir6 = (/ jpso, jpno, jpsw, jpse, jpnw, jpne /) 
     99               llsend3(idir6) = llsend3(idir6) .OR. lsend_bdyint(ib_bdy,3,idir6,ir)   ! north/south, V points 
     100               idir3 = (/ jpso, jpsw, jpse /) 
     101               llsend3(idir3) = llsend3(idir3) .OR. lsend_bdyext(ib_bdy,3,idir3,ir)   ! nei might search point towards its north bdy 
     102               llrecv3(idir6) = llrecv3(idir6) .OR. lrecv_bdyint(ib_bdy,3,idir6,ir)   ! north/south, V points 
     103               idir3 = (/ jpno, jpnw, jpne /) 
     104               llrecv3(idir3) = llrecv3(idir3) .OR. lrecv_bdyext(ib_bdy,3,idir3,ir)   ! might search point towards bdy on the north 
    97105            CASE('orlanski', 'orlanski_npo') 
    98                llsend2(:) = llsend2(:) .OR. lsend_bdy(ib_bdy,2,:,ir)   ! possibly every direction, U points 
    99                llrecv2(:) = llrecv2(:) .OR. lrecv_bdy(ib_bdy,2,:,ir)   ! possibly every direction, U points 
    100                llsend3(:) = llsend3(:) .OR. lsend_bdy(ib_bdy,3,:,ir)   ! possibly every direction, V points 
    101                llrecv3(:) = llrecv3(:) .OR. lrecv_bdy(ib_bdy,3,:,ir)   ! possibly every direction, V points 
     106               llsend2(:) = llsend2(:) .OR. lsend_bdyolr(ib_bdy,2,:,ir)   ! possibly every direction, U points 
     107               llrecv2(:) = llrecv2(:) .OR. lrecv_bdyolr(ib_bdy,2,:,ir)   ! possibly every direction, U points 
     108               llsend3(:) = llsend3(:) .OR. lsend_bdyolr(ib_bdy,3,:,ir)   ! possibly every direction, V points 
     109               llrecv3(:) = llrecv3(:) .OR. lrecv_bdyolr(ib_bdy,3,:,ir)   ! possibly every direction, V points 
    102110            END SELECT 
    103111         END DO 
    … …  
    310318      INTEGER ::   ibeg, iend      ! length of rim to be treated (rim 0 or rim 1) 
    311319      LOGICAL ::   llrim0          ! indicate if rim 0 is treated 
    312       LOGICAL, DIMENSION(4) :: llsend1, llrecv1  ! indicate how communications are to be carried out 
     320      LOGICAL, DIMENSION(8) :: llsend1, llrecv1  ! indicate how communications are to be carried out 
    313321      !!---------------------------------------------------------------------- 
    314322      llsend1(:) = .false.   ;   llrecv1(:) = .false. 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdydyn3d.F90

    r15127 r15440  
    1515   USE bdy_oce         ! ocean open boundary conditions 
    1616   USE bdylib          ! for orlanski library routines 
    17    USE lib_mpp, ONLY: jpfillnothing 
     17   USE lib_mpp 
    1818   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1919   USE in_out_manager  ! 
    20    USE lib_mpp, ONLY: ctl_stop 
    2120   Use phycst 
    2221 
    … …  
    4544      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
    4645      ! 
    47       INTEGER  ::   ib_bdy, ir     ! BDY set index, rim index 
    48       LOGICAL  ::   llrim0         ! indicate if rim 0 is treated 
    49       LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
    50  
    51       !!---------------------------------------------------------------------- 
     46      INTEGER               ::   ib_bdy, ir     ! BDY set index, rim index 
     47      INTEGER, DIMENSION(6) ::   idir6 
     48      LOGICAL               ::   llrim0         ! indicate if rim 0 is treated 
     49      LOGICAL, DIMENSION(8) ::   llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
     50      !!---------------------------------------------------------------------- 
     51       
    5252      llsend2(:) = .false.   ;   llrecv2(:) = .false. 
    5353      llsend3(:) = .false.   ;   llrecv3(:) = .false. 
    … …  
    8282            SELECT CASE( cn_dyn3d(ib_bdy) ) 
    8383            CASE('orlanski', 'orlanski_npo') 
    84                llsend2(:) = llsend2(:) .OR. lsend_bdy(ib_bdy,2,:,ir)   ! possibly every direction, U points 
    85                llrecv2(:) = llrecv2(:) .OR. lrecv_bdy(ib_bdy,2,:,ir)   ! possibly every direction, U points 
    86                llsend3(:) = llsend3(:) .OR. lsend_bdy(ib_bdy,3,:,ir)   ! possibly every direction, V points 
    87                llrecv3(:) = llrecv3(:) .OR. lrecv_bdy(ib_bdy,3,:,ir)   ! possibly every direction, V points 
     84               llsend2(:) = llsend2(:) .OR. lsend_bdyolr(ib_bdy,2,:,ir)   ! possibly every direction, U points 
     85               llrecv2(:) = llrecv2(:) .OR. lrecv_bdyolr(ib_bdy,2,:,ir)   ! possibly every direction, U points 
     86               llsend3(:) = llsend3(:) .OR. lsend_bdyolr(ib_bdy,3,:,ir)   ! possibly every direction, V points 
     87               llrecv3(:) = llrecv3(:) .OR. lrecv_bdyolr(ib_bdy,3,:,ir)   ! possibly every direction, V points 
    8888            CASE('zerograd') 
    89                llsend2(3:4) = llsend2(3:4) .OR. lsend_bdyint(ib_bdy,2,3:4,ir)   ! north/south, U points 
    90                llrecv2(3:4) = llrecv2(3:4) .OR. lrecv_bdyint(ib_bdy,2,3:4,ir)   ! north/south, U points 
    91                llsend3(1:2) = llsend3(1:2) .OR. lsend_bdyint(ib_bdy,3,1:2,ir)   ! west/east, V points 
    92                llrecv3(1:2) = llrecv3(1:2) .OR. lrecv_bdyint(ib_bdy,3,1:2,ir)   ! west/east, V points 
     89               idir6 = (/ jpso, jpno, jpsw, jpse, jpnw, jpne /) 
     90               llsend2(idir6) = llsend2(idir6) .OR. lsend_bdyint(ib_bdy,2,idir6,ir)   ! north/south, U points 
     91               llrecv2(idir6) = llrecv2(idir6) .OR. lrecv_bdyint(ib_bdy,2,idir6,ir)   ! north/south, U points 
     92               idir6 = (/ jpwe, jpea, jpsw, jpse, jpnw, jpne /) 
     93               llsend3(idir6) = llsend3(idir6) .OR. lsend_bdyint(ib_bdy,3,idir6,ir)   ! west/east, V points 
     94               llrecv3(idir6) = llrecv3(idir6) .OR. lrecv_bdyint(ib_bdy,3,idir6,ir)   ! west/east, V points 
    9395            CASE('neumann') 
    9496               llsend2(:) = llsend2(:) .OR. lsend_bdyint(ib_bdy,2,:,ir)   ! possibly every direction, U points 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdyice.F90

    r15127 r15440  
    5858      INTEGER ::   ibeg, iend                           ! length of rim to be treated (rim 0 or rim 1) 
    5959      LOGICAL ::   llrim0                               ! indicate if rim 0 is treated 
    60       LOGICAL, DIMENSION(4)  :: llsend1, llrecv1        ! indicate how communications are to be carried out 
     60      LOGICAL, DIMENSION(8)  :: llsend1, llrecv1        ! indicate how communications are to be carried out 
    6161      !!---------------------------------------------------------------------- 
    6262      ! controls 
    … …  
    327327      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points 
    328328      ! 
    329       INTEGER  ::   i_bdy, jgrd       ! dummy loop indices 
    330       INTEGER  ::   ji, jj            ! local scalar 
    331       INTEGER  ::   jbdy, ir     ! BDY set index, rim index 
    332       INTEGER  ::   ibeg, iend   ! length of rim to be treated (rim 0 or rim 1) 
    333       REAL(wp) ::   zmsk1, zmsk2, zflag 
    334       LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
     329      INTEGER               ::   i_bdy, jgrd       ! dummy loop indices 
     330      INTEGER               ::   ji, jj            ! local scalar 
     331      INTEGER               ::   jbdy, ir          ! BDY set index, rim index 
     332      INTEGER               ::   ibeg, iend        ! length of rim to be treated (rim 0 or rim 1) 
     333      INTEGER, DIMENSION(3) ::   idir3 
     334      REAL(wp)              ::   zmsk1, zmsk2, zflag 
     335      LOGICAL, DIMENSION(8) :: llsend2, llrecv2, llsend3, llrecv3  ! indicate how communications are to be carried out 
    335336      !!------------------------------------------------------------------------------ 
    336337      IF( ln_timing )   CALL timing_start('bdy_ice_dyn') 
    … …  
    430431            DO jbdy = 1, nb_bdy 
    431432               IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN 
    432                   llsend2(:) = llsend2(:) .OR. lsend_bdyint(jbdy,2,:,ir)   ! possibly every direction, U points 
    433                   llsend2(1) = llsend2(1) .OR. lsend_bdyext(jbdy,2,1,ir)   ! neighbour might search point towards its west bdy 
    434                   llrecv2(:) = llrecv2(:) .OR. lrecv_bdyint(jbdy,2,:,ir)   ! possibly every direction, U points 
    435                   llrecv2(2) = llrecv2(2) .OR. lrecv_bdyext(jbdy,2,2,ir)   ! might search point towards east bdy 
     433                  llsend2(  :  ) = llsend2(  :  ) .OR. lsend_bdyint(jbdy,2,  :  ,ir)   ! possibly every direction, U points 
     434                  idir3 = (/ jpwe, jpsw, jpnw /) 
     435                  llsend2(idir3) = llsend2(idir3) .OR. lsend_bdyext(jbdy,2,idir3,ir)   ! nei might search point towards its ea bdy 
     436                  llrecv2(  :  ) = llrecv2(  :  ) .OR. lrecv_bdyint(jbdy,2,  :  ,ir)   ! possibly every direction, U points 
     437                  idir3 = (/ jpea, jpse, jpne /) 
     438                  llrecv2(idir3) = llrecv2(idir3) .OR. lrecv_bdyext(jbdy,2,idir3,ir)   ! might search point towards east bdy 
    436439               END IF 
    437440            END DO 
    … …  
    444447            DO jbdy = 1, nb_bdy 
    445448               IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN 
    446                   llsend3(:) = llsend3(:) .OR. lsend_bdyint(jbdy,3,:,ir)   ! possibly every direction, V points 
    447                   llsend3(3) = llsend3(3) .OR. lsend_bdyext(jbdy,3,3,ir)   ! neighbour might search point towards its south bdy 
    448                   llrecv3(:) = llrecv3(:) .OR. lrecv_bdyint(jbdy,3,:,ir)   ! possibly every direction, V points 
    449                   llrecv3(4) = llrecv3(4) .OR. lrecv_bdyext(jbdy,3,4,ir)   ! might search point towards north bdy 
     449                  llsend3(  :  ) = llsend3(  :  ) .OR. lsend_bdyint(jbdy,3,  :  ,ir)   ! possibly every direction, V points 
     450                  idir3 = (/ jpso, jpsw, jpse /) 
     451                  llsend3(idir3) = llsend3(idir3) .OR. lsend_bdyext(jbdy,3,idir3,ir)   ! nei might search point towards its no bdy 
     452                  llrecv3(  :  ) = llrecv3(  :  ) .OR. lrecv_bdyint(jbdy,3,  :  ,ir)   ! possibly every direction, V points 
     453                  idir3 = (/ jpno, jpnw, jpne /) 
     454                  llrecv3(idir3) = llrecv3(idir3) .OR. lrecv_bdyext(jbdy,3,idir3,ir)   ! might search point towards north bdy 
    450455               END IF 
    451456            END DO 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdyini.F90

    r15349 r15440  
    146146      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices 
    147147      INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers 
    148       INTEGER  ::   ilen1                                  !   -       - 
     148      INTEGER  ::   ilen1                           !   -       - 
     149      INTEGER  ::   iiRst, iiRnd, iiSst, iiSnd, iiSstdiag, iiSnddiag, iiSstsono, iiSndsono 
     150      INTEGER  ::   ijRst, ijRnd, ijSst, ijSnd, ijSstdiag, ijSnddiag, ijSstsono, ijSndsono 
     151      INTEGER  ::   iiout, ijout, iioutdir, ijoutdir, icnt 
     152      INTEGER  ::   iRnei, iRdiag, iRsono 
     153      INTEGER  ::   iSnei, iSdiag, iSsono                  !   -       - 
    149154      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       - 
    150155      INTEGER  ::   jpbdta                                 !   -       - 
    … …  
    163168      REAL(wp)             , DIMENSION(jpi,jpj) ::   zfmask   ! temporary fmask array excluding coastal boundary condition (shlat) 
    164169      REAL(wp)             , DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array 
     170      REAL(wp)             , DIMENSION(jpi,jpj) ::   zzbdy 
    165171      !!---------------------------------------------------------------------- 
    166172      ! 
    … …  
    562568      ! Initialize array indicating communications in bdy 
    563569      ! ------------------------------------------------- 
    564       ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) ) 
    565       lsend_bdy(:,:,:,:) = .false. 
    566       lrecv_bdy(:,:,:,:) = .false.  
     570      ALLOCATE( lsend_bdyolr(nb_bdy,jpbgrd,8,0:1), lrecv_bdyolr(nb_bdy,jpbgrd,8,0:1) ) 
     571      lsend_bdyolr(:,:,:,:) = .false. 
     572      lrecv_bdyolr(:,:,:,:) = .false.  
    567573 
    568574      DO ib_bdy = 1, nb_bdy 
    … …  
    576582               ! 
    577583               ! check if point has to be sent     to   a neighbour 
    578                ! W neighbour and on the inner left  side 
    579                IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. mpiSnei(nn_hls,jpwe) > -1 )   lsend_bdy(ib_bdy,igrd,jpwe,ir) = .TRUE. 
    580                ! E neighbour and on the inner right side 
    581                IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. mpiSnei(nn_hls,jpea) > -1 )   lsend_bdy(ib_bdy,igrd,jpea,ir) = .TRUE. 
    582                ! S neighbour and on the inner down side 
    583                IF( ij >= Njs0 .AND. ij < Njs0 + nn_hls .AND. mpiSnei(nn_hls,jpso) > -1 )   lsend_bdy(ib_bdy,igrd,jpso,ir) = .TRUE. 
    584                ! N neighbour and on the inner up   side 
    585                IF( ij <= Nje0 .AND. ij > Nje0 - nn_hls .AND. mpiSnei(nn_hls,jpno) > -1 )   lsend_bdy(ib_bdy,igrd,jpno,ir) = .TRUE. 
     584               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij <= Nje0         ) THEN   ! we inner side 
     585                  IF( mpiSnei(nn_hls,jpwe) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE. 
     586               ENDIF 
     587               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij <= Nje0         ) THEN   ! ea inner side 
     588                  IF( mpiSnei(nn_hls,jpea) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE. 
     589               ENDIF 
     590               IF( ii >= Nis0 .AND. ii <= Nie0         .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so inner side 
     591                  IF( mpiSnei(nn_hls,jpso) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE. 
     592               ENDIF 
     593               IF( ii  < Nis0                          .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so side we-halo 
     594                  IF( mpiSnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE. 
     595               ENDIF 
     596               IF( ii  > Nie0                          .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so side ea-halo  
     597                  IF( mpiSnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE. 
     598               ENDIF 
     599               IF( ii >= Nis0 .AND. ii <= Nie0         .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no inner side 
     600                  IF( mpiSnei(nn_hls,jpno) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE. 
     601               ENDIF 
     602               IF( ii  < Nis0                          .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no side we-halo 
     603                  IF( mpiSnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE. 
     604               ENDIF 
     605               IF( ii  > Nie0                          .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no side ea-halo 
     606                  IF( mpiSnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE. 
     607               ENDIF 
     608               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! sw inner corner 
     609                  IF( mpiSnei(nn_hls,jpsw) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE. 
     610               ENDIF 
     611               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! se inner corner 
     612                  IF( mpiSnei(nn_hls,jpse) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE. 
     613               ENDIF 
     614               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! nw inner corner 
     615                  IF( mpiSnei(nn_hls,jpnw) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE. 
     616               ENDIF 
     617               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! ne inner corner 
     618                  IF( mpiSnei(nn_hls,jpne) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE. 
     619               ENDIF 
    586620               ! 
    587621               ! check if point has to be received from a neighbour 
    588                ! W neighbour and on the outter left  side 
    589                IF( ii  < Nis0                          .AND. mpiRnei(nn_hls,jpwe) > -1 )   lrecv_bdy(ib_bdy,igrd,jpwe,ir) = .TRUE. 
    590                ! E neighbour and on the outter right side 
    591                IF( ii  > Nie0                          .AND. mpiRnei(nn_hls,jpea) > -1 )   lrecv_bdy(ib_bdy,igrd,jpea,ir) = .TRUE. 
    592                ! S neighbour and on the outter down side 
    593                IF( ij  < Njs0                          .AND. mpiRnei(nn_hls,jpso) > -1 )   lrecv_bdy(ib_bdy,igrd,jpso,ir) = .TRUE. 
    594                ! N neighbour and on the outter up   side 
    595                IF( ij  > Nje0                          .AND. mpiRnei(nn_hls,jpno) > -1 )   lrecv_bdy(ib_bdy,igrd,jpno,ir) = .TRUE. 
     622               IF( ii  < Nis0                  .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN   ! we side 
     623                  IF( mpiRnei(nn_hls,jpwe) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE. 
     624               ENDIF 
     625               IF( ii  > Nie0                  .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN   ! ea side 
     626                  IF( mpiRnei(nn_hls,jpea) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE. 
     627               ENDIF 
     628               IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij  < Njs0                  ) THEN   ! so side 
     629                  IF( mpiRnei(nn_hls,jpso) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE. 
     630               ENDIF 
     631               IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij  > Nje0                  ) THEN   ! no side 
     632                  IF( mpiRnei(nn_hls,jpno) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE. 
     633               ENDIF 
     634               IF( ii  < Nis0                  .AND. ij  < Njs0                  ) THEN   ! sw corner 
     635                  IF( mpiRnei(nn_hls,jpsw) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE. 
     636                  IF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE. 
     637               ENDIF 
     638               IF( ii  > Nie0                  .AND. ij  < Njs0                  ) THEN   ! se corner 
     639                  IF( mpiRnei(nn_hls,jpse) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE. 
     640                  IF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE. 
     641               ENDIF 
     642               IF( ii  < Nis0                  .AND. ij  > Nje0                  ) THEN   ! nw corner 
     643                  IF( mpiRnei(nn_hls,jpnw) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE. 
     644                  IF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE. 
     645               ENDIF 
     646               IF( ii  > Nie0                  .AND. ij  > Nje0                  ) THEN   ! ne corner 
     647                  IF( mpiRnei(nn_hls,jpne) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE. 
     648                  IF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE. 
     649               ENDIF 
    596650               ! 
    597651            END DO 
    598          END DO  ! igrd 
    599  
     652         END DO   !   igrd 
     653          
     654         ! Comment out for debug 
     655!!$         DO ir = 0,1 
     656!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   & 
     657!!$               &                              lsend = lsend_bdyolr(ib_bdy,1,:,ir), lrecv = lrecv_bdyolr(ib_bdy,1,:,ir) ) 
     658!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr T', ir ; CALL FLUSH(numout) 
     659!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   & 
     660!!$               &                              lsend = lsend_bdyolr(ib_bdy,2,:,ir), lrecv = lrecv_bdyolr(ib_bdy,2,:,ir) ) 
     661!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr U', ir ; CALL FLUSH(numout) 
     662!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   & 
     663!!$               &                              lsend = lsend_bdyolr(ib_bdy,3,:,ir), lrecv = lrecv_bdyolr(ib_bdy,3,:,ir) )     
     664!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr V', ir ; CALL FLUSH(numout) 
     665!!$         END DO 
     666          
    600667         ! Compute rim weights for FRS scheme 
    601668         ! ---------------------------------- 
    … …  
    709776      ! 
    710777      ! Check which boundaries might need communication 
    711       ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) ) 
     778      ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,8,0:1), lrecv_bdyint(nb_bdy,jpbgrd,8,0:1) ) 
    712779      lsend_bdyint(:,:,:,:) = .false. 
    713780      lrecv_bdyint(:,:,:,:) = .false.  
    714       ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) ) 
     781      ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,8,0:1), lrecv_bdyext(nb_bdy,jpbgrd,8,0:1) ) 
    715782      lsend_bdyext(:,:,:,:) = .false. 
    716783      lrecv_bdyext(:,:,:,:) = .false. 
    717784      ! 
    718       DO igrd = 1, jpbgrd 
    719          DO ib_bdy = 1, nb_bdy 
     785      DO ib_bdy = 1, nb_bdy 
     786         DO igrd = 1, jpbgrd 
    720787            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    721788               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE 
    … …  
    731798               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )   ! free ocean neighbours 
    732799               ! 
    733                ! search neighbour in the  west/east  direction 
     800               !  take care of the 4 sides 
    734801               ! 
    735                ! Rim is on the halo and computed ocean is towards exterior of mpi domain : 
    736                !      <--    (o exterior)     -->   
    737                ! (1)  o|x         OR    (2)   x|o 
    738                !       |___                 ___| 
    739                ! ==> cannot compute the point x -> need to receive it 
    740                IF( iibi==0     .OR. ii1==0     .OR. ii2==0     .OR. ii3==0     )   lrecv_bdyint(ib_bdy,igrd,jpwe,ir) = .TRUE. 
    741                IF( iibe==0                                                     )   lrecv_bdyext(ib_bdy,igrd,jpwe,ir) = .TRUE. 
    742                IF( iibi==jpi+1 .OR. ii1==jpi+1 .OR. ii2==jpi+1 .OR. ii3==jpi+1 )   lrecv_bdyint(ib_bdy,igrd,jpea,ir) = .TRUE.   
    743                IF( iibe==jpi+1                                                 )   lrecv_bdyext(ib_bdy,igrd,jpea,ir) = .TRUE.   
    744                ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo. 
    745                ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨:  
    746                ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     : 
    747                ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....: 
    748                ! ==> the neighbour cannot compute the point x -> need to send it 
    749                IF( ii ==     2*nn_hls   .AND. mpiSnei(nn_hls,jpwe) > -1 ) THEN   ! 2*nn_hls      -> ji=jpi of western neighbour 
    750                   IF( iibi==ii+1 .OR. ii1==ii+1 .OR. ii2==ii+1 .OR. ii3==ii+1 )   lsend_bdyint(ib_bdy,igrd,jpwe,ir) = .TRUE. 
    751                   IF( iibe==ii+1                                              )   lsend_bdyext(ib_bdy,igrd,jpwe,ir) = .TRUE. 
    752                ENDIF 
    753                IF( ii == jpi-2*nn_hls+1 .AND. mpiSnei(nn_hls,jpea) > -1 ) THEN   ! jpi-2*nn_hls+1-> ji=1   of eastern neighbour 
    754                   IF( iibi==ii-1 .OR. ii1==ii-1 .OR. ii2==ii-1 .OR. ii3==ii-1 )   lsend_bdyint(ib_bdy,igrd,jpea,ir) = .TRUE. 
    755                   IF( iibe==ii-1                                              )   lsend_bdyext(ib_bdy,igrd,jpea,ir) = .TRUE. 
    756                ENDIF 
     802               DO icnt = 1, 4 
     803                  SELECT CASE( icnt ) 
     804                     !                                           ... _____ 
     805                  CASE( 1 )   ! x: rim on rcvwe/sndea-side         o|  : 
     806                     !          o: potential neighbour(s)          o|x : 
     807                     !             outside of the MPI domain     ..o|__:__ 
     808                     iRnei    = jpwe             ;   iSnei    = jpea 
     809                     iiRst    = 1                ;   ijRst    = Njs0            ! Rcv we-side starting point, excluding sw-corner 
     810                     iiRnd    = nn_hls           ;   ijRnd    = Nje0            ! Rcv we-side   ending point, excluding nw-corner 
     811                     iiSst    = Nie0-nn_hls+1    ;   ijSst    = Njs0            ! Snd ea-side starting point, excluding se-corner 
     812                     iiSnd    = Nie0             ;   ijSnd    = Nje0            ! Snd ea-side   ending point, excluding ne-corner 
     813                     iioutdir = -1               ;   ijoutdir = -999            ! outside MPI domain: westward 
     814                     !                                           ______.... 
     815                  CASE( 2 )   ! x: rim on rcvea/sndwe-side          :  |o 
     816                     !          o: potential neighbour(s)           : x|o 
     817                     !             outside of the MPI domain     ___:__|o.. 
     818                     iRnei    = jpea             ;   iSnei    = jpwe 
     819                     iiRst    = Nie0+1           ;   ijRst    =  Njs0            ! Rcv ea-side starting point, excluding se-corner 
     820                     iiRnd    = jpi              ;   ijRnd    =  Nje0            ! Rcv ea-side   ending point, excluding ne-corner 
     821                     iiSst    = Nis0             ;   ijSst    =  Njs0            ! Snd we-side starting point, excluding sw-corner 
     822                     iiSnd    = Nis0+nn_hls-1    ;   ijSnd    =  Nje0            ! Snd we-side   ending point, excluding nw-corner 
     823                     iioutdir = 1                ;   ijoutdir = -999             ! outside MPI domain: eastward 
     824                     ! 
     825                  CASE( 3 )   ! x: rim on rcvso/sndno-side       |       | 
     826                     !          o: potential neighbour(s)        |¨¨¨¨¨¨¨| 
     827                     !             outside of the MPI domain     |___x___| 
     828                     !                                           : o o o : 
     829                     !                                           :       : 
     830                     iRnei    = jpso             ;   iSnei    = jpno 
     831                     iiRst    = Nis0             ;   ijRst    = 1                ! Rcv so-side starting point, excluding sw-corner 
     832                     iiRnd    = Nie0             ;   ijRnd    = nn_hls           ! Rcv so-side   ending point, excluding se-corner 
     833                     iiSst    = Nis0             ;   ijSst    = Nje0-nn_hls+1    ! Snd no-side starting point, excluding nw-corner 
     834                     iiSnd    = Nie0             ;   ijSnd    = Nje0             ! Snd no-side   ending point, excluding ne-corner 
     835                     iioutdir = -999             ;   ijoutdir = -1               ! outside MPI domain: southward 
     836                     !                                           :       : 
     837                  CASE( 4 )   ! x: rim on rcvno/sndso-side       :_o_o_o_: 
     838                     !          o: potential neighbour(s)        |   x   | 
     839                     !             outside of the MPI domain     |       | 
     840                     !                                           |¨¨¨¨¨¨¨| 
     841                     iRnei    = jpno             ;   iSnei    = jpso 
     842                     iiRst    = Nis0             ;   ijRst    = Nje0+1           ! Rcv no-side starting point, excluding nw-corner 
     843                     iiRnd    = Nie0             ;   ijRnd    = jpj              ! Rcv no-side   ending point, excluding ne-corner 
     844                     iiSst    = Nis0             ;   ijSst    = Njs0             ! Snd so-side starting point, excluding sw-corner 
     845                     iiSnd    = Nie0             ;   ijSnd    = Njs0+nn_hls-1    ! Snd so-side   ending point, excluding se-corner 
     846                     iioutdir = -999             ;   ijoutdir = 1                ! outside MPI domain: northward 
     847                  END SELECT 
     848                  ! 
     849                  IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN   ! rim point in recv side 
     850                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain? 
     851                     ! take care of neighbourg(s) in the interior of the computational domain 
     852                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain 
     853                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it 
     854                        IF( mpiRnei(nn_hls,iRnei) > -1 )   lrecv_bdyint(ib_bdy,igrd,iRnei,ir) = .TRUE. 
     855                     ENDIF 
     856                     ! take care of neighbourg in the exterior of the computational domain 
     857                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it 
     858                        IF( mpiRnei(nn_hls,iRnei) > -1 )   lrecv_bdyext(ib_bdy,igrd,iRnei,ir) = .TRUE. 
     859                     ENDIF 
     860                  ENDIF 
     861                   
     862                  IF( ii >= iiSst .AND. ii <= iiSnd .AND. ij >= ijSst .AND. ij <= ijSnd ) THEN   ! rim point in send side 
     863                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain? 
     864                     ! take care of neighbourg(s) in the interior of the computational domain 
     865                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of nei MPI domain 
     866                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> nei cannot compute it 
     867                        IF( mpiSnei(nn_hls,iSnei) > -1 )   lsend_bdyint(ib_bdy,igrd,iSnei,ir) = .TRUE.   ! -> send to nei 
     868                     ENDIF 
     869                     ! take care of neighbourg in the exterior of the computational domain 
     870                     IF( iibe == iiout .OR. ijbe == ijout ) THEN   ! Neib outside of the nei MPI domain -> nei cannot compute it 
     871                        IF( mpiSnei(nn_hls,iSnei) > -1 )   lsend_bdyext(ib_bdy,igrd,iSnei,ir) = .TRUE.   ! -> send to nei 
     872                     ENDIF 
     873                  END IF 
     874 
     875               END DO   ! 4 sides 
    757876               ! 
    758                ! search neighbour in the north/south direction    
     877               ! specific treatment for the corners 
    759878               ! 
    760                ! Rim is on the halo and computed ocean is towards exterior of mpi domain 
    761                ! ==> cannot compute the point x -> need to receive it 
    762                !(3)   |       |         ^   ___o___      
    763                !  |   |___x___|   OR    |  |   x   | 
    764                !  v       o           (4)  |       | 
    765                IF( ijbi==0     .OR. ij1==0     .OR. ij2==0     .OR. ij3==0     )   lrecv_bdyint(ib_bdy,igrd,jpso,ir) = .TRUE. 
    766                IF( ijbe==0                                                     )   lrecv_bdyext(ib_bdy,igrd,jpso,ir) = .TRUE. 
    767                IF( ijbi==jpj+1 .OR. ij1==jpj+1 .OR. ij2==jpj+1 .OR. ij3==jpj+1 )   lrecv_bdyint(ib_bdy,igrd,jpno,ir) = .TRUE. 
    768                IF( ijbe==jpj+1                                                 )   lrecv_bdyext(ib_bdy,igrd,jpno,ir) = .TRUE. 
    769                ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo 
    770                !   ^  |    o    |                                                :         :  
    771                !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....| 
    772                !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |    
    773                ! ==> the neighbour cannot compute the point x -> need to send it 
    774                IF( ij ==     2*nn_hls   .AND. mpiSnei(nn_hls,jpso) > -1 ) THEN   ! 2*nn_hls      -> jj=jpj of southern neighbour  
    775                   IF( ijbi==ij+1 .OR. ij1==ij+1 .OR. ij2==ij+1 .OR. ij3==ij+1 )   lsend_bdyint(ib_bdy,igrd,jpso,ir) = .TRUE. 
    776                   IF( ijbe==ij+1                                              )   lsend_bdyext(ib_bdy,igrd,jpso,ir) = .TRUE. 
    777                ENDIF 
    778                IF( ij == jpj-2*nn_hls+1 .AND. mpiSnei(nn_hls,jpno) > -1 ) THEN   ! jpj-2*nn_hls+1-> jj=1   of northern neighbour 
    779                   IF( ijbi==ij-1 .OR. ij1==ij-1 .OR. ij2==ij-1 .OR. ij3==ij-1 )   lsend_bdyint(ib_bdy,igrd,jpno,ir) = .TRUE. 
    780                   IF( ijbe==ij-1                                              )   lsend_bdyext(ib_bdy,igrd,jpno,ir) = .TRUE. 
    781                ENDIF 
    782             END DO 
    783          END DO 
    784       END DO 
     879               DO icnt = 1, 4 
     880                  SELECT CASE( icnt ) 
     881                     !                                       ...|.... 
     882                  CASE( 1 )   ! x: rim on sw-corner            o|   : 
     883                     !          o: potential neighbour(s)      o|x__:__ 
     884                     !             outside of the MPI domain   o o o: 
     885                     !                                              : 
     886                     iRdiag    = jpsw            ;   iRsono    = jpso            ! Recv: for sw or so 
     887                     iSdiag    = jpne            ;   iSsono    = jpno            ! Send: to ne or no 
     888                     iiRst     = 1               ;   ijRst     = 1               ! Rcv sw-corner starting point 
     889                     iiRnd     = nn_hls          ;   ijRnd     = nn_hls          ! Rcv sw-corner   ending point 
     890                     iiSstdiag = Nie0-nn_hls+1   ;   ijSstdiag = Nje0-nn_hls+1   ! send to sw-corner of ne neighbourg 
     891                     iiSnddiag = Nie0            ;   ijSnddiag = Nje0            ! send to sw-corner of ne neighbourg 
     892                     iiSstsono = 1               ;   ijSstsono = Nje0-nn_hls+1   ! send to sw-corner of no neighbourg 
     893                     iiSndsono = nn_hls          ;   ijSndsono = Nje0            ! send to sw-corner of no neighbourg 
     894                     iioutdir  = -1              ;   ijoutdir  = -1              ! outside MPI domain: westward or southward 
     895                     !                                          ....|... 
     896                  CASE( 2 )   ! x: rim on se-corner             :   |o 
     897                     !          o: potential neighbour(s)     __:__x|o 
     898                     !             outside of the MPI domain    :o o o 
     899                     !                                          :     
     900                     iRdiag    = jpse            ;   iRsono    = jpso            ! Recv: for se or so 
     901                     iSdiag    = jpnw            ;   iSsono    = jpno            ! Send: to nw or no 
     902                     iiRst     = Nie0+1          ;   ijRst     = 1               ! Rcv se-corner starting point 
     903                     iiRnd     = jpi             ;   ijRnd     = nn_hls          ! Rcv se-corner   ending point 
     904                     iiSstdiag = Nis0            ;   ijSstdiag = Nje0-nn_hls+1   ! send to se-corner of nw neighbourg 
     905                     iiSnddiag = Nis0+nn_hls-1   ;   ijSnddiag = Nje0            ! send to se-corner of nw neighbourg 
     906                     iiSstsono = Nie0+1          ;   ijSstsono = Nje0-nn_hls+1   ! send to se-corner of no neighbourg 
     907                     iiSndsono = jpi             ;   ijSndsono = Nje0            ! send to se-corner of no neighbourg 
     908                     iioutdir  = 1               ;   ijoutdir  = -1              ! outside MPI domain: eastward or southward 
     909                     !                                              :        
     910                     !                                         o o_o:___ 
     911                  CASE( 3 )   ! x: rim on nw-corner            o|x  : 
     912                     !          o: potential neighbour(s)    ..o|...: 
     913                     !             outside of the MPI domain    | 
     914                     iRdiag    = jpnw            ;   iRsono    = jpno            ! Recv: for nw or no 
     915                     iSdiag    = jpse            ;   iSsono    = jpso            ! Send: to se or so 
     916                     iiRst     = 1               ;   ijRst     = Nje0+1          ! Rcv nw-corner starting point 
     917                     iiRnd     = nn_hls          ;   ijRnd     = jpj             ! Rcv nw-corner   ending point 
     918                     iiSstdiag = Nie0-nn_hls+1   ;   ijSstdiag = Njs0            ! send to nw-corner of se neighbourg 
     919                     iiSnddiag = Nie0            ;   ijSnddiag = Njs0+nn_hls-1   ! send to nw-corner of se neighbourg 
     920                     iiSstsono = 1               ;   ijSstsono = Njs0            ! send to nw-corner of so neighbourg 
     921                     iiSndsono = nn_hls          ;   ijSndsono = Njs0+nn_hls-1   ! send to nw-corner of so neighbourg 
     922                     iioutdir  = -1              ;   ijoutdir  =  1              ! outside MPI domain: westward or northward 
     923                     !                                          :        
     924                     !                                       ___:o_o o 
     925                  CASE( 4 )   ! x: rim on ne-corner             :  x|o 
     926                     !          o: potential neighbour(s)       :...|o... 
     927                     !             outside of the MPI domain        | 
     928                     iRdiag    = jpne            ;   iRsono    = jpno            ! Recv: for ne or no 
     929                     iSdiag    = jpsw            ;   iSsono    = jpso            ! Send: to sw or so 
     930                     iiRst     = Nie0+1          ;   ijRst     = Nje0+1          ! Rcv ne-corner starting point 
     931                     iiRnd     = jpi             ;   ijRnd     = jpj             ! Rcv ne-corner   ending point 
     932                     iiSstdiag = Nis0            ;   ijSstdiag = Njs0            ! send to ne-corner of sw neighbourg 
     933                     iiSnddiag = Nis0+nn_hls-1   ;   ijSnddiag = Njs0+nn_hls-1   ! send to ne-corner of sw neighbourg 
     934                     iiSstsono = Nie0+1          ;   ijSstsono = Njs0            ! send to ne-corner of so neighbourg 
     935                     iiSndsono = jpi             ;   ijSndsono = Njs0+nn_hls-1   ! send to ne-corner of so neighbourg 
     936                     iioutdir  = 1               ;   ijoutdir  = 1               ! outside MPI domain: eastward or southward 
     937                  END SELECT 
     938                  ! 
     939                  ! Check if we need to receive data for this rim point 
     940                  IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN   ! rim point on the corner 
     941                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain? 
     942                     ! take care of neighbourg(s) in the interior of the computational domain 
     943                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain 
     944                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it 
     945                        IF( mpiRnei(nn_hls,iRdiag) > -1                    )   lrecv_bdyint(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg 
     946                        IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 )   lrecv_bdyint(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg 
     947                     ENDIF 
     948                     ! take care of neighbourg in the exterior of the computational domain 
     949                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it 
     950                        IF( mpiRnei(nn_hls,iRdiag) > -1                    )   lrecv_bdyext(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg 
     951                        IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 )   lrecv_bdyext(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg 
     952                     ENDIF 
     953                  ENDIF 
     954                  ! 
     955                  ! Check if this rim point corresponds to the corner of one neighbourg. if yes, do we need to send data? 
     956                  ! Direct send to diag: Is this rim point the corner point of a diag neighbour with which we communicate? 
     957                  IF( ii >= iiSstdiag .AND. ii <= iiSnddiag .AND. ij >= ijSstdiag .AND. ij <= ijSnddiag   & 
     958                     &                .AND. mpiSnei(nn_hls,iSdiag) > -1 ) THEN 
     959                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain? 
     960                     ! take care of neighbourg(s) in the interior of the computational domain 
     961                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of diag nei MPI  
     962                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it 
     963                        &    lsend_bdyint(ib_bdy,igrd,iSdiag,ir) = .TRUE.                        ! send rim point data to diag nei 
     964                     ! take care of neighbourg in the exterior of the computational domain 
     965                     IF(  iibe==iiout .OR. ijbe==ijout )   &                                  
     966                        &    lsend_bdyext(ib_bdy,igrd,iSdiag,ir) = .TRUE. 
     967                  ENDIF 
     968                  ! Indirect send to diag (through so/no): rim point is the corner point of a so/no nei with which we communicate 
     969                  IF( ii >= iiSstsono .AND. ii <= iiSndsono .AND. ij >= ijSstsono .AND. ij <= ijSndsono   & 
     970                     &                .AND. mpiSnei(nn_hls,iSsono) > -1 .AND. nn_comm == 1 ) THEN 
     971                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain? 
     972                     ! take care of neighbourg(s) in the interior of the computational domain 
     973                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of so/no nei MPI 
     974                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it 
     975                        &    lsend_bdyint(ib_bdy,igrd,iSsono,ir) = .TRUE.                        ! send rim point data to so/no nei 
     976                     ! take care of neighbourg in the exterior of the computational domain 
     977                     IF(  iibe==iiout .OR. ijbe==ijout )   & 
     978                        &    lsend_bdyext(ib_bdy,igrd,iSsono,ir) = .TRUE. 
     979                  ENDIF 
     980                  ! 
     981               END DO   ! 4 corners 
     982            END DO   ! ib 
     983         END DO   ! igrd 
     984 
     985         ! Comment out for debug 
     986!!$         DO ir = 0,1 
     987!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   & 
     988!!$               &                              lsend = lsend_bdyint(ib_bdy,1,:,ir), lrecv = lrecv_bdyint(ib_bdy,1,:,ir) ) 
     989!!$            IF(lwp) WRITE(numout,*) ' bdy debug int T', ir ; CALL FLUSH(numout) 
     990!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   & 
     991!!$               &                              lsend = lsend_bdyint(ib_bdy,2,:,ir), lrecv = lrecv_bdyint(ib_bdy,2,:,ir) ) 
     992!!$            IF(lwp) WRITE(numout,*) ' bdy debug int U', ir ; CALL FLUSH(numout) 
     993!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   & 
     994!!$               &                              lsend = lsend_bdyint(ib_bdy,3,:,ir), lrecv = lrecv_bdyint(ib_bdy,3,:,ir) )     
     995!!$            IF(lwp) WRITE(numout,*) ' bdy debug int V', ir ; CALL FLUSH(numout) 
     996!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   & 
     997!!$               &                              lsend = lsend_bdyext(ib_bdy,1,:,ir), lrecv = lrecv_bdyext(ib_bdy,1,:,ir) ) 
     998!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext T', ir ; CALL FLUSH(numout) 
     999!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   & 
     1000!!$               &                              lsend = lsend_bdyext(ib_bdy,2,:,ir), lrecv = lrecv_bdyext(ib_bdy,2,:,ir) ) 
     1001!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext U', ir ; CALL FLUSH(numout) 
     1002!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   & 
     1003!!$               &                              lsend = lsend_bdyext(ib_bdy,3,:,ir), lrecv = lrecv_bdyext(ib_bdy,3,:,ir) )     
     1004!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext V', ir ; CALL FLUSH(numout) 
     1005!!$         END DO 
     1006          
     1007      END DO   ! ib_bdy 
    7851008 
    7861009      DO ib_bdy = 1,nb_bdy 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdytra.F90

    r15127 r15440  
    5555      TYPE(ztrabdy), DIMENSION(jpts) :: zdta                   ! Temporary data structure 
    5656      LOGICAL                        :: llrim0                 ! indicate if rim 0 is treated 
    57       LOGICAL, DIMENSION(4)          :: llsend1, llrecv1       ! indicate how communications are to be carried out 
     57      LOGICAL, DIMENSION(8)          :: llsend1, llrecv1       ! indicate how communications are to be carried out 
    5858      !!---------------------------------------------------------------------- 
    5959      igrd = 1 
    … …  
    9696               llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points 
    9797            CASE('orlanski', 'orlanski_npo') 
    98                llsend1(:) = llsend1(:) .OR. lsend_bdy(ib_bdy,1,:,ir)   ! possibly every direction, T points 
    99                llrecv1(:) = llrecv1(:) .OR. lrecv_bdy(ib_bdy,1,:,ir)   ! possibly every direction, T points 
     98               llsend1(:) = llsend1(:) .OR. lsend_bdyolr(ib_bdy,1,:,ir)   ! possibly every direction, T points 
     99               llrecv1(:) = llrecv1(:) .OR. lrecv_bdyolr(ib_bdy,1,:,ir)   ! possibly every direction, T points 
    100100            END SELECT 
    101101         END DO 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DYN/dynspg_ts.F90

    r15127 r15440  
    10511051#if defined key_agrif 
    10521052      ! Restrict the use of Agrif to the forward case only 
    1053 !!!      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 
     1053      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 
    10541054#endif 
    10551055      ! 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ICB/icbini.F90

    r15127 r15440  
    7373      ! 
    7474      IF( .NOT. ln_icebergs )   RETURN 
     75      ! 
     76      ALLOCATE( utau_icb(jpi,jpj), vtau_icb(jpi,jpj) ) 
    7577      ! 
    7678      !                          ! allocate gridded fields 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ICB/icbutl.F90

    r15127 r15440  
    9393      sss_e(1:jpi,1:jpj) = sss_m(:,:) 
    9494      fr_e (1:jpi,1:jpj) = fr_i (:,:) 
    95       ua_e (1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    96       va_e (1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
     95      ua_e (1:jpi,1:jpj) = utau_icb (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
     96      va_e (1:jpi,1:jpj) = vtau_icb (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 
    9797      ff_e(1:jpi,1:jpj) = ff_f (:,:)  
    9898      ! 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbc_lnk_call_generic.h90

    r15127 r15440  
    5252      REAL(PRECISION)      , OPTIONAL        , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries) 
    5353      INTEGER              , OPTIONAL        , INTENT(in   ) ::   khls        ! halo size, default = nn_hls 
    54       LOGICAL, DIMENSION(4), OPTIONAL        , INTENT(in   ) ::   lsend, lrecv   ! indicate how communications are to be carried out 
     54      LOGICAL, DIMENSION(8), OPTIONAL        , INTENT(in   ) ::   lsend, lrecv   ! indicate how communications are to be carried out 
    5555      LOGICAL              , OPTIONAL        , INTENT(in   ) ::   ld4only     ! if .T., do only 4-neighbour comm (ignore corners) 
    5656      !! 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbc_lnk_neicoll_generic.h90

    r15349 r15440  
    99      REAL(PRECISION),      OPTIONAL, INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries) 
    1010      INTEGER ,             OPTIONAL, INTENT(in   ) ::   khls        ! halo size, default = nn_hls 
    11       LOGICAL, DIMENSION(4),OPTIONAL, INTENT(in   ) ::   lsend, lrecv  ! communication with other 4 proc 
     11      LOGICAL, DIMENSION(8),OPTIONAL, INTENT(in   ) ::   lsend, lrecv  ! communication with other 4 proc 
    1212      LOGICAL,              OPTIONAL, INTENT(in   ) ::   ld4only     ! if .T., do only 4-neighbour comm (ignore corners) 
    1313      ! 
    … …  
    7474      ! define llsend and llrecv: logicals which say if mpi-neibourgs for send or receive exist or not. 
    7575      IF     ( PRESENT(lsend) .AND. PRESENT(lrecv) ) THEN   ! localy defined neighbourgs  
    76          CALL ctl_stop( 'STOP', 'mpp_nc_generic+BDY not yet implemented') 
    77 !!$         ---> llsend(:) = lsend(:)   ;   llrecv(:) = lrecv(:) ??? 
     76         CALL ctl_stop( 'STOP', 'mpp_nc_generic+lsend and lrecv not yet implemented') 
    7877      ELSE IF( PRESENT(lsend) .OR.  PRESENT(lrecv) ) THEN 
    7978         WRITE(ctmp1,*) TRIM(cdname), '  is calling lbc_lnk with only one of the two arguments lsend or lrecv' 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbc_lnk_pt2pt_generic.h90

    r15349 r15440  
    1010      REAL(PRECISION),      OPTIONAL, INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries) 
    1111      INTEGER ,             OPTIONAL, INTENT(in   ) ::   khls        ! halo size, default = nn_hls 
    12       LOGICAL, DIMENSION(4),OPTIONAL, INTENT(in   ) ::   lsend, lrecv  ! communication with other 4 proc 
     12      LOGICAL, DIMENSION(8),OPTIONAL, INTENT(in   ) ::   lsend, lrecv  ! communication with other 4 proc 
    1313      LOGICAL,              OPTIONAL, INTENT(in   ) ::   ld4only     ! if .T., do only 4-neighbour comm (ignore corners) 
    1414      ! 
    … …  
    7272      ! define llsend and llrecv: logicals which say if mpi-neibourgs for send or receive exist or not. 
    7373      IF     ( PRESENT(lsend) .AND. PRESENT(lrecv) ) THEN   ! localy defined neighbourgs  
    74          llsend(1:4) = lsend(1:4)   ;   llrecv(1:4) = lrecv(1:4) 
     74         llsend(:) = lsend(:)   ;   llrecv(:) = lrecv(:) 
    7575      ELSE IF( PRESENT(lsend) .OR.  PRESENT(lrecv) ) THEN 
    7676         WRITE(ctmp1,*) TRIM(cdname), ' is calling lbc_lnk with only one of the two arguments lsend or lrecv' 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbc_oce.F90

    r14227 r15440  
    106106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2] 
    107107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2] 
     108   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs   [N/m2] 
    108109   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
    109110   !! wndm is used compute surface gases exchanges in ice-free ocean or leads 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r14072 r15440  
    3131   USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer 
    3232   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 
     33   USE sbcwave, ONLY : charn 
     34   USE sbc_oce, ONLY : ln_charn ! wave module 
    3335 
    3436   IMPLICIT NONE 
    … …  
    229231      u_star = 0.035_wp*Ubzu*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
    230232 
    231       z0     = charn0_ecmwf*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     233      IF (ln_charn) THEN ! Charnock value if wave coupling 
     234         z0 = charn*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     235      ELSE 
     236         z0 = charn0_ecmwf*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     237      ENDIF 
    232238      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    233239 
    … …  
    296302         ztmp2  = u_star*u_star 
    297303         ztmp1  = znu_a/u_star 
     304         IF (ln_charn) THEN ! Charnock value if wave coupling 
     305            z0  = MIN( ABS( alpha_M*ztmp1 + charn*ztmp2/grav ) , 0.001_wp) 
     306         ELSE 
     307            z0  = MIN( ABS( alpha_M*ztmp1 + charn0_ecmwf*ztmp2/grav ) , 0.001_wp) 
     308         ENDIF 
    298309         z0     = MIN( ABS( alpha_M*ztmp1 + charn0_ecmwf*ztmp2/grav ) , 0.001_wp) 
    299310         z0t    = MIN( ABS( alpha_H*ztmp1                           ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcfwb.F90

    r15127 r15440  
    3636 
    3737   REAL(wp) ::   rn_fwb0   ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only) 
    38    REAL(wp) ::   a_fwb     ! annual domain averaged freshwater budget from the 
    39                            ! previous year 
     38   REAL(wp) ::   a_fwb     ! annual domain averaged freshwater budget from the previous year 
     39   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget from the year before or at initial state 
     40   REAL(wp) ::   a_fwb_ini ! initial domain averaged freshwater budget 
    4041   REAL(wp) ::   area      ! global mean ocean surface (interior domain) 
    4142 
    … …  
    129130         ! 
    130131      CASE ( 2 )                             !==  fw adjustment based on fw budget at the end of the previous year  ==! 
    131          ! 
    132          IF( kt == nit000 ) THEN                                                                    ! initialisation 
    133             !                                                                                       ! set the fw adjustment (a_fwb) 
    134             IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb',   ldstop = .FALSE. ) > 0 ) THEN        !    as read from restart file 
    135                IF(lwp) WRITE(numout,*) 'sbc_fwb : reading FW-budget adjustment from restart file' 
    136                CALL iom_get( numror, 'a_fwb',   a_fwb ) 
    137             ELSE                                                                                    !    as specified in namelist 
    138                a_fwb = rn_fwb0 
     132         !                                                simulation is supposed to start 1st of January 
     133         IF( kt == nit000 ) THEN                                                                 ! initialisation 
     134            !                                                                                    ! set the fw adjustment (a_fwb) 
     135            IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0     &      !    as read from restart file 
     136               &           .AND. iom_varid( numror, 'a_fwb',   ldstop = .FALSE. ) > 0 ) THEN 
     137               IF(lwp)   WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file' 
     138               CALL iom_get( numror, 'a_fwb_b', a_fwb_b ) 
     139               CALL iom_get( numror, 'a_fwb'  , a_fwb ) 
     140               ! 
     141               a_fwb_ini = a_fwb_b 
     142            ELSE                                                                                 !    as specified in namelist 
     143               IF(lwp)   WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0' 
     144               a_fwb   = rn_fwb0 
     145               a_fwb_b = 0._wp   ! used only the first year then it is replaced by a_fwb_ini 
     146               ! 
     147               a_fwb_ini = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) & 
     148                  &      * rho0 / ( area * rday * REAL(nyear_len(1), wp) ) 
    139149            END IF 
    140150            ! 
    141             IF(lwp)WRITE(numout,*) 
    142             IF(lwp)WRITE(numout,*)'sbc_fwb : initial freshwater-budget adjustment = ', a_fwb, 'kg/m2/s' 
    143             ! 
    144          ENDIF    
    145          !                                         ! Update a_fwb if new year start 
    146          ikty = 365 * 86400 / rn_Dt                  !!bug  use of 365 days leap year or 360d year !!!!!!! 
    147          IF( MOD( kt, ikty ) == 0 ) THEN 
    148                                                       ! mean sea level taking into account the ice+snow 
    149                                                       ! sum over the global domain 
    150             a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) 
    151             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s 
    152 !!gm        !                                                      !!bug 365d year  
    153          ENDIF 
    154          !  
    155          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes 
    156             zcoef = a_fwb * rcp 
    157             emp(:,:) = emp(:,:) + a_fwb              * tmask(:,:,1) 
    158             qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
     151            IF(lwp)   WRITE(numout,*) 
     152            IF(lwp)   WRITE(numout,*)'sbc_fwb : freshwater-budget at the end of previous year = ', a_fwb    , 'kg/m2/s' 
     153            IF(lwp)   WRITE(numout,*)'          freshwater-budget at initial state            = ', a_fwb_ini, 'kg/m2/s' 
     154            ! 
     155         ELSE 
     156            ! at the end of year n: 
     157            ikty = nyear_len(1) * 86400 / NINT(rn_Dt) 
     158            IF( MOD( kt, ikty ) == 0 ) THEN   ! Update a_fwb at the last time step of a year 
     159               !                                It should be the first time step of a year MOD(kt-1,ikty) but then the restart would be wrong 
     160               !                                Hence, we make a small error here but the code is restartable 
     161               a_fwb_b = a_fwb_ini 
     162               ! mean sea level taking into account ice+snow 
     163               a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) 
     164               a_fwb   = a_fwb * rho0 / ( area * rday * REAL(nyear_len(1), wp) )   ! convert in kg/m2/s 
     165            ENDIF 
     166            ! 
     167         ENDIF 
     168         ! 
     169         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes using previous year budget minus initial state 
     170            zcoef = ( a_fwb - a_fwb_b ) 
     171            emp(:,:) = emp(:,:) + zcoef * tmask(:,:,1) 
     172            qns(:,:) = qns(:,:) - zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
    159173            ! outputs 
    160             IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zcoef * sst_m(:,:) * tmask(:,:,1) ) 
    161             IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -a_fwb              * tmask(:,:,1) ) 
     174            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ) 
     175            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -zcoef * tmask(:,:,1) ) 
    162176         ENDIF 
    163177         ! Output restart information 
    … …  
    166180            IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt 
    167181            IF(lwp) WRITE(numout,*) '~~~~' 
    168             CALL iom_rstput( kt, nitrst, numrow, 'a_fwb',   a_fwb ) 
     182            CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b ) 
     183            CALL iom_rstput( kt, nitrst, numrow, 'a_fwb',   a_fwb   ) 
    169184         END IF 
    170185         ! 
    171          IF( kt == nitend .AND. lwp ) WRITE(numout,*) 'sbc_fwb : final freshwater-budget adjustment = ', a_fwb, 'kg/m2/s' 
     186         IF( kt == nitend .AND. lwp ) THEN 
     187            WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', a_fwb  , 'kg/m2/s' 
     188            WRITE(numout,*) '          freshwater-budget at initial state                    = ', a_fwb_b, 'kg/m2/s' 
     189         ENDIF 
    172190         ! 
    173191      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==! 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcmod.F90

    r15349 r15440  
    466466      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. ) 
    467467      ! 
     468      IF( ln_icebergs ) THEN  ! save pure stresses (with no ice-ocean stress) for use by icebergs 
     469         utau_icb(:,:) = utau(:,:) ; vtau_icb(:,:) = vtau(:,:)  
     470      ENDIF 
     471      ! 
    468472      !                                            !==  Misc. Options  ==! 
    469473      ! 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/lib_fortran.F90

    r15349 r15440  
    3434   PUBLIC   glob_sum_vec 
    3535   PUBLIC   glob_sum_full_vec 
     36   PUBLIC   glob_min_vec, glob_max_vec 
    3637#if defined key_nosignedzero 
    3738   PUBLIC SIGN 
    … …  
    6162   INTERFACE glob_sum_full_vec 
    6263      MODULE PROCEDURE glob_sum_full_vec_3d, glob_sum_full_vec_4d 
     64   END INTERFACE 
     65   INTERFACE glob_min_vec 
     66      MODULE PROCEDURE glob_min_vec_3d, glob_min_vec_4d 
     67   END INTERFACE 
     68   INTERFACE glob_max_vec 
     69      MODULE PROCEDURE glob_max_vec_3d, glob_max_vec_4d 
    6370   END INTERFACE 
    6471 
    … …  
    503510   END FUNCTION glob_sum_full_vec_4d 
    504511 
     512   FUNCTION glob_min_vec_3d( cdname, ptab ) RESULT( ptmp ) 
     513      !!---------------------------------------------------------------------- 
     514      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine 
     515      REAL(wp),          INTENT(in) ::   ptab(:,:,:)   ! array on which operation is applied 
     516      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp 
     517      ! 
     518      INTEGER     ::   jk    ! dummy loop indice & dimension 
     519      INTEGER     ::   ipk   ! dimension 
     520      !!----------------------------------------------------------------------- 
     521      ! 
     522      ipk = SIZE(ptab,3) 
     523      DO jk = 1, ipk 
     524         ptmp(jk) = MINVAL( ptab(:,:,jk) * tmask_i(:,:) ) 
     525      ENDDO 
     526      ! 
     527      CALL mpp_min( cdname, ptmp (:) ) 
     528      ! 
     529   END FUNCTION glob_min_vec_3d 
     530 
     531   FUNCTION glob_min_vec_4d( cdname, ptab ) RESULT( ptmp ) 
     532      !!---------------------------------------------------------------------- 
     533      CHARACTER(len=*),  INTENT(in) ::   cdname          ! name of the calling subroutine 
     534      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:)   ! array on which operation is applied 
     535      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp 
     536      ! 
     537      INTEGER     ::   jk , jl    ! dummy loop indice & dimension 
     538      INTEGER     ::   ipk, ipl   ! dimension 
     539      !!----------------------------------------------------------------------- 
     540      ! 
     541      ipk = SIZE(ptab,3) 
     542      ipl = SIZE(ptab,4) 
     543      DO jl = 1, ipl 
     544            ptmp(jl) = MINVAL( ptab(:,:,1,jl) * tmask_i(:,:) )          
     545         DO jk = 2, ipk 
     546            ptmp(jl) = MIN( ptmp(jl), MINVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) ) 
     547         ENDDO 
     548      ENDDO 
     549      ! 
     550      CALL mpp_min( cdname, ptmp (:) ) 
     551      ! 
     552   END FUNCTION glob_min_vec_4d 
     553    
     554   FUNCTION glob_max_vec_3d( cdname, ptab ) RESULT( ptmp ) 
     555      !!---------------------------------------------------------------------- 
     556      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine 
     557      REAL(wp),          INTENT(in) ::   ptab(:,:,:)   ! array on which operation is applied 
     558      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp 
     559      ! 
     560      INTEGER     ::   jk    ! dummy loop indice & dimension 
     561      INTEGER     ::   ipk   ! dimension 
     562      !!----------------------------------------------------------------------- 
     563      ! 
     564      ipk = SIZE(ptab,3) 
     565      DO jk = 1, ipk 
     566         ptmp(jk) = MAXVAL( ptab(:,:,jk) * tmask_i(:,:) ) 
     567      ENDDO 
     568      ! 
     569      CALL mpp_max( cdname, ptmp (:) ) 
     570      ! 
     571   END FUNCTION glob_max_vec_3d 
     572 
     573   FUNCTION glob_max_vec_4d( cdname, ptab ) RESULT( ptmp ) 
     574      !!---------------------------------------------------------------------- 
     575      CHARACTER(len=*),  INTENT(in) ::   cdname          ! name of the calling subroutine 
     576      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:)   ! array on which operation is applied 
     577      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp 
     578      ! 
     579      INTEGER     ::   jk , jl    ! dummy loop indice & dimension 
     580      INTEGER     ::   ipk, ipl   ! dimension 
     581      !!----------------------------------------------------------------------- 
     582      ! 
     583      ipk = SIZE(ptab,3) 
     584      ipl = SIZE(ptab,4) 
     585      DO jl = 1, ipl 
     586            ptmp(jl) = MAXVAL( ptab(:,:,1,jl) * tmask_i(:,:) )          
     587         DO jk = 2, ipk 
     588            ptmp(jl) = MAX( ptmp(jl), MAXVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) ) 
     589         ENDDO 
     590      ENDDO 
     591      ! 
     592      CALL mpp_max( cdname, ptmp (:) ) 
     593      ! 
     594   END FUNCTION glob_max_vec_4d 
     595    
    505596   SUBROUTINE DDPDD( ydda, yddb ) 
    506597      !!---------------------------------------------------------------------- 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/step.F90

    r15127 r15440  
    229229         IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
    230230                  &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment 
     231         IF( ln_bkgwri )    CALL asm_bkg_wri( kstp, Nnn )     ! output background fields 
    231232         IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
    232233#if defined key_agrif 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/stpmlf.F90

    r15127 r15440  
    241241         IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
    242242                  &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment 
     243         IF( ln_bkgwri )    CALL asm_bkg_wri( kstp, Nnn )     ! output background fields 
    243244         IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
    244245#if defined key_agrif 
    … …  
    507508         ! so that asselin contribution is removed at the same time 
    508509         DO jk = 1, jpkm1 
    509             DO_2D( 0, 0, 0, 0 ) 
    510                puu(ji,jj,jk,Kmm) = ( puu(ji,jj,jk,Kmm) - un_adv(ji,jj)*r1_hu(ji,jj,Kmm) + uu_b(ji,jj,Kmm) )*umask(ji,jj,jk) 
    511                pvv(ji,jj,jk,Kmm) = ( pvv(ji,jj,jk,Kmm) - vn_adv(ji,jj)*r1_hv(ji,jj,Kmm) + vv_b(ji,jj,Kmm) )*vmask(ji,jj,jk) 
    512             END_2D 
     510            puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
     511            pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    513512         END DO 
    514513      ENDIF 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/TRP/trcsbc.F90

    r14215 r15440  
    5151      !!            The surface freshwater flux modify the ocean volume 
    5252      !!         and thus the concentration of a tracer as : 
    53       !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t_   for k=1 
    54       !!         where emp, the surface freshwater budget (evaporation minus 
    55       !!         precipitation ) given in kg/m2/s is divided 
    56       !!         by 1035 kg/m3 (density of ocean water) to obtain m/s. 
     53      !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t_ + fmmflx * tri / e3t   for k=1 
     54      !!          - tr(Kmm) , the concentration of tracer in the ocean 
     55      !!          - tri, the concentration of tracer in the sea-ice 
     56      !!          - emp, the surface freshwater budget (evaporation minus precipitation + fmmflx) 
     57      !!            given in kg/m2/s is divided by 1035 kg/m3 (density of ocean water) to obtain m/s. 
     58      !!          - fmmflx, the flux asscociated to freezing-melting of sea-ice  
     59      !!            In linear free surface case (ln_linssh=T), the volume of the 
     60      !!            ocean does not change with the water exchanges at the (air+ice)-sea 
    5761      !! 
    5862      !! ** Action  : - Update the 1st level of tr(:,:,:,:,Krhs) with the trend associated 
    … …  
    6670      INTEGER  ::   ji, jj, jn                      ! dummy loop indices 
    6771      REAL(wp) ::   zse3t, zrtrn, zfact     ! local scalars 
    68       REAL(wp) ::   zftra, zdtra, ztfx, ztra   !   -      - 
     72      REAL(wp) ::   zdtra          !   -      - 
    6973      CHARACTER (len=22) :: charout 
    70       REAL(wp), DIMENSION(jpi,jpj)   ::   zsfx 
    7174      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd 
    7275      !!--------------------------------------------------------------------- 
    … …  
    106109      ENDIF 
    107110 
    108       ! Coupling online : river runoff is added to the horizontal divergence (hdiv) in the subroutine sbc_rnf_div  
    109       ! one only consider the concentration/dilution effect due to evaporation minus precipitation + freezing/melting of sea-ice 
    110       ! Coupling offline : runoff are in emp which contains E-P-R 
    111       ! 
    112       IF( .NOT.ln_linssh ) THEN  ! online coupling with vvl 
    113          zsfx(:,:) = 0._wp 
    114       ELSE                                      ! online coupling free surface or offline with free surface 
    115          zsfx(:,:) = emp(:,:) 
    116       ENDIF 
    117  
    118111      ! 0. initialization 
    119112      SELECT CASE ( nn_ice_tr ) 
    120113 
    121       CASE ( -1 ) ! No tracers in sea ice (null concentration in sea ice) 
    122          ! 
    123          DO jn = 1, jptra 
    124             DO_2D( 0, 0, 0, 1 ) 
    125                sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    126             END_2D 
    127          END DO 
    128          ! 
    129       CASE ( 0 )  ! Same concentration in sea ice and in the ocean 
    130          ! 
    131          DO jn = 1, jptra 
    132             DO_2D( 0, 0, 0, 1 ) 
    133                sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    134             END_2D 
    135          END DO 
     114      CASE ( -1 ) ! ! No tracers in sea ice ( trc_i = 0 ) 
     115         ! 
     116         DO jn = 1, jptra 
     117            DO_2D( 0, 0, 0, 1 ) 
     118               sbc_trc(ji,jj,jn) = 0._wp 
     119            END_2D 
     120         END DO 
     121         ! 
     122         IF( ln_linssh ) THEN  !* linear free surface   
     123            DO jn = 1, jptra 
     124               DO_2D( 0, 0, 0, 1 ) 
     125                  sbc_trc(ji,jj,jn) = sbc_trc(ji,jj,jn) + r1_rho0 * emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) !==>> add concentration/dilution effect due to constant volume cell 
     126               END_2D 
     127            END DO 
     128         ENDIF 
     129         ! 
     130      CASE ( 0 )  ! Same concentration in sea ice and in the ocean ( trc_i = ptr(...,Kmm)  ) 
     131         ! 
     132         DO jn = 1, jptra 
     133            DO_2D( 0, 0, 0, 1 ) 
     134               sbc_trc(ji,jj,jn) = - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
     135            END_2D 
     136         END DO 
     137         ! 
     138         IF( ln_linssh ) THEN  !* linear free surface   
     139            DO jn = 1, jptra 
     140               DO_2D( 0, 0, 0, 1 ) 
     141                  sbc_trc(ji,jj,jn) = sbc_trc(ji,jj,jn) + r1_rho0 * emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) !==>> add concentration/dilution effect due to constant volume cell 
     142               END_2D 
     143            END DO 
     144         ENDIF 
    136145         ! 
    137146      CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice  
    … …  
    139148         DO jn = 1, jptra 
    140149            DO_2D( 0, 0, 0, 1 ) 
    141                zse3t = 1. / e3t(ji,jj,1,Kmm) 
    142                ! tracer flux at the ice/ocean interface (tracer/m2/s) 
    143                zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
    144                !                                         ! only used in the levitating sea ice case 
    145                ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux 
    146                ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux 
    147                ztfx  = zftra                        ! net tracer flux 
    148                ! 
    149                zdtra = r1_rho0 * ( ztfx + ( zsfx(ji,jj) + fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) )  
    150                IF ( zdtra < 0. ) THEN 
    151                   zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise 
    152                ENDIF 
    153                sbc_trc(ji,jj,jn) =  zdtra  
    154             END_2D 
    155          END DO 
     150               sbc_trc(ji,jj,jn) = - fmmflx(ji,jj) * r1_rho0 * trc_i(ji,jj,jn) 
     151            END_2D 
     152         END DO 
     153         ! 
     154         IF( ln_linssh ) THEN  !* linear free surface   
     155            DO jn = 1, jptra 
     156               DO_2D( 0, 0, 0, 1 ) 
     157                  sbc_trc(ji,jj,jn) = sbc_trc(ji,jj,jn) + r1_rho0 * emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) !==>> add concentration/dilution effect due to constant volume cell 
     158               END_2D 
     159            END DO 
     160         ENDIF 
     161         ! 
     162         DO jn = 1, jptra 
     163            DO_2D( 0, 0, 0, 1 ) 
     164               zse3t = rDt_trc / e3t(ji,jj,1,Kmm) 
     165               zdtra = ptr(ji,jj,1,jn,Kmm) + sbc_trc(ji,jj,jn) * zse3t  
     166               IF( zdtra < 0. ) sbc_trc(ji,jj,jn) = MAX( zdtra, -ptr(ji,jj,1,jn,Kmm) / zse3t  ) ! avoid negative concentration that can occurs if trc_i > ptr  
     167            END_2D 
     168         END DO 
     169         !                              
    156170      END SELECT 
    157171      ! 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/trc.F90

    r15127 r15440  
    116116   ! 
    117117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   trc3d   !: 3D diagnostics for tracers 
    118    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   trc2d   !: 2D diagnostics for tracers 
     118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   trc2d   !: 2D diagnostics for tracers 
    119119 
    120120   !! information for inputs 
    … …  
    179179      IF (jp_dia3d > 0 )   ALLOCATE( trc3d(jpi,jpj,jpk,jp_dia3d), STAT = ierr(3) ) 
    180180      ! 
    181       IF (jp_dia2d > 0 )   ALLOCATE( trc2d(jpi,jpj,jpk,jp_dia2d), STAT = ierr(4) ) 
     181      IF (jp_dia2d > 0 )   ALLOCATE( trc2d(jpi,jpj,jp_dia2d)    , STAT = ierr(4) ) 
    182182      !  
    183183      trc_alloc = MAXVAL( ierr ) 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/trcbdy.F90

    r13527 r15440  
    5050      REAL(wp), POINTER, DIMENSION(:,:) ::  ztrc 
    5151      LOGICAL                           :: llrim0               ! indicate if rim 0 is treated 
    52       LOGICAL, DIMENSION(4)             :: llsend1, llrecv1     ! indicate how communications are to be carried out 
     52      LOGICAL, DIMENSION(8)             :: llsend1, llrecv1     ! indicate how communications are to be carried out 
    5353      !!---------------------------------------------------------------------- 
    5454      ! 
    … …  
    9898               llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points 
    9999            CASE('orlanski','orlanski_npo') 
    100                llsend1(:) = llsend1(:) .OR. lsend_bdy(ib_bdy,1,:,ir)   ! possibly every direction, T points 
    101                llrecv1(:) = llrecv1(:) .OR. lrecv_bdy(ib_bdy,1,:,ir)   ! possibly every direction, T points 
     100               llsend1(:) = llsend1(:) .OR. lsend_bdyolr(ib_bdy,1,:,ir)   ! possibly every direction, T points 
     101               llrecv1(:) = llrecv1(:) .OR. lrecv_bdyolr(ib_bdy,1,:,ir)   ! possibly every direction, T points 
    102102            END SELECT 
    103103         END DO 
Note: See TracChangeset for help on using the changeset viewer.