Changeset 15294

Timestamp:

2021-09-27T18:40:59+02:00 (20 months ago)

Author:

mattmartin

Message:

First version implementing stopack code from the standard GO6 FOAM branch.

Location:

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC

Files:

• DYN/dynldf.F90 (modified) (3 diffs)
• SBC/albedo.F90 (modified) (18 diffs)
• SBC/sbcblk_core.F90 (modified) (42 diffs)
• SBC/sbcrnf.F90 (modified) (3 diffs)
• SBC/sbcssr.F90 (modified) (10 diffs)
• STO/stopack.F90 (added)
• TRA/trabbc.F90 (modified) (9 diffs)
• TRA/trabbl.F90 (modified) (10 diffs)
• TRA/traldf.F90 (modified) (3 diffs)
• TRA/traqsr.F90 (modified) (35 diffs)
• TRD/trddyn.F90 (modified) (2 diffs)
• TRD/trdtra.F90 (modified) (2 diffs)
• ZDF/zdfbfr.F90 (modified) (8 diffs)
• ZDF/zdfevd.F90 (modified) (7 diffs)
• ZDF/zdfgls.F90 (modified) (40 diffs)
• ZDF/zdftke.F90 (modified) (43 diffs)
• nemogcm.F90 (modified) (2 diffs)
• step.F90 (modified) (18 diffs)
• step_oce.F90 (modified) (2 diffs)

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

@@ -30 +30 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -38 +39 @@
 
    INTEGER ::   nldf = -2   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals)
+
+#if defined key_dynldf_c3d
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahm10, ahm20, ahm30, ahm40
+#endif
+#if defined key_dynldf_c2d
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: ahm10, ahm20, ahm30, ahm40
+#endif
 
    !! * Substitutions
@@ -67 +75 @@
          ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
+
+#if defined key_dynldf_c3d
+         IF( kt == nit000 .AND. ln_stopack .AND. &
+            & ( (nn_spp_ahm1 + nn_spp_ahm2) > 0 ) ) THEN
+             ALLOCATE ( ahm10(jpi,jpj,jpk), ahm20(jpi,jpj,jpk) )
+             ALLOCATE ( ahm30(jpi,jpj,jpk), ahm40(jpi,jpj,jpk) )
+             ahm10 = ahm1
+             ahm20 = ahm2
+             ahm30 = ahm3
+             ahm40 = ahm4
+         ENDIF
+#endif
+#if defined key_dynldf_c2d
+         IF( kt == nit000 .AND. ln_stopack .AND. &
+            & ( (nn_spp_ahm1 + nn_spp_ahm2) > 0 ) ) THEN
+             ALLOCATE ( ahm10(jpi,jpj), ahm20(jpi,jpj) )
+             ALLOCATE ( ahm30(jpi,jpj), ahm40(jpi,jpj) )
+             ahm10 = ahm1
+             ahm20 = ahm2
+             ahm30 = ahm3
+             ahm40 = ahm4
+         ENDIF
+#endif
+
+#if defined key_traldf_c3d || defined key_traldf_c2d
+         IF( ln_stopack ) THEN
+            IF( nn_spp_ahm1 > 0 ) THEN
+               IF( ln_dynldf_lap ) THEN
+                  ahm1 = ahm10
+                  CALL spp_ahm(kt, ahm1, nn_spp_ahm1, rn_ahm1_sd, jk_spp_ahm1)
+               ENDIF
+               IF( ln_dynldf_bilap ) THEN
+                  ahm3 = ahm30
+                  CALL spp_ahm(kt, ahm3, nn_spp_ahm1, rn_ahm1_sd, jk_spp_ahm3)
+               ENDIF
+            ENDIF
+            IF( nn_spp_ahm2 > 0 ) THEN
+               IF( ln_dynldf_lap ) THEN
+               ahm2 = ahm20
+               CALL spp_ahm(kt, ahm2, nn_spp_ahm2, rn_ahm2_sd, jk_spp_ahm2)
+               ENDIF
+               IF( ln_dynldf_bilap ) THEN
+                  ahm4 = ahm40
+                  CALL spp_ahm(kt, ahm4, nn_spp_ahm2, rn_ahm2_sd, jk_spp_ahm4)
+               ENDIF
+            ENDIF
+         ENDIF
+#endif
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

@@ -21 +21 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE stopack
 
    IMPLICIT NONE
@@ -30 +31 @@
 
    INTEGER  ::   albd_init = 0      !: control flag for initialization
-  
+
    REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
    REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
@@ -36 +37 @@
    REAL(wp) ::   c2       = 0.10    !  "        "
    REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
- 
+
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
@@ -51 +52 @@
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_ice  ***
-      !!          
-      !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!       
+      !!
+      !! ** Purpose :   Computation of the albedo of the snow/ice system
+      !!
       !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
       !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
@@ -72 +73 @@
       !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
       !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
-      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger
       !!                     under melting conditions than under freezing conditions
-      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows
       !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
       !!
       !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
       !!                Brandt et al. 2005, J. Climate, vol 18
-      !!                Grenfell & Perovich 2004, JGR, vol 109 
+      !!                Grenfell & Perovich 2004, JGR, vol 109
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -96 +97 @@
 
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-      
+
       CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
 
-      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
-
-      
+      IF( albd_init == 0 )   CALL albedo_init      ! initialization
+
+
       SELECT CASE ( nn_ice_alb )
 
@@ -108 +109 @@
       !------------------------------------------
       CASE( 0 )
-       
+
          ralb_sf = 0.80       ! dry snow
          ralb_sm = 0.65       ! melting snow
          ralb_if = 0.72       ! bare frozen ice
-         ralb_im = rn_albice  ! bare puddled ice 
-         
+         ralb_im = rn_albice  ! bare puddled ice
+
          !  Computation of ice albedo (free of snow)
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
          ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
          END  WHERE
-      
+
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
@@ -126 +127 @@
          ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
          END WHERE
-     
+
          DO jl = 1, ijpl
             DO jj = 1, jpj
@@ -132 +133 @@
                   ! freezing snow
                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
-                  !                                        !  freezing snow        
+                  !                                        !  freezing snow
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
                   zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
                      &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
-                     &        +         zswitch   * ralb_sf  
+                     &        +         zswitch   * ralb_sf
 
                   ! melting snow
@@ -142 +143 @@
                   zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
                   zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
-                      &     +         zswitch   *   ralb_sm 
+                      &     +         zswitch   *   ralb_sm
                   !
                   ! snow albedo
-                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )
                   zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-               
+
                   ! Ice/snow albedo
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
@@ -154 +155 @@
                END DO
             END DO
+
+#if defined key_traldf_c2d || key_traldf_c3d
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
          END DO
 
@@ -161 +171 @@
       !  New parameterization (2016)
       !------------------------------------------
-      CASE( 1 ) 
+      CASE( 1 )
 
          ralb_im = rn_albice  ! bare puddled ice
@@ -176 +186 @@
 !         ralb_sm = 0.82      ! melting snow
 !         ralb_if = 0.54      ! bare frozen ice
-! 
+!
          !  Computation of ice albedo (free of snow)
-         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )
          z1_c2 = 1. / 0.05
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
          ELSE WHERE                                              ;   zalb = ralb_if
          END  WHERE
-         
+
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
@@ -200 +210 @@
 
                    ! snow albedo
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )
                   zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
 
-                  ! Ice/snow albedo   
+                  ! Ice/snow albedo
                   zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
                   pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
@@ -209 +219 @@
               END DO
             END DO
+
+#if defined key_traldf_c2d || key_traldf_c3d
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
          END DO
          ! Effect of the clouds (2d order polynomial)
-         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );
 
       END SELECT
-      
+
       CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
       !
@@ -223 +242 @@
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_oce  ***
-      !! 
+      !!
       !! ** Purpose :   Computation of the albedo of the ocean
       !!----------------------------------------------------------------------
@@ -229 +248 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) :: zcoef 
+      REAL(wp) :: zcoef
       !!----------------------------------------------------------------------
       !
       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
-      pa_oce_cs(:,:) = zcoef 
+      pa_oce_cs(:,:) = zcoef
       pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
@@ -248 +267 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice
       !!----------------------------------------------------------------------
       !

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -51 +51 @@
    USE par_ice_2
 #endif
+   USE stopack
 
    IMPLICIT NONE
@@ -89 +90 @@
    REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem)
    REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
+   REAL(wp), ALLOCATABLE, SAVE ::   rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
+   REAL(wp), PUBLIC :: rn_sfac ! multiplication factor for snow precipitation over sea-ice
 
    !! * Substitutions
@@ -151 +154 @@
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
          &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
-         &                  sn_tdif, rn_zqt,  rn_zu
+         &                  sn_tdif, rn_zqt,  rn_zu, rn_sfac
       !!---------------------------------------------------------------------
       !
@@ -158 +161 @@
          !                                      ! ====================== !
          !
+         rn_sfac = 1._wp       ! Default to one if missing from namelist
          REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
          READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
@@ -196 +200 @@
          sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
-      ENDIF
+         ALLOCATE ( rn_vfac0(jpi,jpj) )
+         rn_vfac0(:,:) = rn_vfac
+         !
+      ENDIF
+
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN
+         rn_vfac0(:,:) = rn_vfac
+         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw )
+      ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_relw > 0 ) &
+         & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
 
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
@@ -205 +223 @@
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
-         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
+         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)
          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
          ENDIF
-         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
+         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
@@ -219 +237 @@
       !
    END SUBROUTINE sbc_blk_core
-   
-   
+
+
    SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv )
       !!---------------------------------------------------------------------
@@ -230 +248 @@
       !! ** Method  :   CORE bulk formulea for the ocean using atmospheric
       !!      fields read in sbc_read
-      !! 
+      !!
       !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
       !!              - vtau    : j-component of the stress at V-point  (N/m2)
@@ -268 +286 @@
       ! local scalars ( place there for vector optimisation purposes)
       zcoef_qsatw = 0.98 * 640380. / rhoa
-      
+
       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
 
@@ -276 +294 @@
 
       ! ... components ( U10m - U_oce ) at T-point (unmasked)
-      zwnd_i(:,:) = 0.e0  
+      zwnd_i(:,:) = 0.e0
       zwnd_j(:,:) = 0.e0
 #if defined key_cyclone
@@ -289 +307 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
-            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
+            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
+            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
          END DO
       END DO
@@ -320 +338 @@
       CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
          &               Cd, Ch, Ce, zt_zu, zq_zu )
-    
+
       ! ... tau module, i and j component
       DO jj = 1, jpj
@@ -351 +369 @@
       CALL lbc_lnk( vtau(:,:), 'V', -1. )
 
-    
+
       !  Turbulent fluxes over ocean
       ! -----------------------------
@@ -376 +394 @@
          CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ')
       ENDIF
-       
+
       ! ----------------------------------------------------------------------------- !
       !     III    Total FLUXES                                                       !
@@ -384 +402 @@
          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
       !
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
@@ -425 +443 @@
       !
    END SUBROUTINE blk_oce_core
- 
-   
+
+
 #if defined key_lim2 || defined key_lim3
    SUBROUTINE blk_ice_core_tau
@@ -465 +483 @@
             DO ji = 2, jpim1   ! B grid : NO vector opt
                ! ... scalar wind at I-point (fld being at T-point)
-               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj)
-               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
+               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)    &
+                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) &
+                  &      - rn_vfac0(ji,jj) * u_ice(ji,jj)
+               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)    &
+                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) &
+                  &      - rn_vfac0(ji,jj) * v_ice(ji,jj)
                zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
                ! ... ice stress at I-point
@@ -474 +494 @@
                vtau_ice(ji,jj) = zwnorm_f * zwndj_f
                ! ... scalar wind at T-point (fld being at T-point)
-               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
-                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
-               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
-                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
+               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1)                                         &
+                  &      - rn_vfac0(ji,jj) * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
+                  &                                  + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
+               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1)                                         &
+                  &      - rn_vfac0(ji,jj) * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
+                  &                                  + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
                wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
@@ -488 +510 @@
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vect. opt.
-               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
-               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
+               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
+               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
                wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
@@ -495 +517 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
-               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
-               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
+               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )             &
+                  &              * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) &
+                  &                  - rn_vfac0(ji,jj) * u_ice(ji,jj) )
+               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )             &
+                  &              * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) &
+                  &                  - rn_vfac0(ji,jj) * v_ice(ji,jj) )
             END DO
          END DO
@@ -513 +537 @@
 
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau')
-      
+
    END SUBROUTINE blk_ice_core_tau
 
@@ -526 +550 @@
       !!                between atmosphere and sea-ice using CORE bulk
       !!                formulea, ice variables and read atmmospheric fields.
-      !! 
+      !!
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
@@ -546 +570 @@
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx')
       !
-      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
+      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
 
       ! local scalars ( place there for vector optimisation purposes)
@@ -569 +593 @@
                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
-               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
+               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
 
                ! ----------------------------!
@@ -579 +603 @@
                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &
                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
               ! Latent heat sensitivity for ice (Dqla/Dt)
@@ -610 +634 @@
 
 #if defined  key_lim3
-      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
+      CALL wrk_alloc( jpi,jpj, zevap, zsnw )
 
       ! --- evaporation --- !
@@ -620 +644 @@
       ! --- evaporation minus precipitation --- !
       zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
+      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
@@ -643 +667 @@
       DO jl = 1, jpl
          qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
-                                   ! But we do not have Tice => consider it at 0°C => evap=0 
+                                   ! But we do not have Tice => consider it at 0 degC => evap=0
       END DO
 
-      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
+      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )
 #endif
 
@@ -670 +694 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx')
-      
+
    END SUBROUTINE blk_ice_core_flx
 #endif
@@ -683 +707 @@
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
       !!
-      !! ** Method : Monin Obukhov Similarity Theory 
+      !! ** Method : Monin Obukhov Similarity Theory
       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)
       !!
@@ -693 +717 @@
       !!    - better first guess of stability by checking air-sea difference of virtual temperature
       !!       rather than temperature difference only...
-      !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
+      !!    - added function "cd_neutral_10m" that uses the improved parametrization of
       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!
       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them
@@ -728 +752 @@
 
       IF( nn_timing == 1 )  CALL timing_start('turb_core_2z')
-    
+
       CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd )
       CALL wrk_alloc( jpi,jpj, zeta_u, stab )
@@ -740 +764 @@
       U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
 
-      !! First guess of stability: 
+      !! First guess of stability:
       ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt
       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
@@ -754 +778 @@
       Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 )
       Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
-    
+
       !! Initializing transf. coeff. with their first guess neutral equivalents :
       Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10
@@ -765 +789 @@
          !
          ztmp1 = T_zu - sst   ! Updating air/sea differences
-         ztmp2 = q_zu - q_sat 
+         ztmp2 = q_zu - q_sat
 
          ! Updating turbulent scales :   (L&Y 2004 eq. (7))
          ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta*
          ztmp2  = Ce/sqrt_Cd*ztmp2    ! q*
-       
+
          ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu
 
          ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
-         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 
+         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)
          !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu)
 
@@ -781 +805 @@
          zpsi_h_u = psi_h( zeta_u )
          zpsi_m_u = psi_m( zeta_u )
-       
+
          !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
          IF ( .NOT. l_zt_equal_zu ) THEN
@@ -790 +814 @@
             q_zu = max(0., q_zu)
          END IF
-       
+
          IF( ln_cdgw ) THEN      ! surface wave case
-            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
+            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )
             Cd      = sqrt_Cd * sqrt_Cd
          ELSE
@@ -802 +826 @@
            ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10
            sqrt_Cd_n10 = sqrt(ztmp0)
-       
+
            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b)
            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
@@ -809 +833 @@
            !! Update of transfer coefficients:
            ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a)
-           Cd      = ztmp0 / ( ztmp1*ztmp1 )   
+           Cd      = ztmp0 / ( ztmp1*ztmp1 )
            sqrt_Cd = SQRT( Cd )
          ENDIF
@@ -815 +839 @@
          ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
          ztmp2 = sqrt_Cd / sqrt_Cd_n10
-         ztmp1 = 1. + Ch_n10*ztmp0               
+         ztmp1 = 1. + Ch_n10*ztmp0
          Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
          !
-         ztmp1 = 1. + Ce_n10*ztmp0               
+         ztmp1 = 1. + Ce_n10*ztmp0
          Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
          !
@@ -836 +860 @@
    FUNCTION cd_neutral_10m( zw10 )
       !!----------------------------------------------------------------------
-      !! Estimate of the neutral drag coefficient at 10m as a function 
+      !! Estimate of the neutral drag coefficient at 10m as a function
       !! of neutral wind  speed at 10m
       !!
@@ -842 +866 @@
       !!
       !! Author: L. Brodeau, june 2014
-      !!----------------------------------------------------------------------    
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s)
       REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m
       !
       REAL(wp), DIMENSION(:,:), POINTER ::   rgt33
-      !!----------------------------------------------------------------------    
+      !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj, rgt33 )
       !
       !! When wind speed > 33 m/s => Cyclone conditions => special treatment
-      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1  
+      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1
       cd_neutral_10m = 1.e-3 * ( &
          &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33.

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

@@ -54 +54 @@
    REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s]
    REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_avt_rnf0
 
    LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis
@@ -159 +160 @@
             rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
          ENDIF
- 
+
          IF(lwp .AND. lflush) CALL flush(numout)
 
@@ -465 +466 @@
          !                                      !    - mixed upstream-centered (ln_traadv_cen2=T)
          !
+         ALLOCATE ( rn_avt_rnf0(jpi,jpj) )
+         rn_avt_rnf0(:,:) = rn_avt_rnf
+         !
          IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
             &                                              'be spread through depth by ln_rnf_depth'               )

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

@@ -24 +24 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE stopack
+   USE wrk_nemo       ! Memory Allocation
 
    IMPLICIT NONE
@@ -40 +42 @@
    REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS
    REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS
-   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
+   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term
    REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
 
@@ -75 +77 @@
       REAL(wp) ::   zerp     ! local scalar for evaporation damping
       REAL(wp) ::   zqrp     ! local scalar for heat flux damping
-      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
       REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
+      REAL(wp), POINTER, DIMENSION(:,:) :: rn_dqdt_s, zsrp
       INTEGER  ::   ierror   ! return error code
       !!
@@ -95 +97 @@
             !
             IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
+
+               CALL wrk_alloc( jpi, jpj, rn_dqdt_s )
+               rn_dqdt_s(:,:) = rn_dqdt
+
+#if defined key_traldf_c2d || key_traldf_c3d
+            IF( ln_stopack .AND. nn_spp_dqdt > 0 ) &
+               & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt )
+#else
+            IF ( ln_stopack .AND. nn_spp_dqdt > 0 ) &
+               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
+
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
+                     zqrp = rn_dqdt_s(ji,jj) * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
                      qns(ji,jj) = qns(ji,jj) + zqrp
                      qrp(ji,jj) = zqrp
@@ -103 +118 @@
                END DO
                CALL iom_put( "qrp", qrp )                             ! heat flux damping
+               CALL wrk_dealloc( jpi, jpj, rn_dqdt_s )
             ENDIF
             !
             IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
-               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
+               CALL wrk_alloc( jpi, jpj, zsrp)
+               zsrp(:,:) = rn_deds
+#if defined key_traldf_c2d || key_traldf_c3d
+               IF( ln_stopack .AND. nn_spp_dedt > 0 ) &
+                  & CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
 !CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
-                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
+                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )
                      sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
                      erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
@@ -117 +144 @@
                END DO
                CALL iom_put( "erp", erp )                             ! freshwater flux damping
+               CALL wrk_dealloc( jpi,jpj, zsrp )
                !
             ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
-               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
-               zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
+               CALL wrk_alloc( jpi, jpj, zsrp)
+               zsrp(:,:) = rn_deds
+#if defined key_traldf_c2d || key_traldf_c3d
+               IF( ln_stopack .AND. nn_spp_dedt > 0 ) &
+                  & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
+#else
+               IF ( ln_stopack .AND. nn_spp_dedt > 0 ) &
+                  & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                   'key_traldf_c2d or key_traldf_c3d')
+#endif
+               zerp_bnd = rn_sssr_bnd / rday                          !       -              -
 !CDIR COLLAPSE
                DO jj = 1, jpj
-                  DO ji = 1, jpi                            
-                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                  DO ji = 1, jpi
+                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
                         &        / MAX(  sss_m(ji,jj), 1.e-20   )
@@ -134 +171 @@
                END DO
                CALL iom_put( "erp", erp )                             ! freshwater flux damping
+               CALL wrk_dealloc( jpi,jpj,zsrp )
             ENDIF
             !
@@ -144 +182 @@
    END SUBROUTINE sbc_ssr
 
- 
+
    SUBROUTINE sbc_ssr_init
       !!---------------------------------------------------------------------
@@ -167 +205 @@
       !!----------------------------------------------------------------------
       !
- 
-      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist : 
+
+      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist :
       READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp )
@@ -224 +262 @@
       ENDIF
       !
-      !                            !* Initialize qrp and erp if no restoring 
+      !                            !* Initialize qrp and erp if no restoring
       IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
       IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp
       !
    END SUBROUTINE sbc_ssr_init
-      
+
    !!======================================================================
 END MODULE sbcssr

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

@@ -12 +12 @@
 
    !!----------------------------------------------------------------------
-   !!   tra_bbc      : update the tracer trend at ocean bottom 
+   !!   tra_bbc      : update the tracer trend at ocean bottom
    !!   tra_bbc_init : initialization of geothermal heat flux trend
    !!----------------------------------------------------------------------
@@ -19 +19 @@
    USE phycst          ! physical constants
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O manager
@@ -28 +28 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -41 +42 @@
 
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd0   ! geothermal heating trend
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read)
- 
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -56 +58 @@
       !!                  ***  ROUTINE tra_bbc  ***
       !!
-      !! ** Purpose :   Compute the bottom boundary contition on temperature 
-      !!              associated with geothermal heating and add it to the 
+      !! ** Purpose :   Compute the bottom boundary contition on temperature
+      !!              associated with geothermal heating and add it to the
       !!              general trend of temperature equations.
       !!
-      !! ** Method  :   The geothermal heat flux set to its constant value of 
+      !! ** Method  :   The geothermal heat flux set to its constant value of
       !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999).
       !!       The temperature trend associated to this heat flux through the
@@ -89 +91 @@
       !
       !                             !  Add the geothermal heat flux trend on temperature
+
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack .AND. nn_spp_geot > 0) THEN
+          qgh_trd1(:,:) = qgh_trd0(:,:)
+          CALL spp_gen(kt, qgh_trd1, nn_spp_geot, rn_geot_sd, jk_spp_geot)
+      ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_geot > 0 ) &
+         & CALL ctl_stop( 'tra_bbc: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
             ik = mbkt(ji,jj)
-            zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik)
+            zqgh_trd = qgh_trd1(ji,jj) / fse3t(ji,jj,ik)
             tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd
          END DO
@@ -137 +152 @@
       CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files
       !
-      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
+      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir
       !!----------------------------------------------------------------------
 
@@ -164 +179 @@
          !
          ALLOCATE( qgh_trd0(jpi,jpj) )    ! allocation
+         ALLOCATE( qgh_trd1(jpi,jpj) )    ! allocation
          !
          SELECT CASE ( nn_geoflx )        ! geothermal heat flux / (rauO * Cp)
@@ -193 +209 @@
             !
          END SELECT
+         qgh_trd1(:,:) = qgh_trd0(:,:)
          !
       ELSE

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

@@ -32 +32 @@
    USE trdtra         ! trends: active tracers
    !
-   USE iom            ! IOM library               
+   USE iom            ! IOM library
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions
@@ -38 +38 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE stopack
 
    IMPLICIT NONE
@@ -67 +68 @@
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_1, ahv_bbl_1   ! diffusive bbl flux coefficients at u and v-points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM)
 
@@ -86 +88 @@
          &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
          &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
+         &      ahu_bbl_1(jpi,jpj) , ahv_bbl_1(jpi,jpj) ,                                          &
          &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
          !
@@ -195 +198 @@
       ALLOCATE(zptb(1:jpi, 1:jpj))
       !
+      ahu_bbl_1(:,:) = ahu_bbl(:,:)
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN
+          CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
+      ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) &
+         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
+      ahv_bbl_1(:,:) = ahv_bbl(:,:)
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN
+          CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
+      ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) &
+         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+      !
       DO jn = 1, kjpt                                     ! tracer loop
          !                                                ! ===========
@@ -203 +230 @@
             END DO
          END DO
-         !               
+         !
          DO jj = 2, jpjm1                                    ! Compute the trend
             DO ji = 2, jpim1
@@ -209 +236 @@
                zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
-                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
-                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
-                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
-                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
+                  &               + (   ahu_bbl_1(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
+                  &                   - ahu_bbl_1(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
+                  &                   + ahv_bbl_1(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
+                  &                   - ahv_bbl_1(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
             END DO
          END DO
@@ -415 +442 @@
                   za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
                   zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
-                  !                                                          ! 2*masked bottom density gradient 
+                  !                                                          ! 2*masked bottom density gradient
                   zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
                             - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
@@ -578 +605 @@
                gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
             ENDIF
-            !     
+            !
             IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
                mgrhv(ji,jj) = INT(  SIGN( 1.e0, &
@@ -598 +625 @@
       ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
       ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
-
 
       IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

@@ -32 +32 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -43 +44 @@
    REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile
    !                                                               !  (key_traldf_ano only)
+#if defined key_traldf_c3d
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: ahtu0, ahtv0, ahtw0, ahtt0
+#endif
+#if defined key_traldf_c2d
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:  ) :: ahtu0, ahtv0, ahtw0, ahtt0
+#endif
 
    !! * Substitutions
@@ -75 +82 @@
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
       ENDIF
+
+#if defined key_traldf_c3d
+         IF(  ( kt == nit000 ) .AND. &
+            & ( ln_stopack )   .AND. &
+            & ( ( nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt ) > 0 ) ) THEN
+             ALLOCATE ( ahtu0(jpi,jpj,jpk), ahtv0(jpi,jpj,jpk) )
+             ALLOCATE ( ahtt0(jpi,jpj,jpk), ahtw0(jpi,jpj,jpk) )
+             ahtu0 = ahtu
+             ahtv0 = ahtv
+             ahtw0 = ahtw
+             ahtt0 = ahtt
+         ENDIF
+#endif
+#if defined key_traldf_c2d
+         IF(  ( kt == nit000 ) .AND. &
+            & ( ln_stopack )   .AND. &
+            & ( ( nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt ) > 0 ) ) THEN
+             ALLOCATE ( ahtu0(jpi,jpj), ahtv0(jpi,jpj) )
+             ALLOCATE ( ahtt0(jpi,jpj), ahtw0(jpi,jpj) )
+             ahtu0 = ahtu
+             ahtv0 = ahtv
+             ahtw0 = ahtw
+             ahtt0 = ahtt
+         ENDIF
+#endif
+#if defined key_traldf_c3d || defined key_traldf_c2d
+         IF( ln_stopack .AND. ( nn_spp_ahtu > 0 ) ) THEN
+             ahtu = ahtu0
+             CALL spp_aht(kt, ahtu, nn_spp_ahtu, rn_ahtu_sd, jk_spp_ahtu)
+         ENDIF
+         IF( ln_stopack .AND. ( nn_spp_ahtv > 0 ) ) THEN
+             ahtv = ahtv0
+             CALL spp_aht(kt, ahtv, nn_spp_ahtv, rn_ahtv_sd, jk_spp_ahtv)
+         ENDIF
+         IF( ln_stopack .AND. ( nn_spp_ahtw > 0 ) ) THEN
+             ahtw = ahtw0
+             CALL spp_aht(kt, ahtw, nn_spp_ahtw, rn_ahtw_sd, jk_spp_ahtw)
+         ENDIF
+         IF( ln_stopack .AND. ( nn_spp_ahtt > 0 ) ) THEN
+             ahtt = ahtt0
+             CALL spp_aht(kt, ahtt, nn_spp_ahtt, rn_ahtt_sd, jk_spp_ahtt)
+         ENDIF
+#endif
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -9 +9 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
-   !!            3.2  !  2009-04  (G. Madec & NEMO team) 
-   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.2  !  2009-04  (G. Madec & NEMO team)
+   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
    !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
    !!----------------------------------------------------------------------
@@ -33 +33 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -42 +43 @@
    !                                 !!* Namelist namtra_qsr: penetrative solar radiation
    LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
-   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag  
+   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
@@ -50 +51 @@
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
- 
+
    ! Module variables
-   REAL(wp) ::   xsi0r                           !: inverse of rn_si0
+   REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0
    REAL(wp) ::   xsi1r                           !: inverse of rn_si1
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
@@ -79 +80 @@
       !!      Considering the 2 wavebands case:
       !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
-      !!         The temperature trend associated with the solar radiation penetration 
+      !!         The temperature trend associated with the solar radiation penetration
       !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
       !!         At the bottom, boudary condition for the radiation is no flux :
@@ -85 +86 @@
       !!      in the last ocean level.
       !!         In z-coordinate case, the computation is only done down to the
-      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients 
+      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
       !!      used for the computation are calculated one for once as they
       !!      depends on k only.
@@ -105 +106 @@
       REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
       REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
-      REAL(wp) ::   zlogc, zlogc2, zlogc3 
+      REAL(wp) ::   zlogc, zlogc2, zlogc3
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d
@@ -112 +113 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
       !
-      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
+      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )
       !
       IF( kt == nit000 ) THEN
@@ -124 +125 @@
 
       IF( l_trdtra ) THEN      ! Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdt )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
       ENDIF
@@ -153 +154 @@
       !                                        Compute now qsr tracer content field
       !                                        ************************************
-      
+
       !                                           ! ============================================== !
       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
@@ -162 +163 @@
          !                                        Add to the general trend
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   z1_e3t = zfact / fse3t(ji,jj,jk)
@@ -183 +184 @@
          ENDIF
          !                                        ! ============================================== !
-      ELSE                                        !  Ocean alone : 
+      ELSE                                        !  Ocean alone :
          !                                        ! ============================================== !
          !
-         !                                                ! ------------------------- !
+         !
+#if defined key_traldf_c2d || key_traldf_c3d
+         IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN
+            xsi0r = rn_si0
+            CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 )
+            xsi0r = 1.e0 / xsi0r
+         ENDIF
+#else
+         IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) &
+            & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+         !                                               ! ------------------------- !
          IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
             !                                             ! ------------------------- !
@@ -196 +210 @@
                   CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
+                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)
                   ENDDO
                   !
@@ -217 +231 @@
                         zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
                         zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
-                        zCze    = 1.12  * (zchl)**0.803 
+                        zCze    = 1.12  * (zchl)**0.803
                         DO jk = 1, nksr + 1
                            zpsi = fsdept(ji,jj,jk) / zze
@@ -228 +242 @@
                ELSE                              !* Variable ocean volume but constant chrlorophyll
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = 0.05 
+                     zchl3d(:,:,jk) = 0.05
                   ENDDO
                ENDIF
@@ -253 +267 @@
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
-!CDIR NOVERRCHK   
-                     DO ji = 1, jpi
-                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
+!CDIR NOVERRCHK
+                     DO ji = 1, jpi
+                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) )
                         zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
                         zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
@@ -286 +300 @@
                      END DO
                   END DO
-                  ! 
+                  !
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     )
+                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r(ji,jj) )
                         zc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
                         zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
                         zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
-                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2)
                      END DO
                   END DO
@@ -314 +328 @@
             !                                             ! ------------------------- !
             !
-            IF( lk_vvl ) THEN                                  !* variable volume
+            IF( lk_vvl .OR. ( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) ) THEN        !* variable volume
+
                zz0   =        rn_abs   * r1_rau0_rcp
                zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
-               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m 
+               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
-                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
+                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
+                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
+                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )
                      END DO
                   END DO
@@ -330 +345 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )
+                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )
+                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )
                         fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp
                      END DO
@@ -340 +355 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        ! (ISF) no light penetration below the ice shelves         
+                        ! (ISF) no light penetration below the ice shelves
                         qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1)
                      END DO
@@ -356 +371 @@
          !                                        Add to the general trend
          DO jk = 1, nksr
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   z1_e3t = zfact / fse3t(ji,jj,jk)
@@ -375 +390 @@
             IF(lflush) CALL flush(numout)
          ENDIF
-         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
+         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
-         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm 
+         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm
          IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
          !
@@ -385 +400 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )
       ENDIF
       !                       ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
-      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
+      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
@@ -407 +422 @@
       !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
       !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
-      !!      default values correspond to clear water (type I in Jerlov' 
+      !!      default values correspond to clear water (type I in Jerlov'
       !!      (1968) classification.
       !!         called by tra_qsr at the first timestep (nit000)
@@ -434 +449 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
       !
-      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
       !
 
@@ -465 +480 @@
 
       IF( ln_traqsr ) THEN     ! control consistency
-         !                      
+         !
          IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN
             CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
@@ -480 +495 @@
             &              ' 2 bands, 3 RGB bands or bio-model light penetration' )
          !
-         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
+         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1
          IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
          IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3
@@ -498 +513 @@
       ENDIF
       !                          ! ===================================== !
-      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !  
+      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !
          !                       ! ===================================== !
          !
+         ALLOCATE( xsi0r(jpi,jpj) )
          xsi0r = 1.e0 / rn_si0
          xsi1r = 1.e0 / rn_si1
@@ -539 +555 @@
                   zchl = 0.05                                 ! constant chlorophyll
                   irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll 
+                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll
                   zekg(:,:) = rkrgb(2,irgb)
                   zekr(:,:) = rkrgb(3,irgb)
@@ -546 +562 @@
                   ze0(:,:,1) = rn_abs
                   ze1(:,:,1) = zcoef
-                  ze2(:,:,1) = zcoef 
+                  ze2(:,:,1) = zcoef
                   ze3(:,:,1) = zcoef
                   zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask
-               
+
                   DO jk = 2, nksr+1
 !CDIR NOVERRCHK
                      DO jj = 1, jpj
-!CDIR NOVERRCHK   
+!CDIR NOVERRCHK
                         DO ji = 1, jpi
-                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r     )
+                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) )
                            zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) )
                            zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) )
@@ -566 +582 @@
                         END DO
                      END DO
-                  END DO 
+                  END DO
                   !
                   DO jk = 1, nksr
@@ -600 +616 @@
                   DO jj = 1, jpj                              ! top 400 meters
                      DO ji = 1, jpi
-                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
-                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
-                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
+                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
+                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
+                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)
                      END DO
                   END DO
@@ -611 +627 @@
          ENDIF
          !                       ! ===================================== !
-      ELSE                       !        No light penetration           !                   
+      ELSE                       !        No light penetration           !
          !                       ! ===================================== !
          IF(lwp) THEN
@@ -617 +633 @@
             WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
             WRITE(numout,*) '~~~~~~~~~~~~'
-            IF(lflush) CALL flush(numout)
+            IF(lflush) CALL flush(numout)            
          ENDIF
       ENDIF
@@ -630 +646 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init')

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90

@@ -29 +29 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! Memory allocation
+   USE stopack     
 
    IMPLICIT NONE
@@ -68 +69 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt )
+      IF( ln_dyn_trd .AND. ln_stopack .AND. ln_sppt_dyn ) &
+         & CALL dyn_sppt_collect( putrd, pvtrd, ktrd, kt )
          
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

@@ -32 +32 @@
    USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
+   USE stopack
    USE wrk_nemo       ! Memory allocation
 
@@ -271 +272 @@
       !                   ! 3D output of tracers trends using IOM interface
       IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt )
-
-      !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_tra_trd .AND. ln_stopack .AND. ln_sppt_tra )  &
+         & CALL tra_sppt_collect( ptrdx, ptrdy, ktrd, kt )
+
+      !  Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt )
 

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

@@ -26 +26 @@
    USE wrk_nemo        ! Memory Allocation
    USE phycst, ONLY: vkarmn
+   USE stopack
 
    IMPLICIT NONE
@@ -52 +53 @@
    REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
    LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d,bfrcoef2d0   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
 
    !! * Substitutions
@@ -68 +69 @@
       !!                ***  FUNCTION zdf_bfr_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
+      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), bfrcoef2d0(jpi,jpj),STAT=zdf_bfr_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc )
@@ -107 +108 @@
          IF(lflush) CALL flush(numout)
       ENDIF
+      !
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN
+         bfrcoef2d = bfrcoef2d0
+         CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr)
+      ENDIF
+#else
+      IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) &
+         & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
       !
       IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
@@ -135 +147 @@
                END DO
             END IF
-         !   
+         !
          ELSE
             zbfrt(:,:) = bfrcoef2d(:,:)
@@ -178 +190 @@
                DO ji = 2, jpim1
                   ! (ISF) ========================================================================
-                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
                   ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
                   !
@@ -359 +371 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
-         
+
          IF ( ln_isfcav ) THEN
             IF(ln_tfr2d) THEN
@@ -494 +506 @@
       ENDIF
       !
+      bfrcoef2d0(:,:) = bfrcoef2d(:,:)
       IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init')
       !

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

@@ -20 +20 @@
    USE zdfkpp          ! KPP vertical mixing
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE in_out_manager  ! I/O manager
    USE iom             ! for iom_put
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE timing          ! Timing
+   USE stopack
 
    IMPLICIT NONE
@@ -30 +31 @@
 
    PUBLIC   zdf_evd    ! called by step.F90
+   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: rn_avevd0
 
    !! * Substitutions
@@ -43 +45 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_evd  ***
-      !!                   
+      !!
       !! ** Purpose :   Local increased the vertical eddy viscosity and diffu-
       !!      sivity coefficients when a static instability is encountered.
       !!
       !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients
-      !!      are set to avevd (namelist parameter) if the water column is 
+      !!      are set to avevd (namelist parameter) if the water column is
       !!      statically unstable (i.e. if rn2 < -1.e-12 )
       !!
@@ -70 +72 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp .AND. lflush) CALL flush(numout)
+         ALLOCATE ( rn_avevd0(jpi,jpj) )
+         rn_avevd0(:,:) = rn_avevd
       ENDIF
 
       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application
+
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN
+         rn_avevd0(:,:) = rn_avevd
+         CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd)
+      ENDIF
+#else
+      IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) &
+         & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
 
       SELECT CASE ( nn_evdm )
@@ -80 +95 @@
          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
          !
-         DO jk = 1, jpkm1 
+         DO jk = 1, jpkm1
             DO jj = 2, jpj             ! no vector opt.
                DO ji = 2, jpi
@@ -89 +104 @@
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
 #endif
-                     avt (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk)
-                     avm (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk)
-                     avmu(ji  ,jj  ,jk) = rn_avevd * umask(ji  ,jj  ,jk)
-                     avmu(ji-1,jj  ,jk) = rn_avevd * umask(ji-1,jj  ,jk)
-                     avmv(ji  ,jj  ,jk) = rn_avevd * vmask(ji  ,jj  ,jk)
-                     avmv(ji  ,jj-1,jk) = rn_avevd * vmask(ji  ,jj-1,jk)
+                     avt (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk)
+                     avm (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk)
+                     avmu(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji  ,jj  ,jk)
+                     avmu(ji-1,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji-1,jj  ,jk)
+                     avmv(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj  ,jk)
+                     avmv(ji  ,jj-1,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj-1,jk)
                   ENDIF
                END DO
             END DO
-         END DO 
+         END DO
          CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions
          CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
@@ -105 +120 @@
          CALL iom_put( "avm_evd", zavm_evd )              ! output this change
          !
-      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12)
          DO jk = 1, jpkm1
-!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
+!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!
             DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call)
                DO ji = 1, jpi
 #if defined key_zdfkpp
                   ! no evd mixing in the boundary layer with KPP
-                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &          
+                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &
 #else
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
 #endif
-                     avt(ji,jj,jk) = rn_avevd * tmask(ji,jj,jk)
+                     avt(ji,jj,jk) = rn_avevd0(ji,jj) * tmask(ji,jj,jk)
                END DO
             END DO
          END DO
          !
-      END SELECT 
+      END SELECT
 
       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdfgls  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the gls 
+   !! Ocean physics:  vertical mixing coefficient computed from the gls
    !!                 turbulent closure parameterization
    !!======================================================================
@@ -16 +16 @@
    !!   gls_rst       : read/write gls restart in ocean restart file
    !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and active tracers 
+   USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
    USE domvvl         ! ocean space and time domain : variable volume layer
@@ -31 +31 @@
    USE iom            ! I/O manager library
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE stopack
 
    IMPLICIT NONE
@@ -61 +62 @@
    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
-   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)
 
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
-   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
-   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn    
+   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
+   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn
    REAL(wp) ::   rghmin        = -0.28_wp
    REAL(wp) ::   rgh0          =  0.0329_wp
@@ -72 +73 @@
    REAL(wp) ::   ra2           =  0.74_wp
    REAL(wp) ::   rb1           = 16.60_wp
-   REAL(wp) ::   rb2           = 10.10_wp         
-   REAL(wp) ::   re2           =  1.33_wp         
+   REAL(wp) ::   rb2           = 10.10_wp
+   REAL(wp) ::   re2           =  1.33_wp
    REAL(wp) ::   rl1           =  0.107_wp
    REAL(wp) ::   rl2           =  0.0032_wp
@@ -133 +134 @@
       INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
       REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
-      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      - 
+      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -
       REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      -
       REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
@@ -139 +140 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
@@ -156 +157 @@
       CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
       CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
-      
+
       ! Preliminary computing
 
@@ -165 +166 @@
          avm (:,:,:) = avm_k (:,:,:)
          avmu(:,:,:) = avmu_k(:,:,:)
-         avmv(:,:,:) = avmv_k(:,:,:) 
+         avmv(:,:,:) = avmv_k(:,:,:)
       ENDIF
 
       ! Compute surface and bottom friction at T-points
-!CDIR NOVERRCHK          
-      DO jj = 2, jpjm1          
-!CDIR NOVERRCHK         
-         DO ji = fs_2, fs_jpim1   ! vector opt.         
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
             !
             ! surface friction
             ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
-            !   
-            ! bottom friction (explicit before friction)        
-            ! Note that we chose here not to bound the friction as in dynbfr)   
-            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
-               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )      
-            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
-               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )      
-            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
-         END DO         
-      END DO    
+            !
+            ! bottom friction (explicit before friction)
+            ! Note that we chose here not to bound the friction as in dynbfr)
+            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &
+               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )
+            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &
+               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )
+            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+         END DO
+      END DO
 
       ! Set surface roughness length
       SELECT CASE ( nn_z0_met )
       !
-      CASE ( 0 )             ! Constant roughness          
+      CASE ( 0 )             ! Constant roughness
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
@@ -226 +227 @@
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)
@@ -275 +276 @@
                IF( ln_sigpsi ) THEN
                   zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
-                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
+                  zwall_psi(ji,jj,jk) = rsc_psi /   &
                      &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
                ELSE
@@ -294 +295 @@
                ! diagonal
                z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
+                  &                       + rdt * zdiss * tmask(ji,jj,jk)
                !
                ! right hand side in en
@@ -316 +317 @@
       ! First level
       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
+      en(:,:,1) = MAX(en(:,:,1), rn_emin)
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0._wp
       z_elem_b(:,:,1) = 1._wp
-      ! 
+      !
       ! One level below
       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) &
           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
       en(:,:,2) = MAX(en(:,:,2), rn_emin )
-      z_elem_a(:,:,2) = 0._wp 
+      z_elem_a(:,:,2) = 0._wp
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
@@ -334 +335 @@
       ! Dirichlet conditions at k=1
       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
+      en(:,:,1)       = MAX(en(:,:,1), rn_emin)
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0._wp
@@ -357 +358 @@
       SELECT CASE ( nn_bc_bot )
       !
-      CASE ( 0 )             ! Dirichlet 
+      CASE ( 0 )             ! Dirichlet
          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
@@ -377 +378 @@
                z_elem_c(ji,jj,ibotm1) = 0._wp
                z_elem_b(ji,jj,ibotm1) = 1._wp
-               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
+               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
             END DO
          END DO
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -428 +429 @@
          END DO
       END DO
-      !                                            ! set the minimum value of tke 
+      !                                            ! set the minimum value of tke
       en(:,:,:) = MAX( en(:,:,:), rn_emin )
 
@@ -470 +471 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
+                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn
                END DO
             END DO
@@ -489 +490 @@
                !
                ! psi / k
-               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
+               zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
                !
                ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
@@ -509 +510 @@
                zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
                zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
-               !                                                        
+               !
                ! building the matrix
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
@@ -551 +552 @@
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
-      ! 
+      !
       !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
@@ -575 +576 @@
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
 
-      !   
+      !
       !
       END SELECT
@@ -585 +586 @@
       !
       !
-      CASE ( 0 )             ! Dirichlet 
+      CASE ( 0 )             ! Dirichlet
          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
@@ -610 +611 @@
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -692 +693 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
+                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
                END DO
             END DO
@@ -719 +720 @@
                eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
                mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
-               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
+               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
             END DO
          END DO
-      END DO 
+      END DO
 
       !
@@ -806 +807 @@
             END DO
          END DO
+#if defined key_traldf_c2d || key_traldf_c3d
+         IF( ln_stopack) THEN
+            IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk)
+            IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk)
+         ENDIF
+#else
+         IF ( ln_stopack ) &
+            & CALL ctl_stop( 'zdf_gls: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
       END DO
       !
@@ -846 +857 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_gls_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
       !!      viscosity when using a gls turbulent closure scheme
       !!
@@ -881 +892 @@
       READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp )
-      IF(lwm .AND. nprint > 2) WRITE ( numond, namzdf_gls )
+      IF(lwm) WRITE ( numond, namzdf_gls )
 
       IF(lwp) THEN                     !* Control print
@@ -1069 +1080 @@
       END SELECT
       !
-      IF(lwp .AND. lflush) CALL flush(numout) 
+      IF(lwp .AND. lflush) CALL flush(numout)
       !                                !* Set Schmidt number for psi diffusion in the wave breaking case
       !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
       !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
       IF( ln_sigpsi ) THEN
-         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 
+         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
          ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
          ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
@@ -1082 +1093 @@
          rsc_psi0 = rsc_psi
       ENDIF
- 
+
       !                                !* Shear free turbulence parameters
       !
@@ -1129 +1140 @@
       rc04  = rc03 * rc0
       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
-      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
+      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
       rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
-      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
+      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
-      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
+      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
 
       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
@@ -1150 +1161 @@
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-      !                              
+      !
       CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
       !
@@ -1161 +1172 @@
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE ts_rst  ***
-      !!                     
+      !!
       !! ** Purpose :   Read or write TKE file (en) in restart file
       !!
       !! ** Method  :   use of IOM library
-      !!                if the restart does not contain TKE, en is either 
+      !!                if the restart does not contain TKE, en is either
       !!                set to rn_emin or recomputed (nn_igls/=0)
       !!----------------------------------------------------------------------
@@ -1177 +1188 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
          !                                   ! ---------------
          IF( ln_rstart ) THEN                   !* Read the restart file
@@ -1196 +1207 @@
                CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
                IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
-            ELSE                        
+            ELSE
                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
                IF(lwp .AND. lflush) CALL flush(numout)
                en  (:,:,:) = rn_emin
-               mxln(:,:,:) = 0.05        
+               mxln(:,:,:) = 0.05
                avt_k (:,:,:) = avt (:,:,:)
                avm_k (:,:,:) = avm (:,:,:)
@@ -1211 +1222 @@
             IF(lwp .AND. lflush) CALL flush(numout)
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            mxln(:,:,:) = 0.05
          ENDIF
          !
@@ -1219 +1230 @@
          IF(lwp .AND. lflush) CALL flush(numout)
          IF(nn_timing == 2)  CALL timing_start('iom_rstput')
-         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
          CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
          CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
-         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
          CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
          CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )
@@ -1241 +1252 @@
       WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
    END SUBROUTINE zdf_gls_init
-   
+
    SUBROUTINE zdf_gls( kt )          ! Empty routine
    IMPLICIT NONE
-      INTEGER, INTENT(in) ::   kt ! ocean time step   
+      INTEGER, INTENT(in) ::   kt ! ocean time step
       WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_gls
-   
+
    SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
    IMPLICIT NONE
@@ -1258 +1269 @@
    !!======================================================================
 END MODULE zdfgls
-

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdftke  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the tke 
+   !! Ocean physics:  vertical mixing coefficient computed from the tke
    !!                 turbulent closure parameterization
    !!=====================================================================
@@ -22 +22 @@
    !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
    !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
-   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
-   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
+   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
+   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
@@ -52 +52 @@
    USE wrk_nemo       ! work arrays
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 #if defined key_agrif
    USE agrif_opa_interp
    USE agrif_opa_update
 #endif
-
-
+   USE stopack
 
    IMPLICIT NONE
@@ -75 +74 @@
    INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
    REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
-   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
+   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
    REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
    REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
@@ -94 +93 @@
 
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
 #if defined key_c1d
@@ -102 +99 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
 #endif
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
 
    !! * Substitutions
@@ -118 +118 @@
       !!----------------------------------------------------------------------
       ALLOCATE(                                                                    &
-         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &      
+         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &
 #if defined key_c1d
          &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
@@ -147 +147 @@
       !!         surface: en = max( rn_emin0, rn_ebb * taum )
       !!         bottom : en = rn_emin
-      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
-      !!
-      !!        The now Turbulent kinetic energy is computed using the following 
+      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
+      !!
+      !!        The now Turbulent kinetic energy is computed using the following
       !!      time stepping: implicit for vertical diffusion term, linearized semi
-      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
-      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
+      !!      implicit for kolmogoroff dissipation term, and explicit forward for
+      !!      both buoyancy and shear production terms. Therefore a tridiagonal
       !!      linear system is solved. Note that buoyancy and shear terms are
       !!      discretized in a energy conserving form (Bruchard 2002).
@@ -160 +160 @@
       !!
       !!        The now vertical eddy vicosity and diffusivity coefficients are
-      !!      given by: 
+      !!      given by:
       !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!              avt = max( avmb, pdl * avm                 )  
+      !!              avt = max( avmb, pdl * avm                 )
       !!              eav = max( avmb, avm )
       !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
-      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
+      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
       !!
       !! ** Action  :   compute en (now turbulent kinetic energy)
@@ -180 +180 @@
       !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:) 
-         avm (:,:,:) = avm_k (:,:,:) 
-         avmu(:,:,:) = avmu_k(:,:,:) 
-         avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF 
+         avt (:,:,:) = avt_k (:,:,:)
+         avm (:,:,:) = avm_k (:,:,:)
+         avmu(:,:,:) = avmu_k(:,:,:)
+         avmv(:,:,:) = avmv_k(:,:,:)
+      ENDIF
+      !
+#if defined key_traldf_c2d || key_traldf_c3d
+      IF( ln_stopack) THEN
+         IF( nn_spp_tkelc > 0 ) THEN
+             rn_lc0 = rn_lc
+             CALL spp_gen( kt, rn_lc0, nn_spp_tkelc, rn_tkelc_sd, jk_spp_tkelc )
+         ENDIF
+         IF( nn_spp_tkedf > 0 ) THEN
+             rn_ediff0 = rn_ediff
+             CALL spp_gen( kt, rn_ediff0, nn_spp_tkedf, rn_tkedf_sd, jk_spp_tkedf )
+         ENDIF
+         IF( nn_spp_tkeds > 0 ) THEN
+             rn_ediss0 = rn_ediss
+             CALL spp_gen( kt, rn_ediss0, nn_spp_tkeds, rn_tkeds_sd, jk_spp_tkeds )
+         ENDIF
+         IF( nn_spp_tkebb > 0 ) THEN
+             rn_ebb0 = rn_ebb
+             CALL spp_gen( kt, rn_ebb0, nn_spp_tkebb, rn_tkebb_sd, jk_spp_tkebb )
+        ENDIF
+         IF( nn_spp_tkefr > 0 ) THEN
+             rn_efr0 = rn_efr
+             CALL spp_gen( kt, rn_efr0, nn_spp_tkefr, rn_tkefr_sd, jk_spp_tkefr )
+         ENDIF
+      ENDIF
+#else
+      IF ( ln_stopack ) &
+         & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
       !
       CALL tke_tke      ! now tke (en)
@@ -190 +219 @@
       CALL tke_avn      ! now avt, avm, avmu, avmv
       !
-      avt_k (:,:,:) = avt (:,:,:) 
-      avm_k (:,:,:) = avm (:,:,:) 
-      avmu_k(:,:,:) = avmu(:,:,:) 
-      avmv_k(:,:,:) = avmv(:,:,:) 
+      avt_k (:,:,:) = avt (:,:,:)
+      avm_k (:,:,:) = avm (:,:,:)
+      avmu_k(:,:,:) = avmu(:,:,:)
+      avmv_k(:,:,:) = avmv(:,:,:)
       !
 #if defined key_agrif
-      ! Update child grid f => parent grid 
+      ! Update child grid f => parent grid
       IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
-#endif      
-     ! 
+#endif
+      IF (  kt == nitend ) THEN
+        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 )
+      ENDIF
+      !
+
    END SUBROUTINE zdf_tke
 
@@ -225 +258 @@
       REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
-      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
+      REAL(wp) ::   zesh2                           ! temporary scalars
+      REAL(wp) ::   zfact1                          !    -         -
       REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
       REAL(wp) ::   ztau  , zdif                    !    -         -
@@ -233 +266 @@
 !!bfr      REAL(wp) ::   zebot                           !    -         -
       INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc, zbbrau,zfact2,zfact3
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
@@ -240 +273 @@
       !
       CALL wrk_alloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_alloc( jpi,jpj, zhlc ) 
-      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
-      !
-      zbbrau = rn_ebb / rau0       ! Local constant initialisation
-      zfact1 = -.5_wp * rdt 
-      zfact2 = 1.5_wp * rdt * rn_ediss
-      zfact3 = 0.5_wp       * rn_ediss
+      CALL wrk_alloc( jpi,jpj, zhlc )
+      CALL wrk_alloc( jpi,jpj, zbbrau )
+      CALL wrk_alloc( jpi,jpj, zfact2 )
+      CALL wrk_alloc( jpi,jpj, zfact3 )
+      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
+      !
+      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation
+      zfact1 = -.5_wp * rdt
+      zfact2 = 1.5_wp * rdt * rn_ediss0
+      zfact3 = 0.5_wp       * rn_ediss0
       !
       !
@@ -261 +297 @@
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      
+            en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1)
+         END DO
+      END DO
+
 !!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -303 +339 @@
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
          DO jk = jpkm1, 2, -1
-            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
+            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
                DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
                   zus  = zcof * taum(ji,jj)
@@ -311 +347 @@
          END DO
          !                               ! finite LC depth
-         DO jj = 1, jpj 
+         DO jj = 1, jpj
             DO ji = 1, jpi
                zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
@@ -326 +362 @@
                   !                                           ! vertical velocity due to LC
                   zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
-                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   &
+                  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)
+
                END DO
             END DO
@@ -347 +384 @@
             DO ji = 1, jpi
                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
-                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
+                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
                   &                            / (  fse3uw_n(ji,jj,jk)               &
                   &                              *  fse3uw_b(ji,jj,jk)  )
@@ -376 +413 @@
                   !                                                           ! shear prod. at w-point weightened by mask
                zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
                   !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
                zd_lw(ji,jj,jk) = zzd_lw
-               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
+               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk)
                !
                !                                   ! right hand side in en
                en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
-                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
+                  &                                 + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
                   &                                 * wmask(ji,jj,jk)
             END DO
@@ -433 +470 @@
       END DO
 
-      !                                 ! Save TKE prior to nn_etau addition  
-      e_niw(:,:,:) = en(:,:,:)  
-      !  
+      !                                 ! Save TKE prior to nn_etau addition
+      e_niw(:,:,:) = en(:,:,:)
+      !
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  TKE due to surface and internal wave breaking
@@ -455 +492 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
                      &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -464 +501 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                jk = nmln(ji,jj)
-               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
                   &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
@@ -477 +514 @@
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
-                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
+                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
-                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
                      &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -486 +523 @@
          END DO
       ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up)
-         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
+         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -504 +541 @@
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
-      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term  
-         DO jj = 2, jpjm1  
-            DO ji = fs_2, fs_jpim1   ! vector opt.  
-               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)  
-            END DO  
-         END DO  
-      END DO  
-      !  
-      CALL lbc_lnk( e_niw, 'W', 1. )  
+      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( e_niw, 'W', 1. )
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_dealloc( jpi,jpj, zhlc ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_dealloc( jpi,jpj, zhlc )
+      CALL wrk_dealloc( jpi,jpj, zbbrau )
+      CALL wrk_dealloc( jpi,jpj, zfact2 )
+      CALL wrk_dealloc( jpi,jpj, zfact3 )
+      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
@@ -529 +569 @@
       !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
       !!
-      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
-      !!              the tke_tke routine). First, the now mixing lenth is 
+      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
+      !!              the tke_tke routine). First, the now mixing lenth is
       !!      computed from en and the strafification (N^2), then the mixings
       !!      coefficients are computed.
       !!              - Mixing length : a first evaluation of the mixing lengh
       !!      scales is:
-      !!                      mxl = sqrt(2*en) / N  
+      !!                      mxl = sqrt(2*en) / N
       !!      where N is the brunt-vaisala frequency, with a minimum value set
       !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
-      !!        The mixing and dissipative length scale are bound as follow : 
+      !!        The mixing and dissipative length scale are bound as follow :
       !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
       !!                        zmxld = zmxlm = mxl
       !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
-      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
       !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
       !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
-      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
-      !!                    the surface to obtain ldown. the resulting length 
+      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
+      !!                    the surface to obtain ldown. the resulting length
       !!                    scales are:
-      !!                        zmxld = sqrt( lup * ldown ) 
+      !!                        zmxld = sqrt( lup * ldown )
       !!                        zmxlm = min ( lup , ldown )
       !!              - Vertical eddy viscosity and diffusivity:
       !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!                      avt = max( avmb, pdlr * avm )  
+      !!                      avt = max( avmb, pdlr * avm )
       !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
       !!
@@ -567 +607 @@
       IF( nn_timing == 1 )  CALL timing_start('tke_avn')
 
-      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
 
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -576 +616 @@
       !
       ! initialisation of interior minimum value (avoid a 2d loop with mikt)
-      zmxlm(:,:,:)  = rmxl_min    
+      zmxlm(:,:,:)  = rmxl_min
       zmxld(:,:,:)  = rmxl_min
       !
@@ -586 +626 @@
             END DO
          END DO
-      ELSE 
+      ELSE
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
@@ -604 +644 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
+      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value
       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
       !
@@ -698 +738 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zsqen = SQRT( en(ji,jj,jk) )
-               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
+               zav   = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen
                avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
@@ -704 +744 @@
             END DO
          END DO
+#if defined key_traldf_c2d || key_traldf_c3d
+         IF( ln_stopack) THEN
+            IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk)
+            IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk)
+         ENDIF
+#else
+         IF ( ln_stopack  ) &
+            & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
       END DO
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
@@ -749 +799 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
@@ -759 +809 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_tke_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
       !!              viscosity when using a tke turbulent closure scheme
       !!
@@ -776 +826 @@
          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
          &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
-         &                 nn_etau , nn_htau  , rn_efr , rn_c   
+         &                 nn_etau , nn_htau  , rn_efr , rn_c
       !!----------------------------------------------------------------------
 
@@ -843 +893 @@
          rn_mxl0 = rmxl_min
       ENDIF
-      
+
+      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc
+      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff
+      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss
+      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb
+      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr
+
       IF( nn_etau == 2  ) THEN
           ierr = zdf_mxl_alloc()
@@ -855 +911 @@
 
       !                               !* depth of penetration of surface tke
-      IF( nn_etau /= 0 ) THEN      
+      IF( nn_etau /= 0 ) THEN
          htau(:,:) = 0._wp
          SELECT CASE( nn_htau )             ! Choice of the depth of penetration
@@ -861 +917 @@
             htau(:,:) = 10._wp
          CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
-            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
          CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer
             rhtau = -1._wp / LOG( 1._wp - rn_c )
@@ -904 +960 @@
       END DO
       dissl(:,:,:) = 1.e-12_wp
-      !                              
+      !
       CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
       !
@@ -913 +969 @@
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE tke_rst  ***
-     !!                     
+     !!
      !! ** Purpose :   Read or write TKE file (en) in restart file
      !!
      !! ** Method  :   use of IOM library
-     !!                if the restart does not contain TKE, en is either 
-     !!                set to rn_emin or recomputed 
+     !!                if the restart does not contain TKE, en is either
+     !!                set to rn_emin or recomputed
      !!----------------------------------------------------------------------
      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
@@ -927 +983 @@
      !!----------------------------------------------------------------------
      !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
         !                                   ! ---------------
         IF( ln_rstart ) THEN                   !* Read the restart file
@@ -956 +1012 @@
               CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
               !
-              avt_k (:,:,:) = avt (:,:,:)
-              avm_k (:,:,:) = avm (:,:,:)
-              avmu_k(:,:,:) = avmu(:,:,:)
-              avmv_k(:,:,:) = avmv(:,:,:)
-              !
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
            ENDIF
         ELSE                                   !* Start from rest
            en(:,:,:) = rn_emin * tmask(:,:,:)
-           DO jk = 1, jpk                           ! set the Kz to the background value
-              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
-           END DO
         ENDIF
         !
+        avt_k (:,:,:) = avt (:,:,:)
+        avm_k (:,:,:) = avm (:,:,:)
+        avmu_k(:,:,:) = avmu(:,:,:)
+        avmv_k(:,:,:) = avmv(:,:,:)
+
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -85 +85 @@
    USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 
    USE sbc_oce, ONLY: lk_oasis
+   USE stopack
    USE stopar
    USE stopts
@@ -498 +499 @@
                             CALL dia_hsb_init   ! heat content, salt content and volume budgets
                             CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends
-                            CALL     bias_init  ! Pressure correction bias
+                             CALL stopack_init   ! STOPACK scheme                            
+                            CALL    bias_init   ! Pressure correction bias
       IF( lk_diaobs     ) THEN                  ! Observation & model comparison
                             CALL dia_obs_init            ! Initialize observational data

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -110 +110 @@
       IF( lk_tide    )   CALL sbc_tide( kstp )
       IF( lk_bdy     )  THEN
-         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
+         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib
                          CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
       ENDIF
@@ -119 +119 @@
       END DO
 
+      IF( ln_stopack )   CALL stopack_pert( kstp )
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
+
+      IF( ln_sto_eos )                 CALL sto_par( kstp )   ! Stochastic parameters
+      IF( ln_sto_eos )                 CALL sto_pts( tsn  )   ! Random T/S fluctuations
+      IF( ln_stopack .AND. ln_skeb  )  CALL skeb_comp( kstp ) ! Stochastic Energy Backscatter
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -149 +153 @@
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
+#if defined key_traldf_c2d || key_traldf_c3d
+         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) THEN
+              rn_avt_rnf0 = rn_avt_rnf
+              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf )
+         ENDIF
+#else
+         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) &
+            & CALL ctl_stop( 'stp: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
+         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
       ENDIF
       IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
@@ -163 +177 @@
       IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
       IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
+      IF( lrst_oce .AND. ln_stopack)   CALL stopack_rst( kstp, 'WRITE' )
       !
       !  LATERAL  PHYSICS
@@ -195 +210 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
-      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
-                         CALL wzv           ( kstp )  ! now cross-level velocity 
-
-      IF( lk_dynspg_ts ) THEN 
+      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
+                         CALL wzv           ( kstp )  ! now cross-level velocity
+
+      IF( lk_dynspg_ts ) THEN
           ! In case the time splitting case, update almost all momentum trends here:
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
+            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                               &
@@ -214 +230 @@
                                   va(:,:,:) = 0.e0
           IF(  lk_asminc .AND. ln_asmiau .AND. &
-             & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment
-          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified)
-          IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends
-                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form)
-                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis
-                                  CALL dyn_ldf      ( kstp )   ! lateral mixing
-          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
+             & ln_dyninc       )  CALL dyn_asm_inc   ( kstp )   ! apply dynamics assimilation increment
+          IF( ln_stopack .AND. ln_sppt_dyn ) &
+             &                    CALL dyn_sppt_apply( kstp, 0 )
+          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! subtract Neptune velocities (simplified)
+          IF( lk_bdy           )  CALL bdy_dyn3d_dmp ( kstp )   ! bdy damping trends
+                                  CALL dyn_adv       ( kstp )   ! advection (vector or flux form)
+                                  CALL dyn_vor       ( kstp )   ! vorticity term including Coriolis
+                                  CALL dyn_ldf       ( kstp )   ! lateral mixing
+          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! add Neptune velocities (simplified)
 #if defined key_agrif
           IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
@@ -232 +250 @@
                                   CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
           IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
-                                  CALL wzv           ( kstp )  ! now cross-level velocity 
+                                  CALL wzv           ( kstp )  ! now cross-level velocity
       ENDIF
 
@@ -262 +280 @@
                              tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
 
-      IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
+      IF( ln_stopack .AND. ln_sppt_tra ) &
+         &                   CALL tra_sppt_apply ( kstp, 0 )
+      IF( lk_asminc .AND. ln_asmiau .AND. &
+        & ln_trainc      )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
@@ -280 +300 @@
       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
 #endif
+      IF( ln_stopack .AND. ln_sppt_tra ) &
+         &                   CALL tra_sppt_apply ( kstp, 1 )
                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
 
@@ -285 +307 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+         IF( ln_stopack .AND. ln_sppt_tra ) &
+            &                CALL tra_sppt_apply ( kstp, 2 )
                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                                &
@@ -300 +324 @@
                &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
                &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
-         IF( ln_zps .AND.       ln_isfcav)                                   & 
+         IF( ln_zps .AND.       ln_isfcav)                                   &
                &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
@@ -308 +332 @@
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
          IF( ln_bias )       CALL dyn_bias( kstp )
+         IF( ln_stopack .AND. ln_sppt_tra )   THEN
+                                CALL tra_sppt_apply ( kstp, 2 )
+                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
+                                IF( ln_zps .AND. .NOT. ln_isfcav)                &
+                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+                                IF( ln_zps .AND.       ln_isfcav)                                   &
+                                & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+                                &                   rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+                                &                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+         ENDIF
       ENDIF
 
@@ -318 +353 @@
                                ua(:,:,:) = ua_sv(:,:,:)
                                va(:,:,:) = va_sv(:,:,:)
-                                                             ! Revert now divergence and rotational to previously computed ones 
+                                                             ! Revert now divergence and rotational to previously computed ones
                                                              !(needed because of the time swap in div_cur, at the beginning of each time step)
                                hdivn(:,:,:) = hdivb(:,:,:)
-                               rotn(:,:,:)  = rotb(:,:,:) 
+                               rotn(:,:,:)  = rotb(:,:,:)
 
                                CALL dyn_bfr( kstp )         ! bottom friction
+            IF( ln_stopack .AND. ln_sppt_dyn ) &
+               &               CALL dyn_sppt_apply ( kstp, 1 )
                                CALL dyn_zdf( kstp )         ! vertical diffusion
       ELSE
@@ -331 +368 @@
         IF(  lk_asminc .AND. ln_asmiau .AND. &
            & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
+        IF( ln_stopack .AND. ln_sppt_dyn ) &
+           &                   CALL dyn_sppt_apply ( kstp, 0 )
         IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
@@ -343 +382 @@
                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
                                CALL dyn_bfr( kstp )         ! bottom friction
+            IF( ln_stopack .AND. ln_sppt_dyn ) &
+               &               CALL dyn_sppt_apply ( kstp, 1 )
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient
       ENDIF
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step
+            IF( ln_stopack .AND. ln_sppt_dyn ) &
+               &               CALL dyn_sppt_apply ( kstp, 2 )
+            IF( ln_stopack .AND. ln_skeb ) &
+               &               CALL skeb_apply ( kstp )
 
                                CALL ssh_swp( kstp )         ! swap of sea surface height
@@ -352 +397 @@
       !
 
-      ! Moved bias _wrt to before rst_write as used restart parameters for nn_stocklist option that are changed by rst_wrt   
+      ! Moved bias write to before rst_write
       IF( lrst_bias )          CALL bias_wrt     ( kstp )
-
+      
       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
       IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
@@ -361 +406 @@
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! AGRIF
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
-                               CALL Agrif_Integrate_ChildGrids( stp )  
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                               CALL Agrif_Integrate_ChildGrids( stp )
 
       IF ( Agrif_NbStepint().EQ.0 ) THEN
@@ -393 +438 @@
       !
 #if defined key_iomput
-      IF( kstp == nitend .OR. indic < 0 ) THEN 
+      IF( kstp == nitend .OR. indic < 0 ) THEN
                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
-         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
       ENDIF
 #endif
       !
       IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
-      !     
+      !
       !
    END SUBROUTINE stp

branches/UKMO/dev_r5518_GO6_stophys/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

@@ -96 +96 @@
    USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
    USE diaharm
+   USE diaprod          ! ocean model: product diagnostics
    USE flo_oce          ! floats variables
    USE floats           ! floats computation               (flo_stp routine)
@@ -113 +114 @@
    USE timing           ! Timing
 
+   USE stopack          ! Stochastic physics
+
 #if defined key_agrif
    USE agrif_opa_sponge ! Momemtum and tracers sponges