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

Ignore:
Timestamp:
2016-11-21T10:38:43+01:00 (8 years ago)
Author:
flavoni
Message:

update branch CNRS-2016 to trunk 6720

Location:
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC
Files:
24 edited

Legend:

Unmodified
Added
Removed

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r6140 r7278  
    212212      REAL(wp) ::   zztmp   
    213213      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity 
    214       ! reading initial file 
    215       LOGICAL  ::   ln_tsd_init      !: T & S data flag 
    216       LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag 
    217       CHARACTER(len=100)            ::   cn_dir 
    218       TYPE(FLD_N)                   ::  sn_tem,sn_sal 
    219       INTEGER  ::   ios=0 
    220  
    221       NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal 
    222       ! 
    223  
    224       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
    225       READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
    226 901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp ) 
    227       REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
    228       READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
    229 902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp ) 
    230       IF(lwm) WRITE ( numond, namtsd ) 
    231214      ! 
    232215      !!---------------------------------------------------------------------- 
    … …  
    250233      IF( lk_mpp )   CALL mpp_sum( vol0 ) 
    251234 
    252       CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum ) 
    253       CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  ) 
    254       CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 ) 
     235 
     236      CALL iom_open ( 'sali_ref_clim_monthly', inum ) 
     237      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  ) 
     238      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 ) 
    255239      CALL iom_close( inum ) 
     240 
    256241      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
    257242      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90

    r6075 r7278  
    1717   USE in_out_manager 
    1818   USE sbc_oce 
     19   USE lib_mpp 
    1920   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2021    

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r7277 r7278  
    196196         CALL dom_stiff( zprt ) 
    197197         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )      !    ! Max. grid stiffness ratio 
    198       ENDIF 
     198      ENDIF 
    199199      ! 
    200200      !                                     ! ============================ 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r7277 r7278  
    114114      CASE (30)   ;   CALL xios_set_context_attr(TRIM(clname), calendar_type= "D360") 
    115115      END SELECT 
    116       WRITE(cldate,"(i4.4,'-',i2.2,'-',i2.2,' 00:00:00')") nyear,nmonth,nday  
     116      WRITE(cldate,"(i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':',i2.2,':00')") nyear,nmonth,nday,nhour,nminute 
    117117      CALL xios_set_context_attr(TRIM(clname), start_date=cldate ) 
    118118 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90

    r6140 r7278  
    99   !!            3.5  ! 2012     (S.Mocavero, I. Epicoco) optimization of BDY comm. via lbc_bdy_lnk and lbc_obc_lnk 
    1010   !!            3.4  ! 2012-12  (R. Bourdalle-Badie, G. Reffray)  add a C1D case   
     11   !!            3.6  ! 2015-06  (O. Tintó and M. Castrillo)  add lbc_lnk_multi   
    1112   !!---------------------------------------------------------------------- 
    1213#if defined key_mpp_mpi 
    … …  
    2223 
    2324   INTERFACE lbc_lnk_multi 
    24       MODULE PROCEDURE mpp_lnk_2d_9 
     25      MODULE PROCEDURE mpp_lnk_2d_9, mpp_lnk_2d_multiple 
    2526   END INTERFACE 
    2627   ! 
    … …  
    2930   END INTERFACE 
    3031   ! 
    31 !JMM interface not defined if not key_mpp_mpi : likely do not compile without this CPP key !!!! 
    3232   INTERFACE lbc_sum 
    3333      MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d 
    3434   END INTERFACE 
    35  
     35   ! 
    3636   INTERFACE lbc_bdy_lnk 
    3737      MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d 
    … …  
    8383   ! 
    8484   INTERFACE lbc_sum 
    85       MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d 
     85      MODULE PROCEDURE lbc_lnk_sum_3d, lbc_lnk_sum_2d 
    8686   END INTERFACE 
    8787 
    … …  
    9090   END INTERFACE 
    9191   ! 
     92   INTERFACE lbc_lnk_multi 
     93      MODULE PROCEDURE lbc_lnk_2d_9, lbc_lnk_2d_multiple 
     94   END INTERFACE 
     95 
    9296   INTERFACE lbc_bdy_lnk 
    9397      MODULE PROCEDURE lbc_bdy_lnk_2d, lbc_bdy_lnk_3d 
    … …  
    97101      MODULE PROCEDURE lbc_lnk_2d_e 
    98102   END INTERFACE 
     103    
     104   TYPE arrayptr 
     105      REAL , DIMENSION (:,:),  POINTER :: pt2d 
     106   END TYPE arrayptr 
     107   PUBLIC   arrayptr 
    99108 
    100109   PUBLIC   lbc_lnk       ! ocean/ice  lateral boundary conditions 
     110   PUBLIC   lbc_sum       ! ocean/ice  lateral boundary conditions (sum of the overlap region) 
    101111   PUBLIC   lbc_lnk_e     ! 
     112   PUBLIC   lbc_lnk_multi ! modified ocean lateral boundary conditions 
    102113   PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions 
    103114   PUBLIC   lbc_lnk_icb   ! 
    … …  
    181192      ! 
    182193   END SUBROUTINE lbc_lnk_2d 
     194    
     195   SUBROUTINE lbc_lnk_2d_multiple( pt2d_array , type_array , psgn_array , num_fields ) 
     196      !! 
     197      INTEGER :: num_fields 
     198      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array 
     199      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   type_array   ! define the nature of ptab array grid-points 
     200      !                                                               ! = T , U , V , F , W and I points 
     201      REAL(wp)        , DIMENSION(:), INTENT(in   ) ::   psgn_array   ! =-1 the sign change across the north fold boundary 
     202      !                                                               ! =  1. , the sign is kept 
     203      ! 
     204      INTEGER  ::   ii    !!MULTI SEND DUMMY LOOP INDICES 
     205      ! 
     206      DO ii = 1, num_fields 
     207        CALL lbc_lnk_2d( pt2d_array(ii)%pt2d, type_array(ii), psgn_array(ii) ) 
     208      END DO      
     209      ! 
     210   END SUBROUTINE lbc_lnk_2d_multiple 
     211 
     212   SUBROUTINE lbc_lnk_2d_9( pt2dA, cd_typeA, psgnA, pt2dB, cd_typeB, psgnB, pt2dC, cd_typeC, psgnC   & 
     213      &                   , pt2dD, cd_typeD, psgnD, pt2dE, cd_typeE, psgnE, pt2dF, cd_typeF, psgnF   & 
     214      &                   , pt2dG, cd_typeG, psgnG, pt2dH, cd_typeH, psgnH, pt2dI, cd_typeI, psgnI, cd_mpp, pval) 
     215      !!--------------------------------------------------------------------- 
     216      ! Second 2D array on which the boundary condition is applied 
     217      REAL(wp), DIMENSION(jpi,jpj), TARGET          , INTENT(inout) ::   pt2dA 
     218      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dB , pt2dC , pt2dD , pt2dE 
     219      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dF , pt2dG , pt2dH , pt2dI 
     220      ! define the nature of ptab array grid-points 
     221      CHARACTER(len=1)                              , INTENT(in   ) ::   cd_typeA 
     222      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeB , cd_typeC , cd_typeD , cd_typeE 
     223      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeF , cd_typeG , cd_typeH , cd_typeI 
     224      ! =-1 the sign change across the north fold boundary 
     225      REAL(wp)                                      , INTENT(in   ) ::   psgnA 
     226      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnB , psgnC , psgnD , psgnE 
     227      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnF , psgnG , psgnH , psgnI 
     228      CHARACTER(len=3)                    , OPTIONAL, INTENT(in   ) ::   cd_mpp   ! fill the overlap area only 
     229      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   pval     ! background value (used at closed boundaries) 
     230      !! 
     231      !!--------------------------------------------------------------------- 
     232 
     233      !!The first array 
     234      CALL lbc_lnk( pt2dA, cd_typeA, psgnA )  
     235 
     236      !! Look if more arrays to process 
     237      IF(PRESENT (psgnB) )CALL lbc_lnk( pt2dA, cd_typeA, psgnA )  
     238      IF(PRESENT (psgnC) )CALL lbc_lnk( pt2dC, cd_typeC, psgnC )  
     239      IF(PRESENT (psgnD) )CALL lbc_lnk( pt2dD, cd_typeD, psgnD )  
     240      IF(PRESENT (psgnE) )CALL lbc_lnk( pt2dE, cd_typeE, psgnE )  
     241      IF(PRESENT (psgnF) )CALL lbc_lnk( pt2dF, cd_typeF, psgnF )  
     242      IF(PRESENT (psgnG) )CALL lbc_lnk( pt2dG, cd_typeG, psgnG )  
     243      IF(PRESENT (psgnH) )CALL lbc_lnk( pt2dH, cd_typeH, psgnH )  
     244      IF(PRESENT (psgnI) )CALL lbc_lnk( pt2dI, cd_typeI, psgnI )  
     245 
     246   END SUBROUTINE lbc_lnk_2d_9 
     247 
     248 
     249 
     250 
    183251 
    184252#else 
    … …  
    379447      !     
    380448   END SUBROUTINE lbc_lnk_2d 
     449    
     450   SUBROUTINE lbc_lnk_2d_multiple( pt2d_array , type_array , psgn_array , num_fields ) 
     451      !! 
     452      INTEGER :: num_fields 
     453      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array 
     454      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   type_array   ! define the nature of ptab array grid-points 
     455      !                                                               ! = T , U , V , F , W and I points 
     456      REAL(wp)        , DIMENSION(:), INTENT(in   ) ::   psgn_array   ! =-1 the sign change across the north fold boundary 
     457      !                                                               ! =  1. , the sign is kept 
     458      ! 
     459      INTEGER  ::   ii    !!MULTI SEND DUMMY LOOP INDICES 
     460      ! 
     461      DO ii = 1, num_fields 
     462        CALL lbc_lnk_2d( pt2d_array(ii)%pt2d, type_array(ii), psgn_array(ii) ) 
     463      END DO      
     464      ! 
     465   END SUBROUTINE lbc_lnk_2d_multiple 
     466 
     467   SUBROUTINE lbc_lnk_2d_9( pt2dA, cd_typeA, psgnA, pt2dB, cd_typeB, psgnB, pt2dC, cd_typeC, psgnC   & 
     468      &                   , pt2dD, cd_typeD, psgnD, pt2dE, cd_typeE, psgnE, pt2dF, cd_typeF, psgnF   & 
     469      &                   , pt2dG, cd_typeG, psgnG, pt2dH, cd_typeH, psgnH, pt2dI, cd_typeI, psgnI, cd_mpp, pval) 
     470      !!--------------------------------------------------------------------- 
     471      ! Second 2D array on which the boundary condition is applied 
     472      REAL(wp), DIMENSION(jpi,jpj), TARGET          , INTENT(inout) ::   pt2dA 
     473      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dB , pt2dC , pt2dD , pt2dE 
     474      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dF , pt2dG , pt2dH , pt2dI 
     475      ! define the nature of ptab array grid-points 
     476      CHARACTER(len=1)                              , INTENT(in   ) ::   cd_typeA 
     477      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeB , cd_typeC , cd_typeD , cd_typeE 
     478      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeF , cd_typeG , cd_typeH , cd_typeI 
     479      ! =-1 the sign change across the north fold boundary 
     480      REAL(wp)                                      , INTENT(in   ) ::   psgnA 
     481      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnB , psgnC , psgnD , psgnE 
     482      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnF , psgnG , psgnH , psgnI 
     483      CHARACTER(len=3)                    , OPTIONAL, INTENT(in   ) ::   cd_mpp   ! fill the overlap area only 
     484      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   pval     ! background value (used at closed boundaries) 
     485      !! 
     486      !!--------------------------------------------------------------------- 
     487 
     488      !!The first array 
     489      CALL lbc_lnk( pt2dA, cd_typeA, psgnA )  
     490 
     491      !! Look if more arrays to process 
     492      IF(PRESENT (psgnB) )CALL lbc_lnk( pt2dA, cd_typeA, psgnA )  
     493      IF(PRESENT (psgnC) )CALL lbc_lnk( pt2dC, cd_typeC, psgnC )  
     494      IF(PRESENT (psgnD) )CALL lbc_lnk( pt2dD, cd_typeD, psgnD )  
     495      IF(PRESENT (psgnE) )CALL lbc_lnk( pt2dE, cd_typeE, psgnE )  
     496      IF(PRESENT (psgnF) )CALL lbc_lnk( pt2dF, cd_typeF, psgnF )  
     497      IF(PRESENT (psgnG) )CALL lbc_lnk( pt2dG, cd_typeG, psgnG )  
     498      IF(PRESENT (psgnH) )CALL lbc_lnk( pt2dH, cd_typeH, psgnH )  
     499      IF(PRESENT (psgnI) )CALL lbc_lnk( pt2dI, cd_typeI, psgnI )  
     500 
     501   END SUBROUTINE lbc_lnk_2d_9 
     502 
     503   SUBROUTINE lbc_lnk_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval ) 
     504      !!--------------------------------------------------------------------- 
     505      !!                 ***  ROUTINE lbc_lnk_sum_2d  *** 
     506      !! 
     507      !! ** Purpose :   set lateral boundary conditions on a 2D array (non mpp case) 
     508      !! 
     509      !! ** Comments:   compute the sum of the common cell (overlap region) for the ice sheet/ocean  
     510      !!                coupling if conservation option activated. As no ice shelf are present along 
     511      !!                this line, nothing is done along the north fold. 
     512      !!---------------------------------------------------------------------- 
     513      CHARACTER(len=1)            , INTENT(in   )           ::   cd_type   ! nature of pt3d grid-points 
     514      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   pt2d      ! 2D array on which the lbc is applied 
     515      REAL(wp)                    , INTENT(in   )           ::   psgn      ! control of the sign  
     516      CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   cd_mpp    ! MPP only (here do nothing) 
     517      REAL(wp)                    , INTENT(in   ), OPTIONAL ::   pval      ! background value (for closed boundaries) 
     518      !! 
     519      REAL(wp) ::   zland 
     520      !!---------------------------------------------------------------------- 
     521 
     522      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value (zero by default) 
     523      ELSE                         ;   zland = 0._wp 
     524      ENDIF 
     525 
     526      IF (PRESENT(cd_mpp)) THEN 
     527         ! only fill the overlap area and extra allows  
     528         ! this is in mpp case. In this module, just do nothing 
     529      ELSE 
     530         !                                     ! East-West boundaries 
     531         !                                     ! ==================== 
     532         SELECT CASE ( nperio ) 
     533         ! 
     534         CASE ( 1 , 4 , 6 )                       !** cyclic east-west 
     535            pt2d(jpim1,:) = pt2d(jpim1,:) + pt2d( 1 ,:) 
     536            pt2d(  2  ,:) = pt2d(  2  ,:) + pt2d(jpi,:) 
     537            pt2d( 1 ,:) = 0.0_wp               ! all points 
     538            pt2d(jpi,:) = 0.0_wp 
     539            ! 
     540         CASE DEFAULT                             !** East closed  --  West closed 
     541            SELECT CASE ( cd_type ) 
     542            CASE ( 'T' , 'U' , 'V' , 'W' )            ! T-, U-, V-, W-points 
     543               pt2d( 1 ,:) = zland 
     544               pt2d(jpi,:) = zland 
     545            CASE ( 'F' )                              ! F-point 
     546               pt2d(jpi,:) = zland 
     547            END SELECT 
     548            ! 
     549         END SELECT 
     550         !                                     ! North-South boundaries 
     551         !                                     ! ====================== 
     552         ! Nothing to do for the north fold, there is no ice shelf along this line. 
     553         ! 
     554      END IF 
     555 
     556   END SUBROUTINE 
     557 
     558   SUBROUTINE lbc_lnk_sum_3d( pt3d, cd_type, psgn, cd_mpp, pval ) 
     559      !!--------------------------------------------------------------------- 
     560      !!                 ***  ROUTINE lbc_lnk_sum_3d  *** 
     561      !! 
     562      !! ** Purpose :   set lateral boundary conditions on a 3D array (non mpp case) 
     563      !! 
     564      !! ** Comments:   compute the sum of the common cell (overlap region) for the ice sheet/ocean  
     565      !!                coupling if conservation option activated. As no ice shelf are present along 
     566      !!                this line, nothing is done along the north fold. 
     567      !!---------------------------------------------------------------------- 
     568      CHARACTER(len=1)                , INTENT(in   )           ::   cd_type   ! nature of pt3d grid-points 
     569      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   pt3d      ! 3D array on which the lbc is applied 
     570      REAL(wp)                        , INTENT(in   )           ::   psgn      ! control of the sign  
     571      CHARACTER(len=3)                , INTENT(in   ), OPTIONAL ::   cd_mpp    ! MPP only (here do nothing) 
     572      REAL(wp)                        , INTENT(in   ), OPTIONAL ::   pval      ! background value (for closed boundaries) 
     573      !! 
     574      REAL(wp) ::   zland 
     575      !!---------------------------------------------------------------------- 
     576 
     577      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value (zero by default) 
     578      ELSE                         ;   zland = 0._wp 
     579      ENDIF 
     580 
     581 
     582      IF( PRESENT( cd_mpp ) ) THEN 
     583         ! only fill the overlap area and extra allows  
     584         ! this is in mpp case. In this module, just do nothing 
     585      ELSE 
     586         !                                     !  East-West boundaries 
     587         !                                     ! ====================== 
     588         SELECT CASE ( nperio ) 
     589         ! 
     590         CASE ( 1 , 4 , 6 )                       !**  cyclic east-west 
     591            pt3d(jpim1,:,:) = pt3d(jpim1,:,:) + pt3d( 1 ,:,:) 
     592            pt3d(  2  ,:,:) = pt3d(  2  ,:,:) + pt3d(jpi,:,:)  
     593            pt3d( 1 ,:,:) = 0.0_wp            ! all points 
     594            pt3d(jpi,:,:) = 0.0_wp 
     595            ! 
     596         CASE DEFAULT                             !**  East closed  --  West closed 
     597            SELECT CASE ( cd_type ) 
     598            CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points 
     599               pt3d( 1 ,:,:) = zland 
     600               pt3d(jpi,:,:) = zland 
     601            CASE ( 'F' )                               ! F-point 
     602               pt3d(jpi,:,:) = zland 
     603            END SELECT 
     604            ! 
     605         END SELECT 
     606         !                                     ! North-South boundaries 
     607         !                                     ! ====================== 
     608         ! Nothing to do for the north fold, there is no ice shelf along this line. 
     609         ! 
     610      END IF 
     611   END SUBROUTINE 
     612 
    381613 
    382614#endif 
    … …  
    448680   !!====================================================================== 
    449681END MODULE lbclnk 
     682 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

    r6140 r7278  
    2424   !!            3.5  !  2013  ( C. Ethe, G. Madec ) message passing arrays as local variables  
    2525   !!            3.5  !  2013  (S.Mocavero, I.Epicoco - CMCC) north fold optimizations 
     26   !!            3.6  !  2015 (O. Tintó and M. Castrillo - BSC) Added 'mpp_lnk_2d_multiple', 'mpp_lbc_north_2d_multiple', 'mpp_max_multiple'  
    2627   !!---------------------------------------------------------------------- 
    2728 
    … …  
    6263   USE lbcnfd         ! north fold treatment 
    6364   USE in_out_manager ! I/O manager 
     65   USE wrk_nemo       ! work arrays 
    6466 
    6567   IMPLICIT NONE 
    … …  
    7072   PUBLIC   mpp_ini_north, mpp_lbc_north, mpp_lbc_north_e 
    7173   PUBLIC   mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc 
     74   PUBLIC   mpp_max_multiple 
    7275   PUBLIC   mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e 
    73    PUBLIC   mpp_lnk_2d_9  
     76   PUBLIC   mpp_lnk_2d_9 , mpp_lnk_2d_multiple  
    7477   PUBLIC   mpp_lnk_sum_3d, mpp_lnk_sum_2d 
    7578   PUBLIC   mppscatter, mppgather 
    … …  
    7982   PUBLIC   mpp_lnk_bdy_2d, mpp_lnk_bdy_3d 
    8083   PUBLIC   mpp_lbc_north_icb, mpp_lnk_2d_icb 
     84   PUBLIC   mpprank 
    8185 
    8286   TYPE arrayptr 
    8387      REAL , DIMENSION (:,:),  POINTER :: pt2d 
    8488   END TYPE arrayptr 
     89   PUBLIC   arrayptr 
    8590    
    8691   !! * Interfaces 
    … …  
    106111   INTERFACE mpp_maxloc 
    107112      MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d 
     113   END INTERFACE 
     114 
     115   INTERFACE mpp_max_multiple 
     116      MODULE PROCEDURE mppmax_real_multiple 
    108117   END INTERFACE 
    109118 
    … …  
    726735      ! ----------------------- 
    727736      ! 
    728       DO ii = 1 , num_fields 
    729737         !First Array 
    730          IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 
    731             ! 
    732             SELECT CASE ( jpni ) 
    733             CASE ( 1 )     ;   CALL lbc_nfd      ( pt2d_array(ii)%pt2d( : , : ), type_array(ii) , psgn_array(ii) )   ! only 1 northern proc, no mpp 
    734             CASE DEFAULT   ;   CALL mpp_lbc_north( pt2d_array(ii)%pt2d( : , : ), type_array(ii), psgn_array(ii) )   ! for all northern procs. 
    735             END SELECT 
    736             ! 
    737          ENDIF 
    738          ! 
    739       END DO 
     738      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN 
     739         ! 
     740         SELECT CASE ( jpni ) 
     741         CASE ( 1 )     ;    
     742             DO ii = 1 , num_fields   
     743                       CALL lbc_nfd      ( pt2d_array(ii)%pt2d( : , : ), type_array(ii) , psgn_array(ii) )   ! only 1 northern proc, no mpp 
     744             END DO 
     745         CASE DEFAULT   ;   CALL mpp_lbc_north_2d_multiple( pt2d_array, type_array, psgn_array, num_fields )   ! for all northern procs. 
     746         END SELECT 
     747         ! 
     748      ENDIF 
     749        ! 
    740750      ! 
    741751      DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we ) 
    … …  
    20192029   END SUBROUTINE mppmax_real 
    20202030 
     2031   SUBROUTINE mppmax_real_multiple( ptab, NUM , kcom  ) 
     2032      !!---------------------------------------------------------------------- 
     2033      !!                  ***  routine mppmax_real  *** 
     2034      !! 
     2035      !! ** Purpose :   Maximum 
     2036      !! 
     2037      !!---------------------------------------------------------------------- 
     2038      REAL(wp), DIMENSION(:) ,  INTENT(inout)           ::   ptab   ! ??? 
     2039      INTEGER , INTENT(in   )           ::   NUM 
     2040      INTEGER , INTENT(in   ), OPTIONAL ::   kcom   ! ??? 
     2041      !! 
     2042      INTEGER  ::   ierror, localcomm 
     2043      REAL(wp) , POINTER , DIMENSION(:) ::   zwork 
     2044      !!---------------------------------------------------------------------- 
     2045      ! 
     2046      CALL wrk_alloc(NUM , zwork) 
     2047      localcomm = mpi_comm_opa 
     2048      IF( PRESENT(kcom) )   localcomm = kcom 
     2049      ! 
     2050      CALL mpi_allreduce( ptab, zwork, NUM, mpi_double_precision, mpi_max, localcomm, ierror ) 
     2051      ptab = zwork 
     2052      CALL wrk_dealloc(NUM , zwork) 
     2053      ! 
     2054   END SUBROUTINE mppmax_real_multiple 
     2055 
    20212056 
    20222057   SUBROUTINE mppmin_a_real( ptab, kdim, kcom ) 
    … …  
    29122947   END SUBROUTINE mpp_lbc_north_2d 
    29132948 
     2949   SUBROUTINE mpp_lbc_north_2d_multiple( pt2d_array, cd_type, psgn, num_fields) 
     2950      !!--------------------------------------------------------------------- 
     2951      !!                   ***  routine mpp_lbc_north_2d  *** 
     2952      !! 
     2953      !! ** Purpose :   Ensure proper north fold horizontal bondary condition 
     2954      !!              in mpp configuration in case of jpn1 > 1 
     2955      !!              (for multiple 2d arrays ) 
     2956      !! 
     2957      !! ** Method  :   North fold condition and mpp with more than one proc 
     2958      !!              in i-direction require a specific treatment. We gather 
     2959      !!              the 4 northern lines of the global domain on 1 processor 
     2960      !!              and apply lbc north-fold on this sub array. Then we 
     2961      !!              scatter the north fold array back to the processors. 
     2962      !! 
     2963      !!---------------------------------------------------------------------- 
     2964      INTEGER ,  INTENT (in   ) ::   num_fields  ! number of variables contained in pt2d 
     2965      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array 
     2966      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   cd_type   ! nature of pt2d grid-points 
     2967      !                                                          !   = T ,  U , V , F or W  gridpoints 
     2968      REAL(wp), DIMENSION(:), INTENT(in   ) ::   psgn      ! = -1. the sign change across the north fold  
     2969      !!                                                             ! =  1. , the sign is kept 
     2970      INTEGER ::   ji, jj, jr, jk 
     2971      INTEGER ::   ierr, itaille, ildi, ilei, iilb 
     2972      INTEGER ::   ijpj, ijpjm1, ij, iproc 
     2973      INTEGER, DIMENSION (jpmaxngh)      ::   ml_req_nf          !for mpi_isend when avoiding mpi_allgather 
     2974      INTEGER                            ::   ml_err             ! for mpi_isend when avoiding mpi_allgather 
     2975      INTEGER, DIMENSION(MPI_STATUS_SIZE)::   ml_stat            ! for mpi_isend when avoiding mpi_allgather 
     2976      !                                                              ! Workspace for message transfers avoiding mpi_allgather 
     2977      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE   :: ztab 
     2978      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE   :: znorthloc, zfoldwk 
     2979      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE   :: znorthgloio 
     2980      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE   :: ztabl, ztabr 
     2981      INTEGER :: istatus(mpi_status_size) 
     2982      INTEGER :: iflag 
     2983      !!---------------------------------------------------------------------- 
     2984      ! 
     2985      ALLOCATE( ztab(jpiglo,4,num_fields), znorthloc(jpi,4,num_fields), zfoldwk(jpi,4,num_fields), znorthgloio(jpi,4,num_fields,jpni) )   ! expanded to 3 dimensions 
     2986      ALLOCATE( ztabl(jpi,4,num_fields), ztabr(jpi*jpmaxngh, 4,num_fields) ) 
     2987      ! 
     2988      ijpj   = 4 
     2989      ijpjm1 = 3 
     2990      ! 
     2991       
     2992      DO jk = 1, num_fields 
     2993         DO jj = nlcj-ijpj+1, nlcj             ! put in znorthloc the last 4 jlines of pt2d (for every variable) 
     2994            ij = jj - nlcj + ijpj 
     2995            znorthloc(:,ij,jk) = pt2d_array(jk)%pt2d(:,jj) 
     2996         END DO 
     2997      END DO 
     2998      !                                     ! Build in procs of ncomm_north the znorthgloio 
     2999      itaille = jpi * ijpj 
     3000                                                                   
     3001      IF ( l_north_nogather ) THEN 
     3002         ! 
     3003         ! Avoid the use of mpi_allgather by exchanging only with the processes already identified  
     3004         ! (in nemo_northcomms) as being  involved in this process' northern boundary exchange 
     3005         ! 
     3006         ztabr(:,:,:) = 0 
     3007         ztabl(:,:,:) = 0 
     3008 
     3009         DO jk = 1, num_fields 
     3010            DO jj = nlcj-ijpj+1, nlcj          ! First put local values into the global array 
     3011               ij = jj - nlcj + ijpj 
     3012               DO ji = nfsloop, nfeloop 
     3013                  ztabl(ji,ij,jk) = pt2d_array(jk)%pt2d(ji,jj) 
     3014               END DO 
     3015            END DO 
     3016         END DO 
     3017 
     3018         DO jr = 1,nsndto 
     3019            IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN 
     3020               CALL mppsend(5, znorthloc, itaille*num_fields, nfipproc(isendto(jr),jpnj), ml_req_nf(jr)) ! Buffer expanded "num_fields" times 
     3021            ENDIF 
     3022         END DO 
     3023         DO jr = 1,nsndto 
     3024            iproc = nfipproc(isendto(jr),jpnj) 
     3025            IF(iproc .ne. -1) THEN 
     3026               ilei = nleit (iproc+1) 
     3027               ildi = nldit (iproc+1) 
     3028               iilb = nfiimpp(isendto(jr),jpnj) - nfiimpp(isendto(1),jpnj) 
     3029            ENDIF 
     3030            IF((iproc .ne. (narea-1)) .and. (iproc .ne. -1)) THEN 
     3031              CALL mpprecv(5, zfoldwk, itaille*num_fields, iproc) ! Buffer expanded "num_fields" times 
     3032              DO jk = 1 , num_fields 
     3033                 DO jj = 1, ijpj 
     3034                    DO ji = ildi, ilei 
     3035                       ztabr(iilb+ji,jj,jk) = zfoldwk(ji,jj,jk)       ! Modified to 3D 
     3036                    END DO 
     3037                 END DO 
     3038              END DO 
     3039            ELSE IF (iproc .eq. (narea-1)) THEN 
     3040              DO jk = 1, num_fields 
     3041                 DO jj = 1, ijpj 
     3042                    DO ji = ildi, ilei 
     3043                          ztabr(iilb+ji,jj,jk) = pt2d_array(jk)%pt2d(ji,nlcj-ijpj+jj)       ! Modified to 3D 
     3044                    END DO 
     3045                 END DO 
     3046              END DO 
     3047            ENDIF 
     3048         END DO 
     3049         IF (l_isend) THEN 
     3050            DO jr = 1,nsndto 
     3051               IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN 
     3052                  CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err) 
     3053               ENDIF 
     3054            END DO 
     3055         ENDIF 
     3056         ! 
     3057         DO ji = 1, num_fields     ! Loop to manage 3D variables 
     3058            CALL mpp_lbc_nfd( ztabl(:,:,ji), ztabr(:,:,ji), cd_type(ji), psgn(ji) )  ! North fold boundary condition 
     3059         END DO 
     3060         ! 
     3061         DO jk = 1, num_fields 
     3062            DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt2d 
     3063               ij = jj - nlcj + ijpj 
     3064               DO ji = 1, nlci 
     3065                  pt2d_array(jk)%pt2d(ji,jj) = ztabl(ji,ij,jk)       ! Modified to 3D 
     3066               END DO 
     3067            END DO 
     3068         END DO 
     3069          
     3070         ! 
     3071      ELSE 
     3072         ! 
     3073         CALL MPI_ALLGATHER( znorthloc  , itaille*num_fields, MPI_DOUBLE_PRECISION,        & 
     3074            &                znorthgloio, itaille*num_fields, MPI_DOUBLE_PRECISION, ncomm_north, ierr ) 
     3075         ! 
     3076         ztab(:,:,:) = 0.e0 
     3077         DO jk = 1, num_fields 
     3078            DO jr = 1, ndim_rank_north            ! recover the global north array 
     3079               iproc = nrank_north(jr) + 1 
     3080               ildi = nldit (iproc) 
     3081               ilei = nleit (iproc) 
     3082               iilb = nimppt(iproc) 
     3083               DO jj = 1, ijpj 
     3084                  DO ji = ildi, ilei 
     3085                     ztab(ji+iilb-1,jj,jk) = znorthgloio(ji,jj,jk,jr) 
     3086                  END DO 
     3087               END DO 
     3088            END DO 
     3089         END DO 
     3090          
     3091         DO ji = 1, num_fields 
     3092            CALL lbc_nfd( ztab(:,:,ji), cd_type(ji), psgn(ji) )   ! North fold boundary condition 
     3093         END DO 
     3094         ! 
     3095         DO jk = 1, num_fields 
     3096            DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt2d 
     3097               ij = jj - nlcj + ijpj 
     3098               DO ji = 1, nlci 
     3099                  pt2d_array(jk)%pt2d(ji,jj) = ztab(ji+nimpp-1,ij,jk) 
     3100               END DO 
     3101            END DO 
     3102         END DO 
     3103         ! 
     3104         ! 
     3105      ENDIF 
     3106      DEALLOCATE( ztab, znorthloc, zfoldwk, znorthgloio ) 
     3107      DEALLOCATE( ztabl, ztabr ) 
     3108      ! 
     3109   END SUBROUTINE mpp_lbc_north_2d_multiple 
    29143110 
    29153111   SUBROUTINE mpp_lbc_north_e( pt2d, cd_type, psgn) 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90

    r6140 r7278  
    198198       
    199199#endif 
    200       IF(lwp) THEN 
    201          WRITE(numout,*) 
    202          WRITE(numout,*) '           defines mpp subdomains' 
    203          WRITE(numout,*) '           ----------------------' 
    204          WRITE(numout,*) '           iresti=',iresti,' irestj=',irestj 
    205          WRITE(numout,*) '           jpni  =',jpni  ,' jpnj  =',jpnj 
    206          ifreq = 4 
    207          il1   = 1 
    208          DO jn = 1, (jpni-1)/ifreq+1 
    209             il2 = MIN( jpni, il1+ifreq-1 ) 
    210             WRITE(numout,*) 
    211             WRITE(numout,9200) ('***',ji = il1,il2-1) 
    212             DO jj = jpnj, 1, -1 
    213                WRITE(numout,9203) ('   ',ji = il1,il2-1) 
    214                WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 ) 
    215                WRITE(numout,9203) ('   ',ji = il1,il2-1) 
    216                WRITE(numout,9200) ('***',ji = il1,il2-1) 
    217             END DO 
    218             WRITE(numout,9201) (ji,ji = il1,il2) 
    219             il1 = il1+ifreq 
    220          END DO 
    221  9200    FORMAT('     ***',20('*************',a3)) 
    222  9203    FORMAT('     *     ',20('         *   ',a3)) 
    223  9201    FORMAT('        ',20('   ',i3,'          ')) 
    224  9202    FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   ')) 
    225       ENDIF 
    226  
    227       zidom = nreci 
    228       DO ji = 1, jpni 
    229          zidom = zidom + ilcit(ji,1) - nreci 
    230       END DO 
    231       IF(lwp) WRITE(numout,*) 
    232       IF(lwp) WRITE(numout,*)' sum ilcit(i,1) = ', zidom, ' jpiglo = ', jpiglo 
    233        
    234       zjdom = nrecj 
    235       DO jj = 1, jpnj 
    236          zjdom = zjdom + ilcjt(1,jj) - nrecj 
    237       END DO 
    238       IF(lwp) WRITE(numout,*)' sum ilcit(1,j) = ', zjdom, ' jpjglo = ', jpjglo 
    239       IF(lwp) WRITE(numout,*) 
    240        
    241200 
    242201      !  2. Index arrays for subdomains 
    … …  
    301260         nlejt(jn) = nlej 
    302261      END DO 
    303        
    304  
    305       ! 4. From global to local 
     262 
     263      ! 4. Subdomain print 
     264      ! ------------------ 
     265       
     266      IF(lwp) WRITE(numout,*) 
     267      IF(lwp) WRITE(numout,*) ' mpp_init: defines mpp subdomains' 
     268      IF(lwp) WRITE(numout,*) ' ~~~~~~  ----------------------' 
     269      IF(lwp) WRITE(numout,*) 
     270      IF(lwp) WRITE(numout,*) 'iresti=',iresti,' irestj=',irestj 
     271      IF(lwp) WRITE(numout,*) 
     272      IF(lwp) WRITE(numout,*) 'jpni=',jpni,' jpnj=',jpnj 
     273      zidom = nreci 
     274      DO ji = 1, jpni 
     275         zidom = zidom + ilcit(ji,1) - nreci 
     276      END DO 
     277      IF(lwp) WRITE(numout,*) 
     278      IF(lwp) WRITE(numout,*)' sum ilcit(i,1)=', zidom, ' jpiglo=', jpiglo 
     279 
     280      zjdom = nrecj 
     281      DO jj = 1, jpnj 
     282         zjdom = zjdom + ilcjt(1,jj) - nrecj 
     283      END DO 
     284      IF(lwp) WRITE(numout,*)' sum ilcit(1,j)=', zjdom, ' jpjglo=', jpjglo 
     285      IF(lwp) WRITE(numout,*) 
     286 
     287      IF(lwp) THEN 
     288         ifreq = 4 
     289         il1   = 1 
     290         DO jn = 1, (jpni-1)/ifreq+1 
     291            il2 = MIN( jpni, il1+ifreq-1 ) 
     292            WRITE(numout,*) 
     293            WRITE(numout,9200) ('***',ji = il1,il2-1) 
     294            DO jj = jpnj, 1, -1 
     295               WRITE(numout,9203) ('   ',ji = il1,il2-1) 
     296               WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 ) 
     297               WRITE(numout,9204) (nfipproc(ji,jj),ji=il1,il2) 
     298               WRITE(numout,9203) ('   ',ji = il1,il2-1) 
     299               WRITE(numout,9200) ('***',ji = il1,il2-1) 
     300            END DO 
     301            WRITE(numout,9201) (ji,ji = il1,il2) 
     302            il1 = il1+ifreq 
     303         END DO 
     304 9200     FORMAT('     ***',20('*************',a3)) 
     305 9203     FORMAT('     *     ',20('         *   ',a3)) 
     306 9201     FORMAT('        ',20('   ',i3,'          ')) 
     307 9202     FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   ')) 
     308 9204     FORMAT('     *  ',20('      ',i3,'   *   ')) 
     309      ENDIF 
     310 
     311      ! 5. From global to local 
    306312      ! ----------------------- 
    307313 
    … …  
    310316 
    311317 
    312       ! 5. Subdomain neighbours 
     318      ! 6. Subdomain neighbours 
    313319      ! ---------------------- 
    314320 
    … …  
    433439         WRITE(numout,*) ' nimpp  = ', nimpp 
    434440         WRITE(numout,*) ' njmpp  = ', njmpp 
    435          WRITE(numout,*) ' nbse   = ', nbse  , ' npse   = ', npse 
    436          WRITE(numout,*) ' nbsw   = ', nbsw  , ' npsw   = ', npsw 
    437          WRITE(numout,*) ' nbne   = ', nbne  , ' npne   = ', npne 
    438          WRITE(numout,*) ' nbnw   = ', nbnw  , ' npnw   = ', npnw 
     441         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse 
     442         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw 
     443         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne 
     444         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw 
     445         WRITE(numout,*) 
    439446      ENDIF 
    440447 
    … …  
    443450      ! Prepare mpp north fold 
    444451 
    445       IF (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
     452      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
    446453         CALL mpp_ini_north 
    447       END IF 
     454         IF(lwp) WRITE(numout,*) ' mpp_init : North fold boundary prepared for jpni >1' 
     455      ENDIF 
    448456 
    449457      ! Prepare NetCDF output file (if necessary) 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90

    r7277 r7278  
    276276         ENDIF 
    277277 
     278         ! Check wet points over the entire domain to preserve the MPI communication stencil 
    278279         isurf = 0 
    279          DO jj = 1+jprecj, ilj-jprecj 
    280             DO  ji = 1+jpreci, ili-jpreci 
     280         DO jj = 1, ilj 
     281            DO  ji = 1, ili 
    281282               IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1)   isurf = isurf+1 
    282283            END DO 
    283284         END DO 
     285 
    284286         IF(isurf /= 0) THEN 
    285287            icont = icont + 1 
    … …  
    291293 
    292294      nfipproc(:,:) = ipproc(:,:) 
    293  
    294295 
    295296      ! Control 
    … …  
    399400      ii = iin(narea) 
    400401      ij = ijn(narea) 
     402 
     403      ! set default neighbours 
     404      noso = ioso(ii,ij) 
     405      nowe = iowe(ii,ij) 
     406      noea = ioea(ii,ij) 
     407      nono = iono(ii,ij)  
     408      npse = iose(ii,ij) 
     409      npsw = iosw(ii,ij) 
     410      npne = ione(ii,ij) 
     411      npnw = ionw(ii,ij) 
     412 
     413      ! check neighbours location 
    401414      IF( ioso(ii,ij) >= 0 .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN  
    402415         iiso = 1 + MOD(ioso(ii,ij),jpni) 
    … …  
    469482      IF (lwp) THEN 
    470483         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea ) 
     484         WRITE(inum,'(a)') '   jpnij     jpi     jpj     jpk  jpiglo  jpjglo' 
    471485         WRITE(inum,'(6i8)') jpnij,jpi,jpj,jpk,jpiglo,jpjglo 
    472486         WRITE(inum,'(a)') 'NAREA nlci nlcj nldi nldj nlei nlej nimpp njmpp' 
    … …  
    481495      END IF 
    482496 
    483       IF( nperio == 1 .AND.jpni /= 1 ) CALL ctl_stop( ' mpp_init2:  error on cyclicity' ) 
    484  
    485       ! Prepare mpp north fold 
    486  
    487       IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
    488          CALL mpp_ini_north 
    489          IF(lwp) WRITE(numout,*) ' mpp_init2 : North fold boundary prepared for jpni >1' 
    490       ENDIF 
    491  
    492497      ! Defined npolj, either 0, 3 , 4 , 5 , 6 
    493498      ! In this case the important thing is that npolj /= 0 
    … …  
    506511      ENDIF 
    507512 
     513      ! Periodicity : no corner if nbondi = 2 and nperio != 1 
     514 
     515      IF(lwp) THEN 
     516         WRITE(numout,*) ' nproc  = ', nproc 
     517         WRITE(numout,*) ' nowe   = ', nowe  , ' noea   =  ', noea 
     518         WRITE(numout,*) ' nono   = ', nono  , ' noso   =  ', noso 
     519         WRITE(numout,*) ' nbondi = ', nbondi 
     520         WRITE(numout,*) ' nbondj = ', nbondj 
     521         WRITE(numout,*) ' npolj  = ', npolj 
     522         WRITE(numout,*) ' nperio = ', nperio 
     523         WRITE(numout,*) ' nlci   = ', nlci 
     524         WRITE(numout,*) ' nlcj   = ', nlcj 
     525         WRITE(numout,*) ' nimpp  = ', nimpp 
     526         WRITE(numout,*) ' njmpp  = ', njmpp 
     527         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse 
     528         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw 
     529         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne 
     530         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw 
     531         WRITE(numout,*) 
     532      ENDIF 
     533 
     534      IF( nperio == 1 .AND. jpni /= 1 ) CALL ctl_stop( ' mpp_init2: error on cyclicity' ) 
     535 
     536      ! Prepare mpp north fold 
     537 
     538      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
     539         CALL mpp_ini_north 
     540         IF(lwp) WRITE(numout,*) ' mpp_init2 : North fold boundary prepared for jpni >1' 
     541      ENDIF 
     542 
    508543      ! Prepare NetCDF output file (if necessary) 
    509544      CALL mpp_init_ioipsl 
    510545 
    511       ! Periodicity : no corner if nbondi = 2 and nperio != 1 
    512  
    513       IF(lwp) THEN 
    514          WRITE(numout,*) ' nproc=  ',nproc 
    515          WRITE(numout,*) ' nowe=   ',nowe 
    516          WRITE(numout,*) ' noea=   ',noea 
    517          WRITE(numout,*) ' nono=   ',nono 
    518          WRITE(numout,*) ' noso=   ',noso 
    519          WRITE(numout,*) ' nbondi= ',nbondi 
    520          WRITE(numout,*) ' nbondj= ',nbondj 
    521          WRITE(numout,*) ' npolj=  ',npolj 
    522          WRITE(numout,*) ' nperio= ',nperio 
    523          WRITE(numout,*) ' nlci=   ',nlci 
    524          WRITE(numout,*) ' nlcj=   ',nlcj 
    525          WRITE(numout,*) ' nimpp=  ',nimpp 
    526          WRITE(numout,*) ' njmpp=  ',njmpp 
    527          WRITE(numout,*) ' nbse=   ',nbse,' npse= ',npse 
    528          WRITE(numout,*) ' nbsw=   ',nbsw,' npsw= ',npsw 
    529          WRITE(numout,*) ' nbne=   ',nbne,' npne= ',npne 
    530          WRITE(numout,*) ' nbnw=   ',nbnw,' npnw= ',npnw 
    531       ENDIF 
    532546 
    533547   END SUBROUTINE mpp_init2 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r4624 r7278  
    99   !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0 
    1010   !!             -   ! 2006-08  (G. Madec)  cleaning for surface module 
     11   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo 
    1112   !!---------------------------------------------------------------------- 
    1213 
    … …  
    2930 
    3031   INTEGER  ::   albd_init = 0      !: control flag for initialization 
    31    REAL(wp) ::   zzero     = 0.e0   ! constant values 
    32    REAL(wp) ::   zone      = 1.e0   !    "       " 
    33  
    34    REAL(wp) ::   c1     = 0.05    ! constants values 
    35    REAL(wp) ::   c2     = 0.10    !    "        " 
    36    REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude 
    37  
     32   
     33   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude 
     34   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     35   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0) 
     36   REAL(wp) ::   c2       = 0.10    !  "        " 
     37   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
     38  
    3839   !                             !!* namelist namsbc_alb 
    39    REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984) 
    40 #if defined key_lim3 
    41    REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) 
    42 #else 
    43    REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) 
    44 #endif 
    45    REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute 
    46    REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972) 
    47    REAL(wp) ::   rn_alphc         !  
     40   INTEGER  ::   nn_ice_alb 
     41   REAL(wp) ::   rn_albice 
    4842 
    4943   !!---------------------------------------------------------------------- 
    … …  
    5953      !!           
    6054      !! ** Purpose :   Computation of the albedo of the snow/ice system  
    61       !!                as well as the ocean one 
    6255      !!        
    63       !! ** Method  : - Computation of the albedo of snow or ice (choose the  
    64       !!                rignt one by a large number of tests 
    65       !!              - Computation of the albedo of the ocean 
    66       !! 
    67       !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250. 
     56      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb) 
     57      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies 
     58      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005) 
     59      !!                                                                           and Grenfell & Perovich (JGR 2004) 
     60      !!                Description of scheme 1: 
     61      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005) 
     62      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999 
     63      !!                     0-5cm  : linear function of ice thickness 
     64      !!                     5-150cm: log    function of ice thickness 
     65      !!                     > 150cm: constant 
     66      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004) 
     67      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting 
     68      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004) 
     69      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law 
     70      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice 
     71      !! 
     72      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions: 
     73      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo 
     74      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger  
     75      !!                     under melting conditions than under freezing conditions 
     76      !!                  3) the evolution of ice albedo as a function of ice thickness shows   
     77      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 
     78      !! 
     79      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 
     80      !!                Brandt et al. 2005, J. Climate, vol 18 
     81      !!                Grenfell & Perovich 2004, JGR, vol 109  
    6882      !!---------------------------------------------------------------------- 
    6983      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin) 
    … …  
    7387      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky 
    7488      !! 
    75       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    76       INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays) 
    77       REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting 
    78       REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing 
    79       REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow 
    80       REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow 
    81       REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow) 
    82       REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow) 
    83       REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1 
    84       REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2 
    85       !! 
    86       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice 
    87       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness 
     89      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
     90      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays) 
     91      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if 
     92      REAL(wp)            ::   zswitch, z1_c1, z1_c2 
     93      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total 
     94      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free) 
    8895      !!--------------------------------------------------------------------- 
    89        
     96 
    9097      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories 
    91  
    92       CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth ) 
     98       
     99      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    93100 
    94101      IF( albd_init == 0 )   CALL albedo_init      ! initialization  
    95102 
    96       !--------------------------- 
    97       !  Computation of  zficeth 
    98       !--------------------------- 
    99       ! ice free of snow and melts 
    100       WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice 
    101       ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi 
    102       END  WHERE 
    103  
    104       WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz 
    105       ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 ) 
    106       ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              & 
    107          &                                                                 - 0.8608 * ph_ice * ph_ice     & 
    108          &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
    109       ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice 
    110       END WHERE 
    111  
    112 !!gm old code 
    113 !      DO jl = 1, ijpl 
    114 !         DO jj = 1, jpj 
    115 !            DO ji = 1, jpi 
    116 !               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN 
    117 !                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl) 
    118 !               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN 
    119 !                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 ) 
    120 !               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN 
    121 !                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               & 
    122 !                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 & 
    123 !                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl) 
    124 !               ELSE 
    125 !                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl)  
    126 !               ENDIF 
    127 !            END DO 
    128 !         END DO 
    129 !      END DO 
    130 !!gm end old code 
    131        
    132       !-----------------------------------------------  
    133       !    Computation of the snow/ice albedo system  
    134       !-------------------------- --------------------- 
    135        
    136       !    Albedo of snow-ice for clear sky. 
    137       !-----------------------------------------------     
    138       DO jl = 1, ijpl 
    139          DO jj = 1, jpj 
    140             DO ji = 1, jpi 
    141                !  Case of ice covered by snow.              
    142                !                                        !  freezing snow         
    143                zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
    144                zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             & 
    145                   &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   & 
    146                   &     +         zihsc1   * rn_alphd   
    147                !                                        !  melting snow                 
    148                zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) ) 
    149                zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   & 
    150                   &     +         zihsc2   *   rn_alphc  
    151                ! 
    152                zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )    
    153                zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf 
    154              
    155                !  Case of ice free of snow. 
    156                zalbpic  = zficeth(ji,jj,jl)  
    157              
    158                ! albedo of the system    
    159                zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) ) 
    160                pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic 
     103       
     104      SELECT CASE ( nn_ice_alb ) 
     105 
     106      !------------------------------------------ 
     107      !  Shine and Henderson-Sellers (1985) 
     108      !------------------------------------------ 
     109      CASE( 0 ) 
     110        
     111         ralb_sf = 0.80       ! dry snow 
     112         ralb_sm = 0.65       ! melting snow 
     113         ralb_if = 0.72       ! bare frozen ice 
     114         ralb_im = rn_albice  ! bare puddled ice  
     115          
     116         !  Computation of ice albedo (free of snow) 
     117         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im 
     118         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if 
     119         END  WHERE 
     120       
     121         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
     122         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
     123         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              & 
     124            &                                                                 - 0.8608 * ph_ice * ph_ice     & 
     125            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
     126         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
     127         END WHERE 
     128      
     129         DO jl = 1, ijpl 
     130            DO jj = 1, jpj 
     131               DO ji = 1, jpi 
     132                  ! freezing snow 
     133                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
     134                  !                                        !  freezing snow         
     135                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
     136                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
     137                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   & 
     138                     &        +         zswitch   * ralb_sf   
     139 
     140                  ! melting snow 
     141                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2 
     142                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
     143                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   & 
     144                      &     +         zswitch   *   ralb_sm  
     145                  ! 
     146                  ! snow albedo 
     147                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     148                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
     149                
     150                  ! Ice/snow albedo 
     151                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     152                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl) 
     153                  ! 
     154               END DO 
    161155            END DO 
    162156         END DO 
    163       END DO 
    164        
    165       !    Albedo of snow-ice for overcast sky. 
    166       !----------------------------------------------   
    167       pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction 
    168       ! 
    169       CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth ) 
     157 
     158         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky 
     159 
     160      !------------------------------------------ 
     161      !  New parameterization (2016) 
     162      !------------------------------------------ 
     163      CASE( 1 )  
     164 
     165         ralb_im = rn_albice  ! bare puddled ice 
     166! compilation of values from literature 
     167         ralb_sf = 0.85      ! dry snow 
     168         ralb_sm = 0.75      ! melting snow 
     169         ralb_if = 0.60      ! bare frozen ice 
     170! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
     171!         ralb_sf = 0.85       ! dry snow 
     172!         ralb_sm = 0.72       ! melting snow 
     173!         ralb_if = 0.65       ! bare frozen ice 
     174! Brandt et al 2005 (East Antarctica) 
     175!         ralb_sf = 0.87      ! dry snow 
     176!         ralb_sm = 0.82      ! melting snow 
     177!         ralb_if = 0.54      ! bare frozen ice 
     178!  
     179         !  Computation of ice albedo (free of snow) 
     180         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
     181         z1_c2 = 1. / 0.05 
     182         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im 
     183         ELSE WHERE                                              ;   zalb = ralb_if 
     184         END  WHERE 
     185          
     186         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
     187         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
     188            &                                                                     ( LOG(1.5) - LOG(ph_ice) ) 
     189         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice 
     190         END WHERE 
     191 
     192         z1_c1 = 1. / 0.02 
     193         z1_c2 = 1. / 0.03 
     194         !  Computation of the snow/ice albedo 
     195         DO jl = 1, ijpl 
     196            DO jj = 1, jpj 
     197               DO ji = 1, jpi 
     198                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 ); 
     199                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 ); 
     200 
     201                   ! snow albedo 
     202                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     203                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
     204 
     205                  ! Ice/snow albedo    
     206                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     207                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
     208 
     209              END DO 
     210            END DO 
     211         END DO 
     212         ! Effect of the clouds (2d order polynomial) 
     213         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );  
     214 
     215      END SELECT 
     216       
     217      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    170218      ! 
    171219   END SUBROUTINE albedo_ice 
    … …  
    181229      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky 
    182230      !! 
    183       REAL(wp) ::   zcoef   ! local scalar 
    184       !!---------------------------------------------------------------------- 
    185       ! 
    186       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982  
    187       pa_oce_cs(:,:) = zcoef                
    188       pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972 
     231      REAL(wp) :: zcoef  
     232      !!---------------------------------------------------------------------- 
     233      ! 
     234      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982 
     235      pa_oce_cs(:,:) = zcoef  
     236      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972 
    189237      ! 
    190238   END SUBROUTINE albedo_oce 
    … …  
    200248      !!---------------------------------------------------------------------- 
    201249      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    202       NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc 
     250      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice  
    203251      !!---------------------------------------------------------------------- 
    204252      ! 
    … …  
    219267         WRITE(numout,*) '~~~~~~~' 
    220268         WRITE(numout,*) '   Namelist namsbc_alb : albedo ' 
    221          WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud 
    222          WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice 
    223          WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd 
    224          WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi 
    225          WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc 
     269         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb 
     270         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice 
    226271      ENDIF 
    227272      ! 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r5407 r7278  
    8080   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_oce       !: heat flux of precip and evap over ocean     [W/m2] 
    8181   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_ice       !: heat flux of precip and evap over ice       [W/m2] 
    82    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: heat flux of precip over ice                [J/m3] 
     82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qevap_ice      !: heat flux of evap over ice                  [W/m2] 
     83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: enthalpy of precip over ice                 [J/m3] 
    8384   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s] 
    8485#endif 
    … …  
    144145#endif 
    145146#if defined key_lim3 
    146          &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,  & 
    147          &      qemp_ice(jpi,jpj)     , qemp_oce(jpi,jpj)      ,                       & 
    148          &      qns_oce (jpi,jpj)     , qsr_oce (jpi,jpj)      , emp_oce (jpi,jpj)  ,  & 
     147         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   & 
     148         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   & 
     149         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce (jpi,jpj)  ,   & 
    149150#endif 
    150151         &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) ) 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r7277 r7278  
    668668      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
    669669 
     670      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
     671      DO jl = 1, jpl 
     672         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     673                                   ! but then qemp_ice should also include sublimation  
     674      END DO 
     675 
    670676      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
    671677#endif 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r7277 r7278  
    612612      ! --- evaporation --- ! 
    613613      z1_lsub = 1._wp / Lsub 
    614       evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation 
    615       devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub 
    616       zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean 
     614      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation 
     615      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT 
     616      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean 
    617617 
    618618      ! --- evaporation minus precipitation --- ! 
    … …  
    637637      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    638638      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     639 
     640      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
     641      DO jl = 1, jpl 
     642         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
     643                                   ! But we do not have Tice => consider it at 0°C => evap=0  
     644      END DO 
    639645 
    640646      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7277 r7278  
    10061006      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling 
    10071007         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1) 
    1008          IF( srcv(jpr_soce)%laction .AND. ln_useCT ) THEN    ! make sure that sst_m is the potential temperature 
     1008         IF( srcv(jpr_soce)%laction .AND. l_useCT ) THEN    ! make sure that sst_m is the potential temperature 
    10091009            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) ) 
    10101010         ENDIF 
    … …  
    13701370      ! 
    13711371      INTEGER ::   jl   ! dummy loop index 
    1372       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk 
    1373       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot 
    1374       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice 
    1375       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3 
     1372      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw 
     1373      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice 
     1374      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
     1375      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 
    13761376      !!---------------------------------------------------------------------- 
    13771377      ! 
    13781378      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_ice_flx') 
    13791379      ! 
    1380       CALL wrk_alloc( jpi,jpj,       zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot ) 
    1381       CALL wrk_alloc( jpi,jpj,jpl,   zqns_ice, zqsr_ice, zdqns_ice ) 
     1380      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw ) 
     1381      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice ) 
     1382      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
     1383      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
    13821384 
    13831385      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
    … …  
    14141416         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:) 
    14151417      END SELECT 
    1416  
    1417       IF( iom_use('subl_ai_cea') )   & 
    1418          CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average) 
    1419       !    
    1420       !                                                           ! runoffs and calving (put in emp_tot) 
     1418#if defined key_lim3 
     1419      ! zsnw = snow percentage over ice after wind blowing 
     1420      zsnw(:,:) = 0._wp 
     1421      CALL lim_thd_snwblow( p_frld, zsnw ) 
     1422       
     1423      ! --- evaporation (used later in sbccpl) --- ! 
     1424      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) 
     1425 
     1426      ! --- evaporation over ice (kg/m2/s) --- ! 
     1427      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) 
     1428      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
     1429      ! therefore, sublimation is not redistributed over the ice categories in case no subgrid scale fluxes are provided by atm. 
     1430      zdevap_ice(:,:) = 0._wp 
     1431       
     1432      ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- ! 
     1433      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw) 
     1434      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)     
     1435 
     1436      ! --- runoffs (included in emp later on) --- ! 
     1437      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
     1438 
     1439      ! --- calving (put in emp_tot and emp_oce) --- ! 
     1440      IF( srcv(jpr_cal)%laction ) THEN  
     1441         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) 
     1442         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1) 
     1443         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) ) 
     1444      ENDIF 
     1445       
     1446      IF( ln_mixcpl ) THEN 
     1447         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) 
     1448         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:) 
     1449         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:) 
     1450         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:) 
     1451         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:) 
     1452         DO jl=1,jpl 
     1453            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:) 
     1454            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:) 
     1455         ENDDO 
     1456      ELSE 
     1457         emp_tot(:,:) =         zemp_tot(:,:) 
     1458         emp_ice(:,:) =         zemp_ice(:,:) 
     1459         emp_oce(:,:) =         zemp_oce(:,:)      
     1460         sprecip(:,:) =         zsprecip(:,:) 
     1461         tprecip(:,:) =         ztprecip(:,:) 
     1462         DO jl=1,jpl 
     1463            evap_ice (:,:,jl) = zevap_ice (:,:) 
     1464            devap_ice(:,:,jl) = zdevap_ice(:,:) 
     1465         ENDDO 
     1466      ENDIF 
     1467       
     1468      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)    )  ! Sublimation over sea-ice (cell average) 
     1469                                     CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow 
     1470      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average) 
     1471      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)     
     1472#else 
     1473      ! Sublimation over sea-ice (cell average) 
     1474      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) 
     1475      ! runoffs and calving (put in emp_tot) 
    14211476      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
    14221477      IF( srcv(jpr_cal)%laction ) THEN  
    … …  
    14421497      IF( iom_use('snow_ai_cea') )   & 
    14431498         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average) 
     1499#endif 
    14441500 
    14451501      !                                                      ! ========================= ! 
    … …  
    14971553      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
    14981554 
    1499 #if defined key_lim3 
    1500       CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce )  
    1501  
    1502       ! --- evaporation --- ! 
    1503       ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation 
    1504       ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice 
    1505       !                 but it is incoherent WITH the ice model   
    1506       DO jl=1,jpl 
    1507          evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1) 
    1508       ENDDO 
    1509       zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean 
    1510  
    1511       ! --- evaporation minus precipitation --- ! 
    1512       emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:) 
    1513  
     1555#if defined key_lim3       
    15141556      ! --- non solar flux over ocean --- ! 
    15151557      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax 
    … …  
    15171559      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:) 
    15181560 
    1519       ! --- heat flux associated with emp --- ! 
    1520       zsnw(:,:) = 0._wp 
    1521       CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing 
     1561      ! --- heat flux associated with emp (W/m2) --- ! 
    15221562      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap 
    15231563         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip 
    15241564         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean 
    1525       qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap 
    1526          &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice 
    1527  
     1565!      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap 
     1566!         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice 
     1567      zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only) 
     1568                                                                                                       ! qevap_ice=0 since we consider Tice=0°C 
     1569       
    15281570      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    15291571      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus ) 
    15301572 
    1531       ! --- total non solar flux --- ! 
    1532       zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:) 
     1573      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
     1574      DO jl = 1, jpl 
     1575         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0°C 
     1576      END DO 
     1577 
     1578      ! --- total non solar flux (including evap/precip) --- ! 
     1579      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:) 
    15331580 
    15341581      ! --- in case both coupled/forced are active, we must mix values --- !  
    … …  
    15371584         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:) 
    15381585         DO jl=1,jpl 
    1539             qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:) 
     1586            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:) 
     1587            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:) 
    15401588         ENDDO 
    15411589         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:) 
    15421590         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:) 
    1543 !!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0) 
     1591         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:) 
    15441592      ELSE 
    15451593         qns_tot  (:,:  ) = zqns_tot  (:,:  ) 
    15461594         qns_oce  (:,:  ) = zqns_oce  (:,:  ) 
    15471595         qns_ice  (:,:,:) = zqns_ice  (:,:,:) 
    1548          qprec_ice(:,:)   = zqprec_ice(:,:) 
    1549          qemp_oce (:,:)   = zqemp_oce (:,:) 
    1550       ENDIF 
    1551  
    1552       CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce )  
     1596         qevap_ice(:,:,:) = zqevap_ice(:,:,:) 
     1597         qprec_ice(:,:  ) = zqprec_ice(:,:  ) 
     1598         qemp_oce (:,:  ) = zqemp_oce (:,:  ) 
     1599         qemp_ice (:,:  ) = zqemp_ice (:,:  ) 
     1600      ENDIF 
    15531601#else 
    1554       ! 
    15551602      ! clem: this formulation is certainly wrong... but better than it was before... 
    15561603      zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with: 
    … …  
    16191666 
    16201667#if defined key_lim3 
    1621       CALL wrk_alloc( jpi,jpj, zqsr_oce )  
    16221668      ! --- solar flux over ocean --- ! 
    16231669      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax 
    … …  
    16271673      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:) 
    16281674      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF 
    1629  
    1630       CALL wrk_dealloc( jpi,jpj, zqsr_oce )  
    16311675#endif 
    16321676 
    … …  
    16791723      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    16801724 
    1681       CALL wrk_dealloc( jpi,jpj,       zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot ) 
    1682       CALL wrk_dealloc( jpi,jpj,jpl,   zqns_ice, zqsr_ice, zdqns_ice ) 
     1725      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw ) 
     1726      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice ) 
     1727      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
     1728      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
    16831729      ! 
    16841730      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx') 
    … …  
    17191765          
    17201766         IF ( nn_components == jp_iam_opa ) THEN 
    1721             ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of ln_useCT on the received part 
     1767            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part 
    17221768         ELSE 
    17231769            ! we must send the surface potential temperature  
    1724             IF( ln_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) ) 
     1770            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) ) 
    17251771            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem) 
    17261772            ENDIF 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6403 r7278  
    106106      INTEGER  ::    jl                 ! dummy loop index 
    107107      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    108       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled) 
    109108      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice  
    110109      !!---------------------------------------------------------------------- 
    … …  
    193192         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%] 
    194193         !---------------------------------------------------------------------------------------- 
    195          CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     194         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    196195         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    197196 
    … …  
    199198         CASE( jp_clio )                                       ! CLIO bulk formulation 
    200199            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
    201             ! (zalb_ice) is computed within the bulk routine 
    202             CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice ) 
    203             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
    204             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     200            ! (alb_ice) is computed within the bulk routine 
     201                                 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
     202            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     203            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    205204         CASE( jp_core )                                       ! CORE bulk formulation 
    206205            ! albedo depends on cloud fraction because of non-linear spectral effects 
    207             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    208             CALL blk_ice_core_flx( t_su, zalb_ice ) 
    209             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
    210             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     206            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     207                                 CALL blk_ice_core_flx( t_su, alb_ice ) 
     208            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     209            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    211210         CASE ( jp_purecpl ) 
    212211            ! albedo depends on cloud fraction because of non-linear spectral effects 
    213             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    214                                  CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
    215             ! clem: evap_ice is forced to 0 in coupled mode for now  
    216             !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models 
    217             evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp 
    218             IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     212            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     213                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     214            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    219215         END SELECT 
    220          CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     216         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    221217 
    222218         !----------------------------! 
    … …  
    577573      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    578574      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    579       sfx_res(:,:) = 0._wp 
     575      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
    580576      ! 
    581577      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    … …  
    593589      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
    594590      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
    595       hfx_err_dif(:,:) = 0._wp   ; 
     591      hfx_err_dif(:,:) = 0._wp 
     592      wfx_err_sub(:,:) = 0._wp 
    596593      ! 
    597594      afx_tot(:,:) = 0._wp   ; 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7277 r7278  
    323323         emp_b (:,:) = emp (:,:) 
    324324         sfx_b (:,:) = sfx (:,:) 
     325         IF ( ln_rnf ) THEN 
     326            rnf_b    (:,:  ) = rnf    (:,:  ) 
     327            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 
     328         ENDIF 
    325329      ENDIF 
    326330      !                                            ! ---------------------------------------- ! 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r7277 r7278  
    109109      ! 
    110110      CALL wrk_alloc( jpi,jpj, ztfrz) 
    111  
    112       !                                            ! ---------------------------------------- ! 
    113       IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          ! 
    114          !                                         ! ---------------------------------------- ! 
    115          rnf_b    (:,:  ) = rnf    (:,:  )               ! Swap the ocean forcing fields except at nit000 
    116          rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)               ! where before fields are set at the end of the routine 
    117          ! 
    118       ENDIF 
    119  
     111      ! 
    120112      !                                            !-------------------! 
    121113      !                                            !   Update runoff   ! 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r7277 r7278  
    7070         ssu_m(:,:) = ub(:,:,1) 
    7171         ssv_m(:,:) = vb(:,:,1) 
    72          IF( ln_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
     72         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
    7373         ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem) 
    7474         ENDIF 
    … …  
    9292            ssu_m(:,:) = zcoef * ub(:,:,1) 
    9393            ssv_m(:,:) = zcoef * vb(:,:,1) 
    94             IF( ln_useCT )  THEN    ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
     94            IF( l_useCT )  THEN    ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
    9595            ELSE                    ;   sst_m(:,:) = zcoef * zts(:,:,jp_tem) 
    9696            ENDIF 
    … …  
    120120         ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1) 
    121121         ssv_m(:,:) = ssv_m(:,:) + vb(:,:,1) 
    122          IF( ln_useCT )  THEN    ;   sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
     122         IF( l_useCT )  THEN    ;   sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) ) 
    123123         ELSE                    ;   sst_m(:,:) = sst_m(:,:) + zts(:,:,jp_tem) 
    124124         ENDIF 
    … …  
    241241         ssu_m(:,:) = ub(:,:,1) 
    242242         ssv_m(:,:) = vb(:,:,1) 
    243          IF( ln_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) ) 
     243         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) ) 
    244244         ELSE                    ;   sst_m(:,:) = tsn(:,:,1,jp_tem) 
    245245         ENDIF 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r7277 r7278  
    7575 
    7676   !                               !!** Namelist nameos ** 
    77    INTEGER , PUBLIC ::   nn_eos     ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
    78    LOGICAL , PUBLIC ::   ln_useCT   ! determine if eos_pt_from_ct is used to compute sst_m 
     77   LOGICAL , PUBLIC ::   ln_TEOS10   ! determine if eos_pt_from_ct is used to compute sst_m 
     78   LOGICAL , PUBLIC ::   ln_EOS80   ! determine if eos_pt_from_ct is used to compute sst_m 
     79   LOGICAL , PUBLIC ::   ln_SEOS   ! determine if eos_pt_from_ct is used to compute sst_m 
     80 
     81   ! Parameters 
     82   LOGICAL , PUBLIC    ::   l_useCT         ! =T in ln_TEOS10=T (i.e. use eos_pt_from_ct to compute sst_m), =F otherwise 
     83   INTEGER , PUBLIC    ::   neos            ! Identifier for equation of state used 
     84 
     85   INTEGER , PARAMETER ::   np_teos10 = -1  ! parameter for using TEOS10 
     86   INTEGER , PARAMETER ::   np_eos80  =  0  ! parameter for using EOS80 
     87   INTEGER , PARAMETER ::   np_seos   = 1   ! parameter for using Simplified Equation of state 
    7988 
    8089   !                               !!!  simplified eos coefficients (default value: Vallis 2006) 
    … …  
    184193      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    185194      !!       potential temperature and salinity using an equation of state 
    186       !!       defined through the namelist parameter nn_eos. 
     195      !!       selected in the nameos namelist 
    187196      !! 
    188197      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
    … …  
    194203      !!                rau0   reference density            kg/m^3 
    195204      !! 
    196       !!     nn_eos = -1 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
     205      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
    197206      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg 
    198207      !! 
    199       !!     nn_eos =  0 : polynomial EOS-80 equation of state is used for rho(t,s,z). 
     208      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z). 
    200209      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu 
    201210      !! 
    202       !!     nn_eos =  1 : simplified equation of state 
     211      !!     ln_seos : simplified equation of state 
    203212      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0 
    204213      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
    … …  
    224233      IF( nn_timing == 1 )   CALL timing_start('eos-insitu') 
    225234      ! 
    226       SELECT CASE( nn_eos ) 
    227       ! 
    228       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     235      SELECT CASE( neos ) 
     236      ! 
     237      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    229238         ! 
    230239         DO jk = 1, jpkm1 
    … …  
    266275         END DO 
    267276         ! 
    268       CASE( 1 )                !==  simplified EOS  ==! 
     277      CASE( np_seos )                !==  simplified EOS  ==! 
    269278         ! 
    270279         DO jk = 1, jpkm1 
    … …  
    300309      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
    301310      !!      potential volumic mass (Kg/m3) from potential temperature and 
    302       !!      salinity fields using an equation of state defined through the 
    303       !!     namelist parameter nn_eos. 
     311      !!      salinity fields using an equation of state selected in the 
     312      !!     namelist. 
    304313      !! 
    305314      !! ** Action  : - prd  , the in situ density (no units) 
    … …  
    322331      IF( nn_timing == 1 )   CALL timing_start('eos-pot') 
    323332      ! 
    324       SELECT CASE ( nn_eos ) 
    325       ! 
    326       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     333      SELECT CASE ( neos ) 
     334      ! 
     335      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    327336         ! 
    328337         ! Stochastic equation of state 
    … …  
    430439         ENDIF 
    431440          
    432       CASE( 1 )                !==  simplified EOS  ==! 
     441      CASE( np_seos )                !==  simplified EOS  ==! 
    433442         ! 
    434443         DO jk = 1, jpkm1 
    … …  
    467476      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    468477      !!      potential temperature and salinity using an equation of state 
    469       !!      defined through the namelist parameter nn_eos. * 2D field case 
     478      !!      selected in the nameos namelist. * 2D field case 
    470479      !! 
    471480      !! ** Action  : - prd , the in situ density (no units) (unmasked) 
    … …  
    486495      prd(:,:) = 0._wp 
    487496      ! 
    488       SELECT CASE( nn_eos ) 
    489       ! 
    490       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     497      SELECT CASE( neos ) 
     498      ! 
     499      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    491500         ! 
    492501         DO jj = 1, jpjm1 
    … …  
    527536         CALL lbc_lnk( prd, 'T', 1. )                    ! Lateral boundary conditions 
    528537         ! 
    529       CASE( 1 )                !==  simplified EOS  ==! 
     538      CASE( np_seos )                !==  simplified EOS  ==! 
    530539         ! 
    531540         DO jj = 1, jpjm1 
    … …  
    576585      IF( nn_timing == 1 )   CALL timing_start('rab_3d') 
    577586      ! 
    578       SELECT CASE ( nn_eos ) 
    579       ! 
    580       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     587      SELECT CASE ( neos ) 
     588      ! 
     589      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    581590         ! 
    582591         DO jk = 1, jpkm1 
    … …  
    635644         END DO 
    636645         ! 
    637       CASE( 1 )                  !==  simplified EOS  ==! 
     646      CASE( np_seos )                  !==  simplified EOS  ==! 
    638647         ! 
    639648         DO jk = 1, jpkm1 
    … …  
    657666      CASE DEFAULT 
    658667         IF(lwp) WRITE(numout,cform_err) 
    659          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     668         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos 
    660669         nstop = nstop + 1 
    661670         ! 
    … …  
    668677      ! 
    669678   END SUBROUTINE rab_3d 
     679 
    670680 
    671681   SUBROUTINE rab_2d( pts, pdep, pab ) 
    … …  
    690700      pab(:,:,:) = 0._wp 
    691701      ! 
    692       SELECT CASE ( nn_eos ) 
    693       ! 
    694       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     702      SELECT CASE ( neos ) 
     703      ! 
     704      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    695705         ! 
    696706         DO jj = 1, jpjm1 
    … …  
    750760         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                     
    751761         ! 
    752       CASE( 1 )                  !==  simplified EOS  ==! 
     762      CASE( np_seos )                  !==  simplified EOS  ==! 
    753763         ! 
    754764         DO jj = 1, jpjm1 
    … …  
    773783      CASE DEFAULT 
    774784         IF(lwp) WRITE(numout,cform_err) 
    775          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     785         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos 
    776786         nstop = nstop + 1 
    777787         ! 
    … …  
    806816      pab(:) = 0._wp 
    807817      ! 
    808       SELECT CASE ( nn_eos ) 
    809       ! 
    810       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     818      SELECT CASE ( neos ) 
     819      ! 
     820      CASE( np_teos10, np_eos80 )      !==  polynomial TEOS-10 / EOS-80 ==! 
    811821         ! 
    812822         ! 
    … …  
    859869         ! 
    860870         ! 
    861       CASE( 1 )                  !==  simplified EOS  ==! 
     871      CASE( np_seos )                  !==  simplified EOS  ==! 
    862872         ! 
    863873         zt    = pts(jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    864874         zs    = pts(jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    865          zh    = pdep                    ! depth at the partial step level 
     875         zh    = pdep                   ! depth at the partial step level 
    866876         ! 
    867877         zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    … …  
    873883      CASE DEFAULT 
    874884         IF(lwp) WRITE(numout,cform_err) 
    875          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     885         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos 
    876886         nstop = nstop + 1 
    877887         ! 
    … …  
    10051015      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celsius] 
    10061016      ! 
    1007       INTEGER  ::   ji, jj   ! dummy loop indices 
    1008       REAL(wp) ::   zt, zs   ! local scalars 
    1009       !!---------------------------------------------------------------------- 
    1010       ! 
    1011       SELECT CASE ( nn_eos ) 
    1012       ! 
    1013       CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==! 
    1014          ! 
     1017      INTEGER  ::   ji, jj          ! dummy loop indices 
     1018      REAL(wp) ::   zt, zs, z1_S0   ! local scalars 
     1019      !!---------------------------------------------------------------------- 
     1020      ! 
     1021      SELECT CASE ( neos ) 
     1022      ! 
     1023      CASE ( np_teos10, np_seos )      !==  CT,SA (TEOS-10 and S-EOS formulations) ==! 
     1024         ! 
     1025         z1_S0 = 1._wp / 35.16504_wp 
    10151026         DO jj = 1, jpj 
    10161027            DO ji = 1, jpi 
    1017                zs= SQRT( ABS( psal(ji,jj) ) * r1_S0 )           ! square root salinity 
     1028               zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
    10181029               ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
    10191030                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
    … …  
    10241035         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 
    10251036         ! 
    1026       CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==! 
     1037      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==! 
    10271038         ! 
    10281039         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
    … …  
    10331044      CASE DEFAULT 
    10341045         IF(lwp) WRITE(numout,cform_err) 
    1035          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     1046         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos 
    10361047         nstop = nstop + 1 
    10371048         ! 
    … …  
    10391050      ! 
    10401051  END SUBROUTINE eos_fzp_2d 
     1052 
    10411053 
    10421054  SUBROUTINE eos_fzp_0d( psal, ptf, pdep ) 
    … …  
    10591071      !!---------------------------------------------------------------------- 
    10601072      ! 
    1061       SELECT CASE ( nn_eos ) 
    1062       ! 
    1063       CASE ( -1, 1 )                !==  CT,SA (TEOS-10 formulation) ==! 
    1064          ! 
    1065          zs  = SQRT( ABS( psal ) * r1_S0 )           ! square root salinity 
     1073      SELECT CASE ( neos ) 
     1074      ! 
     1075      CASE ( np_teos10, np_seos )      !==  CT,SA (TEOS-10 and S-EOS formulations) ==! 
     1076         ! 
     1077         zs  = SQRT( ABS( psal ) / 35.16504_wp )           ! square root salinity 
    10661078         ptf = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
    10671079                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
    … …  
    10701082         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep 
    10711083         ! 
    1072       CASE ( 0 )                     !==  PT,SP (UNESCO formulation)  ==! 
     1084      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==! 
    10731085         ! 
    10741086         ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   & 
    … …  
    10791091      CASE DEFAULT 
    10801092         IF(lwp) WRITE(numout,cform_err) 
    1081          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     1093         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos 
    10821094         nstop = nstop + 1 
    10831095         ! 
    … …  
    11181130      IF( nn_timing == 1 )   CALL timing_start('eos_pen') 
    11191131      ! 
    1120       SELECT CASE ( nn_eos ) 
    1121       ! 
    1122       CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     1132      SELECT CASE ( neos ) 
     1133      ! 
     1134      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    11231135         ! 
    11241136         DO jk = 1, jpkm1 
    … …  
    11831195         END DO 
    11841196         ! 
    1185       CASE( 1 )                !==  Vallis (2006) simplified EOS  ==! 
     1197      CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==! 
    11861198         ! 
    11871199         DO jk = 1, jpkm1 
    … …  
    12051217      CASE DEFAULT 
    12061218         IF(lwp) WRITE(numout,cform_err) 
    1207          IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     1219         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos 
    12081220         nstop = nstop + 1 
    12091221         ! 
    … …  
    12241236      !!---------------------------------------------------------------------- 
    12251237      INTEGER  ::   ios   ! local integer 
    1226       !! 
    1227       NAMELIST/nameos/ nn_eos, ln_useCT, rn_a0, rn_b0, rn_lambda1, rn_mu1,   & 
     1238      INTEGER  ::   ioptio   ! local integer 
     1239      !! 
     1240      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_a0, rn_b0, rn_lambda1, rn_mu1,   & 
    12281241         &                                             rn_lambda2, rn_mu2, rn_nu 
    12291242      !!---------------------------------------------------------------------- 
    … …  
    12451258         WRITE(numout,*) 'eos_init : equation of state' 
    12461259         WRITE(numout,*) '~~~~~~~~' 
    1247          WRITE(numout,*) '          Namelist nameos : set eos parameters' 
    1248          WRITE(numout,*) '             flag for eq. of state and N^2  nn_eos   = ', nn_eos 
    1249          IF( ln_useCT )   THEN 
    1250             WRITE(numout,*) '             model uses Conservative Temperature' 
    1251             WRITE(numout,*) '             Important: model must be initialized with CT and SA fields' 
    1252          ELSE 
    1253             WRITE(numout,*) '             model does not use Conservative Temperature' 
    1254          ENDIF 
     1260         WRITE(numout,*) '   Namelist nameos : Chosen the Equation Of Seawater (EOS)' 
     1261         WRITE(numout,*) '      TEOS-10 : rho=F(Conservative Temperature, Absolute  Salinity, depth)   ln_TEOS10 = ', ln_TEOS10 
     1262         WRITE(numout,*) '      EOS-80  : rho=F(Potential    Temperature, Practical Salinity, depth)   ln_EOS80  = ', ln_EOS80 
     1263         WRITE(numout,*) '      S-EOS   : rho=F(Conservative Temperature, Absolute  Salinity, depth)   ln_SEOS   = ', ln_SEOS 
    12551264      ENDIF 
    1256       ! 
    1257       SELECT CASE( nn_eos )         ! check option 
    1258       ! 
    1259       CASE( -1 )                       !==  polynomial TEOS-10  ==! 
     1265 
     1266      ! Check options for equation of state & set neos based on logical flags 
     1267      ioptio = 0 
     1268      IF( ln_TEOS10 ) THEN   ;   ioptio = ioptio+1   ;   neos = np_teos10   ;   ENDIF 
     1269      IF( ln_EOS80  ) THEN   ;   ioptio = ioptio+1   ;   neos = np_eos80    ;   ENDIF 
     1270      IF( ln_SEOS   ) THEN   ;   ioptio = ioptio+1   ;   neos = np_seos     ;   ENDIF 
     1271      IF( ioptio /= 1 )   CALL ctl_stop("Exactly one equation of state option must be selected") 
     1272      ! 
     1273      SELECT CASE( neos )         ! check option 
     1274      ! 
     1275      CASE( np_teos10 )                       !==  polynomial TEOS-10  ==! 
    12601276         IF(lwp) WRITE(numout,*) 
    12611277         IF(lwp) WRITE(numout,*) '          use of TEOS-10 equation of state (cons. temp. and abs. salinity)' 
     1278         ! 
     1279         l_useCT = .TRUE.                          ! model temperature is Conservative temperature  
    12621280         ! 
    12631281         rdeltaS = 32._wp 
    … …  
    14461464         BPE002 = 1.7269476440e-04_wp 
    14471465         ! 
    1448       CASE( 0 )                        !==  polynomial EOS-80 formulation  ==! 
     1466      CASE( np_eos80 )                        !==  polynomial EOS-80 formulation  ==! 
    14491467         ! 
    14501468         IF(lwp) WRITE(numout,*) 
    14511469         IF(lwp) WRITE(numout,*) '          use of EOS-80 equation of state (pot. temp. and pract. salinity)' 
    14521470         ! 
     1471         l_useCT = .FALSE.                         ! model temperature is Potential temperature 
    14531472         rdeltaS = 20._wp 
    14541473         r1_S0  = 1._wp/40._wp 
    … …  
    16361655         BPE002 = 5.3661089288e-04_wp 
    16371656         ! 
    1638       CASE( 1 )                        !==  Simplified EOS     ==! 
     1657      CASE( np_seos )                        !==  Simplified EOS     ==! 
    16391658         IF(lwp) THEN 
    16401659            WRITE(numout,*) 
    … …  
    16511670            WRITE(numout,*) '               Caution: rn_beta0=0 incompatible with ddm parameterization ' 
    16521671         ENDIF 
    1653          ! 
    1654       CASE DEFAULT                     !==  ERROR in nn_eos  ==! 
    1655          WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos 
     1672         l_useCT = .TRUE.          ! Use conservative temperature 
     1673         ! 
     1674      CASE DEFAULT                     !==  ERROR in neos  ==! 
     1675         WRITE(ctmp1,*) '          bad flag value for neos = ', neos, '. You should never see this error' 
    16561676         CALL ctl_stop( ctmp1 ) 
    16571677         ! 
    … …  
    16621682      r1_rcp      = 1._wp / rcp 
    16631683      r1_rau0_rcp = 1._wp / rau0_rcp  
     1684      ! 
     1685      IF(lwp) THEN 
     1686         IF( l_useCT )   THEN 
     1687            WRITE(numout,*) '             model uses Conservative Temperature' 
     1688            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields' 
     1689         ELSE 
     1690            WRITE(numout,*) '             model does not use Conservative Temperature' 
     1691         ENDIF 
     1692      ENDIF 
    16641693      ! 
    16651694      IF(lwp) WRITE(numout,*) 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r6140 r7278  
    111111         ELSE                                   ! No restart or restart not found: Euler forward time stepping 
    112112            zfact = 1._wp 
     113            sbc_tsc(:,:,:) = 0._wp 
    113114            sbc_tsc_b(:,:,:) = 0._wp 
    114115         ENDIF 
    … …  
    207208         END DO   
    208209      ENDIF 
     210 
     211      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst 
     212      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss 
     213 
    209214      ! 
    210215      !---------------------------------------- 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r6140 r7278  
    174174                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    175175               ! add to the eddy viscosity coef. previously computed 
     176# if defined key_zdftmx_new 
     177               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx 
     178               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds 
     179# else 
    176180               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 
     181# endif 
    177182               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt 
    178183               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r6140 r7278  
    3131   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3232 
    33    USE eosbn2, ONLY : nn_eos 
     33   USE eosbn2, ONLY : neos 
    3434 
    3535   IMPLICIT NONE 
    … …  
    175175      !  Compute Ekman depth from wind stress forcing. 
    176176      ! ------------------------------------------------------- 
    177       zflageos = ( 0.5 + SIGN( 0.5, nn_eos - 1. ) ) * rau0 
     177      zflageos = ( 0.5 + SIGN( 0.5, neos - 1. ) ) * rau0 
    178178      DO jj = 2, jpjm1 
    179179            DO ji = fs_2, fs_jpim1 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r6140 r7278  
    323323                  zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) ) 
    324324                  !                                           ! TKE Langmuir circulation source term 
    325                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     325                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   & 
     326                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    326327               END DO 
    327328            END DO 
    … …  
    375376            DO ji = fs_2, fs_jpim1   ! vector opt. 
    376377               zcof   = zfact1 * tmask(ji,jj,jk) 
     378# if defined key_zdftmx_new 
     379               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability) 
     380               zzd_up = zcof * MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp )   &  ! upper diagonal 
     381                  &          / (  e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  ) 
     382               zzd_lw = zcof * MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp )   &  ! lower diagonal 
     383                  &          / (  e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  ) 
     384# else 
    377385               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
    378386                  &          / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  ) ) 
    379387               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal 
    380388                  &          / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) ) 
     389# endif 
    381390               !                                   ! shear prod. at w-point weightened by mask 
    382391               zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    … …  
    732741      ! 
    733742      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number 
     743# if defined key_zdftmx_new 
     744      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used 
     745      rn_emin  = 1.e-10_wp 
     746      rmxl_min = 1.e-03_wp 
     747      IF(lwp) THEN                  ! Control print 
     748         WRITE(numout,*) 
     749         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 
     750         WRITE(numout,*) '~~~~~~~~~~~~' 
     751      ENDIF 
     752# else 
    734753      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity 
     754# endif 
    735755      ! 
    736756      IF(lwp) THEN                    !* Control print 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r6140 r7278  
    541541   END SUBROUTINE zdf_tmx_init 
    542542 
     543#elif defined key_zdftmx_new 
     544   !!---------------------------------------------------------------------- 
     545   !!   'key_zdftmx_new'               Internal wave-driven vertical mixing 
     546   !!---------------------------------------------------------------------- 
     547   !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz 
     548   !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz 
     549   !!---------------------------------------------------------------------- 
     550   USE oce            ! ocean dynamics and tracers variables 
     551   USE dom_oce        ! ocean space and time domain variables 
     552   USE zdf_oce        ! ocean vertical physics variables 
     553   USE zdfddm         ! ocean vertical physics: double diffusive mixing 
     554   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     555   USE eosbn2         ! ocean equation of state 
     556   USE phycst         ! physical constants 
     557   USE prtctl         ! Print control 
     558   USE in_out_manager ! I/O manager 
     559   USE iom            ! I/O Manager 
     560   USE lib_mpp        ! MPP library 
     561   USE wrk_nemo       ! work arrays 
     562   USE timing         ! Timing 
     563   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     564 
     565   IMPLICIT NONE 
     566   PRIVATE 
     567 
     568   PUBLIC   zdf_tmx         ! called in step module  
     569   PUBLIC   zdf_tmx_init    ! called in nemogcm module  
     570   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module 
     571 
     572   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: wave-driven mixing flag 
     573 
     574   !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing * 
     575   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 
     576   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency 
     577   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F) 
     578 
     579   REAL(wp) ::  r1_6 = 1._wp / 6._wp 
     580 
     581   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_tmx     ! power available from high-mode wave breaking (W/m2) 
     582   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_tmx     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 
     583   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_tmx     ! power available from low-mode, critical slope wave breaking (W/m2) 
     584   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_tmx     ! WKB decay scale for high-mode energy dissipation (m) 
     585   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_tmx     ! decay scale for low-mode critical slope dissipation (m) 
     586   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_tmx     ! local energy density available for mixing (W/kg) 
     587   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_tmx     ! buoyancy flux Kz * N^2 (W/kg) 
     588   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_tmx    ! vertically integrated buoyancy flux (W/m2) 
     589   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T) 
     590   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity 
     591 
     592   !! * Substitutions 
     593#  include "zdfddm_substitute.h90" 
     594#  include "vectopt_loop_substitute.h90" 
     595   !!---------------------------------------------------------------------- 
     596   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     597   !! $Id$ 
     598   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     599   !!---------------------------------------------------------------------- 
     600CONTAINS 
     601 
     602   INTEGER FUNCTION zdf_tmx_alloc() 
     603      !!---------------------------------------------------------------------- 
     604      !!                ***  FUNCTION zdf_tmx_alloc  *** 
     605      !!---------------------------------------------------------------------- 
     606      ALLOCATE(     ebot_tmx(jpi,jpj),  epyc_tmx(jpi,jpj),  ecri_tmx(jpi,jpj)    ,   & 
     607      &             hbot_tmx(jpi,jpj),  hcri_tmx(jpi,jpj),  emix_tmx(jpi,jpj,jpk),   & 
     608      &         bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk),   &  
     609      &         zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc     ) 
     610      ! 
     611      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc ) 
     612      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays') 
     613   END FUNCTION zdf_tmx_alloc 
     614 
     615 
     616   SUBROUTINE zdf_tmx( kt ) 
     617      !!---------------------------------------------------------------------- 
     618      !!                  ***  ROUTINE zdf_tmx  *** 
     619      !!                    
     620      !! ** Purpose :   add to the vertical mixing coefficients the effect of 
     621      !!              breaking internal waves. 
     622      !! 
     623      !! ** Method  : - internal wave-driven vertical mixing is given by: 
     624      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_tmx /( Nu * N^2 )  ) 
     625      !!              where emix_tmx is the 3D space distribution of the wave-breaking  
     626      !!              energy and Nu the molecular kinematic viscosity. 
     627      !!              The function f(Reb) is linear (constant mixing efficiency) 
     628      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 
     629      !! 
     630      !!              - Compute emix_tmx, the 3D power density that allows to compute 
     631      !!              Reb and therefrom the wave-induced vertical diffusivity. 
     632      !!              This is divided into three components: 
     633      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
     634      !!                     emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx ) 
     635      !!                                   / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx 
     636      !!              where hcri_tmx is the characteristic length scale of the bottom  
     637      !!              intensification, ecri_tmx a map of available power, and H the ocean depth. 
     638      !!                 2. Pycnocline-intensified low-mode dissipation 
     639      !!                     emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     640      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
     641      !!              where epyc_tmx is a map of available power, and nn_zpyc 
     642      !!              is the chosen stratification-dependence of the internal wave 
     643      !!              energy dissipation. 
     644      !!                 3. WKB-height dependent high mode dissipation 
     645      !!                     emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx) 
     646      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) ) 
     647      !!              where hbot_tmx is the characteristic length scale of the WKB bottom  
     648      !!              intensification, ebot_tmx is a map of available power, and z_wkb is the 
     649      !!              WKB-stretched height above bottom defined as 
     650      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) 
     651      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    ) 
     652      !! 
     653      !!              - update the model vertical eddy viscosity and diffusivity:  
     654      !!                     avt  = avt  +    av_wave 
     655      !!                     avm  = avm  +    av_wave 
     656      !!                     avmu = avmu + mi(av_wave) 
     657      !!                     avmv = avmv + mj(av_wave) 
     658      !! 
     659      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing: 
     660      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb) 
     661      !! 
     662      !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing 
     663      !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing     
     664      !! 
     665      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
     666      !!---------------------------------------------------------------------- 
     667      INTEGER, INTENT(in) ::   kt   ! ocean time-step  
     668      ! 
     669      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     670      REAL(wp) ::   ztpc         ! scalar workspace 
     671      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure 
     672      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth 
     673      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom 
     674      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution 
     675      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid) 
     676      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid) 
     677      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter 
     678      !!---------------------------------------------------------------------- 
     679      ! 
     680      IF( nn_timing == 1 )   CALL timing_start('zdf_tmx') 
     681      ! 
     682      CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
     683      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
     684 
     685      !                          ! ----------------------------- ! 
     686      !                          !  Internal wave-driven mixing  !  (compute zav_wave) 
     687      !                          ! ----------------------------- ! 
     688      !                              
     689      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
     690      !                                                 using an exponential decay from the seafloor. 
     691      DO jj = 1, jpj                ! part independent of the level 
     692         DO ji = 1, jpi 
     693            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
     694            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) )  ) 
     695            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj) 
     696         END DO 
     697      END DO 
     698 
     699      DO jk = 2, jpkm1              ! complete with the level-dependent part 
     700         emix_tmx(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_tmx(:,:) )                      & 
     701            &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) )  ) * wmask(:,:,jk)   & 
     702            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     703      END DO 
     704 
     705      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
     706      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
     707 
     708      SELECT CASE ( nn_zpyc ) 
     709 
     710      CASE ( 1 )               ! Dissipation scales as N (recommended) 
     711 
     712         zfact(:,:) = 0._wp 
     713         DO jk = 2, jpkm1              ! part independent of the level 
     714            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     715         END DO 
     716 
     717         DO jj = 1, jpj 
     718            DO ji = 1, jpi 
     719               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     720            END DO 
     721         END DO 
     722 
     723         DO jk = 2, jpkm1              ! complete with the level-dependent part 
     724            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     725         END DO 
     726 
     727      CASE ( 2 )               ! Dissipation scales as N^2 
     728 
     729         zfact(:,:) = 0._wp 
     730         DO jk = 2, jpkm1              ! part independent of the level 
     731            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     732         END DO 
     733 
     734         DO jj= 1, jpj 
     735            DO ji = 1, jpi 
     736               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     737            END DO 
     738         END DO 
     739 
     740         DO jk = 2, jpkm1              ! complete with the level-dependent part 
     741            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     742         END DO 
     743 
     744      END SELECT 
     745 
     746      !                        !* WKB-height dependent mixing: distribute energy over the time-varying  
     747      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
     748       
     749      zwkb(:,:,:) = 0._wp 
     750      zfact(:,:) = 0._wp 
     751      DO jk = 2, jpkm1 
     752         zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     753         zwkb(:,:,jk) = zfact(:,:) 
     754      END DO 
     755 
     756      DO jk = 2, jpkm1 
     757         DO jj = 1, jpj 
     758            DO ji = 1, jpi 
     759               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
     760                                            &           * tmask(ji,jj,jk) / zfact(ji,jj) 
     761            END DO 
     762         END DO 
     763      END DO 
     764      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
     765 
     766      zweight(:,:,:) = 0._wp 
     767      DO jk = 2, jpkm1 
     768         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk)                    & 
     769            &   * (  EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) )  ) 
     770      END DO 
     771 
     772      zfact(:,:) = 0._wp 
     773      DO jk = 2, jpkm1              ! part independent of the level 
     774         zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
     775      END DO 
     776 
     777      DO jj = 1, jpj 
     778         DO ji = 1, jpi 
     779            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     780         END DO 
     781      END DO 
     782 
     783      DO jk = 2, jpkm1              ! complete with the level-dependent part 
     784         emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
     785            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     786      END DO 
     787 
     788 
     789      ! Calculate molecular kinematic viscosity 
     790      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
     791         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
     792      DO jk = 2, jpkm1 
     793         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
     794      END DO 
     795 
     796      ! Calculate turbulence intensity parameter Reb 
     797      DO jk = 2, jpkm1 
     798         zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
     799      END DO 
     800 
     801      ! Define internal wave-induced diffusivity 
     802      DO jk = 2, jpkm1 
     803         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     804      END DO 
     805 
     806      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
     807         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
     808            DO jj = 1, jpj 
     809               DO ji = 1, jpi 
     810                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
     811                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     812                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 
     813                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     814                  ENDIF 
     815               END DO 
     816            END DO 
     817         END DO 
     818      ENDIF 
     819 
     820      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
     821         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
     822      END DO 
     823 
     824      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
     825         ztpc = 0._wp 
     826         DO jk = 2, jpkm1 
     827            DO jj = 1, jpj 
     828               DO ji = 1, jpi 
     829                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
     830                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     831               END DO 
     832            END DO 
     833         END DO 
     834         IF( lk_mpp )   CALL mpp_sum( ztpc ) 
     835         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     836  
     837         IF(lwp) THEN 
     838            WRITE(numout,*) 
     839            WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)' 
     840            WRITE(numout,*) '~~~~~~~ ' 
     841            WRITE(numout,*) 
     842            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW' 
     843         ENDIF 
     844      ENDIF 
     845 
     846      !                          ! ----------------------- ! 
     847      !                          !   Update  mixing coefs  !                           
     848      !                          ! ----------------------- ! 
     849      !       
     850      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
     851         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb 
     852            DO jj = 1, jpj 
     853               DO ji = 1, jpi 
     854                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  & 
     855                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   & 
     856                      &                 ) * wmask(ji,jj,jk) 
     857               END DO 
     858            END DO 
     859         END DO 
     860         CALL iom_put( "av_ratio", zav_ratio ) 
     861         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
     862            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
     863            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     864            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     865         END DO 
     866         ! 
     867      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
     868         DO jk = 2, jpkm1 
     869            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     870            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     871            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     872         END DO 
     873      ENDIF 
     874 
     875      DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points 
     876         DO jj = 2, jpjm1 
     877            DO ji = fs_2, fs_jpim1  ! vector opt. 
     878               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     879               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
     880            END DO 
     881         END DO 
     882      END DO 
     883      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition 
     884 
     885      !                             !* output internal wave-driven mixing coefficient 
     886      CALL iom_put( "av_wave", zav_wave ) 
     887                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx),  
     888                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx) 
     889      IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN 
     890         bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
     891         pcmap_tmx(:,:) = 0._wp 
     892         DO jk = 2, jpkm1 
     893            pcmap_tmx(:,:) = pcmap_tmx(:,:) + e3w_n(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk) 
     894         END DO 
     895         pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:) 
     896         CALL iom_put( "bflx_tmx", bflx_tmx ) 
     897         CALL iom_put( "pcmap_tmx", pcmap_tmx ) 
     898      ENDIF 
     899      CALL iom_put( "bn2", rn2 ) 
     900      CALL iom_put( "emix_tmx", emix_tmx ) 
     901       
     902      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
     903      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
     904 
     905      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
     906      ! 
     907      IF( nn_timing == 1 )   CALL timing_stop('zdf_tmx') 
     908      ! 
     909   END SUBROUTINE zdf_tmx 
     910 
     911 
     912   SUBROUTINE zdf_tmx_init 
     913      !!---------------------------------------------------------------------- 
     914      !!                  ***  ROUTINE zdf_tmx_init  *** 
     915      !!                      
     916      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading 
     917      !!              of input power maps and decay length scales in netcdf files. 
     918      !! 
     919      !! ** Method  : - Read the namzdf_tmx namelist and check the parameters 
     920      !! 
     921      !!              - Read the input data in NetCDF files : 
     922      !!              power available from high-mode wave breaking (mixing_power_bot.nc) 
     923      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc) 
     924      !!              power available from critical slope wave-breaking (mixing_power_cri.nc) 
     925      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc) 
     926      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc) 
     927      !! 
     928      !! ** input   : - Namlist namzdf_tmx 
     929      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc, 
     930      !!              decay_scale_bot.nc decay_scale_cri.nc 
     931      !! 
     932      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter 
     933      !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx 
     934      !! 
     935      !! References : de Lavergne et al. 2015, JPO; 2016, in prep. 
     936      !!          
     937      !!---------------------------------------------------------------------- 
     938      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     939      INTEGER  ::   inum         ! local integer 
     940      INTEGER  ::   ios 
     941      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars 
     942      !! 
     943      NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff 
     944      !!---------------------------------------------------------------------- 
     945      ! 
     946      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init') 
     947      ! 
     948      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing 
     949      READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901) 
     950901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp ) 
     951      ! 
     952      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing 
     953      READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 ) 
     954902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp ) 
     955      IF(lwm) WRITE ( numond, namzdf_tmx_new ) 
     956      ! 
     957      IF(lwp) THEN                  ! Control print 
     958         WRITE(numout,*) 
     959         WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing' 
     960         WRITE(numout,*) '~~~~~~~~~~~~' 
     961         WRITE(numout,*) '   Namelist namzdf_tmx_new : set wave-driven mixing parameters' 
     962         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc 
     963         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar 
     964         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff 
     965      ENDIF 
     966       
     967      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and 
     968      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should  
     969      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 
     970      avmb(:) = 1.4e-6_wp        ! viscous molecular value 
     971      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_tmx)     
     972      avtb_2d(:,:) = 1.e0_wp     ! uniform  
     973      IF(lwp) THEN                  ! Control print 
     974         WRITE(numout,*) 
     975         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   & 
     976            &               'the viscous molecular value & a very small diffusive value, resp.' 
     977      ENDIF 
     978       
     979      IF( .NOT.lk_zdfddm )   CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' ) 
     980       
     981      !                             ! allocate tmx arrays 
     982      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' ) 
     983      ! 
     984      !                             ! read necessary fields 
     985      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2] 
     986      CALL iom_get  (inum, jpdom_data, 'field', ebot_tmx, 1 )  
     987      CALL iom_close(inum) 
     988      ! 
     989      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2] 
     990      CALL iom_get  (inum, jpdom_data, 'field', epyc_tmx, 1 ) 
     991      CALL iom_close(inum) 
     992      ! 
     993      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2] 
     994      CALL iom_get  (inum, jpdom_data, 'field', ecri_tmx, 1 ) 
     995      CALL iom_close(inum) 
     996      ! 
     997      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m] 
     998      CALL iom_get  (inum, jpdom_data, 'field', hbot_tmx, 1 ) 
     999      CALL iom_close(inum) 
     1000      ! 
     1001      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m] 
     1002      CALL iom_get  (inum, jpdom_data, 'field', hcri_tmx, 1 ) 
     1003      CALL iom_close(inum) 
     1004 
     1005      ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:) 
     1006      epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:) 
     1007      ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:) 
     1008 
     1009      ! Set once for all to zero the first and last vertical levels of appropriate variables 
     1010      emix_tmx (:,:, 1 ) = 0._wp 
     1011      emix_tmx (:,:,jpk) = 0._wp 
     1012      zav_ratio(:,:, 1 ) = 0._wp 
     1013      zav_ratio(:,:,jpk) = 0._wp 
     1014      zav_wave (:,:, 1 ) = 0._wp 
     1015      zav_wave (:,:,jpk) = 0._wp 
     1016 
     1017      zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) ) 
     1018      zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) ) 
     1019      zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) ) 
     1020      IF(lwp) THEN 
     1021         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW' 
     1022         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW' 
     1023         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW' 
     1024      ENDIF 
     1025      ! 
     1026      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init') 
     1027      ! 
     1028   END SUBROUTINE zdf_tmx_init 
     1029 
    5431030#else 
    5441031   !!---------------------------------------------------------------------- 

  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7277 r7278  
    112112      ! Update stochastic parameters and random T/S fluctuations 
    113113      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    114                          CALL sto_par( kstp )          ! Stochastic parameters 
     114      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters 
     115      IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations 
    115116 
    116117      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    154155      ! 
    155156      IF( l_ldfslp ) THEN                             ! slope of lateral mixing 
    156 !!gm : why this here ???? 
    157          IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations 
    158 !!gm 
    159157                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density 
    160158 
    … …  
    183181      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    184182                            CALL wzv           ( kstp )  ! now cross-level velocity  
    185 !!gm : why also here ???? 
    186       IF( ln_sto_eos    )   CALL sto_pts( tsn )                             ! Random T/S fluctuations 
    187 !!gm 
    188183                            CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
    189184                             
    … …  
    305300!!jc: That would be better, but see comment above 
    306301!! 
    307       IF( lrst_oce   )   CALL rst_write     ( kstp )  ! write output ocean restart file 
     302      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file 
     303      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters 
    308304 
    309305#if defined key_agrif 
Note: See TracChangeset for help on using the changeset viewer.