Changeset 11384


Ignore:
Timestamp:
2019-07-31T18:05:50+02:00 (14 months ago)
Author:
mattmartin
Message:

Included Andrea Storto's STOPACK code into NEMO3.6 branch.

Location:
branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
18 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r4990 r11384  
    3030   USE wrk_nemo        ! Memory Allocation 
    3131   USE timing          ! Timing 
     32   USE stopack 
    3233 
    3334   IMPLICIT NONE 
     
    3839 
    3940   INTEGER ::   nldf = -2   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals) 
     41 
     42#if defined key_dynldf_c3d 
     43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahm10,ahm20,ahm30,ahm40 
     44#endif 
     45#if defined key_dynldf_c2d 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: ahm10,ahm20,ahm30,ahm40 
     47#endif 
    4048 
    4149   !! * Substitutions 
     
    6775         ztrdv(:,:,:) = va(:,:,:)  
    6876      ENDIF 
     77 
     78#if defined key_dynldf_c3d 
     79         IF( kt .eq. nit000 .and. (nn_spp_ahm1+nn_spp_ahm2) .gt. 0 ) THEN 
     80             ALLOCATE ( ahm10(jpi,jpj,jpk), ahm20(jpi,jpj,jpk) ) 
     81             ALLOCATE ( ahm30(jpi,jpj,jpk), ahm40(jpi,jpj,jpk) ) 
     82             ahm10 = ahm1 
     83             ahm20 = ahm2 
     84             ahm30 = ahm3 
     85             ahm40 = ahm4 
     86         ENDIF 
     87#endif 
     88#if defined key_dynldf_c2d 
     89         IF( kt .eq. nit000 .and. (nn_spp_ahm1+nn_spp_ahm2) .gt. 0 ) THEN 
     90             ALLOCATE ( ahm10(jpi,jpj), ahm20(jpi,jpj) ) 
     91             ALLOCATE ( ahm30(jpi,jpj), ahm40(jpi,jpj) ) 
     92             ahm10 = ahm1 
     93             ahm20 = ahm2 
     94             ahm30 = ahm3 
     95             ahm40 = ahm4 
     96         ENDIF 
     97#endif 
     98 
     99#if defined key_traldf_c3d || defined key_traldf_c2d 
     100         IF( nn_spp_ahm1 .GT. 0) THEN 
     101           IF( ln_dynldf_lap ) THEN 
     102             ahm1 = ahm10 
     103             CALL spp_ahm(kt,ahm1,nn_spp_ahm1,rn_ahm1_sd,jk_spp_ahm1) 
     104           ENDIF 
     105           IF( ln_dynldf_bilap ) THEN 
     106             ahm3 = ahm30 
     107             CALL spp_ahm(kt,ahm3,nn_spp_ahm1,rn_ahm1_sd,jk_spp_ahm3) 
     108           ENDIF 
     109         ENDIF 
     110         IF( nn_spp_ahm2 .GT. 0) THEN 
     111           IF( ln_dynldf_lap ) THEN 
     112             ahm2 = ahm20 
     113             CALL spp_ahm(kt,ahm2,nn_spp_ahm2,rn_ahm2_sd,jk_spp_ahm2) 
     114           ENDIF 
     115           IF( ln_dynldf_bilap ) THEN 
     116             ahm4 = ahm40 
     117             CALL spp_ahm(kt,ahm4,nn_spp_ahm2,rn_ahm2_sd,jk_spp_ahm4) 
     118           ENDIF 
     119         ENDIF 
     120#endif 
    69121 
    70122      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r4624 r11384  
    2121   USE wrk_nemo       ! work arrays 
    2222   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     23   USE stopack 
    2324 
    2425   IMPLICIT NONE 
     
    161162            END DO 
    162163         END DO 
     164 
     165         IF ( nn_spp_icealb > 0 ) CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 
     166 
    163167      END DO 
    164168       
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5487 r11384  
    5151   USE par_ice_2 
    5252#endif 
     53   USE stopack 
    5354 
    5455   IMPLICIT NONE 
     
    8990   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem) 
    9091   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
     92   REAL(wp), ALLOCATABLE, SAVE ::   rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
    9193   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9294   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     
    196198         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
    197199         ! 
     200         ALLOCATE ( rn_vfac0(jpi,jpj) ) 
     201         rn_vfac0(:,:) = rn_vfac 
     202         ! 
     203      ENDIF 
     204 
     205      IF( nn_spp_relw > 0 ) THEN 
     206         rn_vfac0 = rn_vfac 
     207         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 
    198208      ENDIF 
    199209 
     
    287297      DO jj = 2, jpjm1 
    288298         DO ji = fs_2, fs_jpim1   ! vect. opt. 
    289             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    290             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     299            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     300            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    291301         END DO 
    292302      END DO 
     
    459469            DO ji = 2, jpim1   ! B grid : NO vector opt 
    460470               ! ... scalar wind at I-point (fld being at T-point) 
    461                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    462                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    463                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    464                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
     471               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)    & 
     472                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) & 
     473                  &      - rn_vfac0(ji,jj) * u_ice(ji,jj) 
     474               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)    & 
     475                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) & 
     476                  &      - rn_vfac0(ji,jj) * v_ice(ji,jj) 
    465477               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    466478               ! ... ice stress at I-point 
     
    468480               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    469481               ! ... scalar wind at T-point (fld being at T-point) 
    470                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
    471                   &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
    472                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    473                   &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     482               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1)                                         & 
     483                  &      - rn_vfac0(ji,jj) * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     484                  &                                  + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     485               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1)                                         & 
     486                  &      - rn_vfac0(ji,jj) * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     487                  &                                  + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    474488               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    475489            END DO 
     
    482496         DO jj = 2, jpj 
    483497            DO ji = fs_2, jpi   ! vect. opt. 
    484                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    485                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     498               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     499               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
    486500               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    487501            END DO 
     
    489503         DO jj = 2, jpjm1 
    490504            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    491                utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    492                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    493                vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    494                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
     505               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )             & 
     506                  &              * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) & 
     507                  &                  - rn_vfac0(ji,jj) * u_ice(ji,jj) ) 
     508               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )             & 
     509                  &              * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) & 
     510                  &                  - rn_vfac0(ji,jj) * v_ice(ji,jj) ) 
    495511            END DO 
    496512         END DO 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r5503 r11384  
    5252   REAL(wp)                   ::   rn_hrnf         !: runoffs, depth over which enhanced vertical mixing is used 
    5353   REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s] 
    54    REAL(wp)                   ::   rn_rfact        !: multiplicative factor for runoff 
     54   REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff 
     55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_avt_rnf0 
    5556 
    5657   LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis 
     
    462463         !                                      !    - mixed upstream-centered (ln_traadv_cen2=T) 
    463464         ! 
     465         ALLOCATE ( rn_avt_rnf0(jpi,jpj) ) 
     466         rn_avt_rnf0(:,:) = rn_avt_rnf 
     467         ! 
    464468         IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   & 
    465469            &                                              'be spread through depth by ln_rnf_depth'               ) 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r4990 r11384  
    2525   USE timing         ! Timing 
    2626   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     27   USE stopack 
     28   USE wrk_nemo       ! Memory Allocation 
    2729 
    2830   IMPLICIT NONE 
     
    7577      REAL(wp) ::   zerp     ! local scalar for evaporation damping 
    7678      REAL(wp) ::   zqrp     ! local scalar for heat flux damping 
    77       REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor 
    7879      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor 
     80      REAL(wp), POINTER, DIMENSION(:,:) :: rn_dqdt_s, zsrp 
    7981      INTEGER  ::   ierror   ! return error code 
    8082      !! 
     
    9597            ! 
    9698            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term 
     99 
     100               CALL wrk_alloc( jpi, jpj, rn_dqdt_s) 
     101               rn_dqdt_s=rn_dqdt 
     102 
     103               IF( nn_spp_dqdt > 0 ) CALL spp_gen(kt, rn_dqdt_s,nn_spp_dqdt,rn_dqdt_sd,jk_spp_dqdt ) 
    97104               DO jj = 1, jpj 
    98105                  DO ji = 1, jpi 
    99                      zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
     106                     zqrp = rn_dqdt_s(ji,jj) * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
    100107                     qns(ji,jj) = qns(ji,jj) + zqrp 
    101108                     qrp(ji,jj) = zqrp 
     
    103110               END DO 
    104111               CALL iom_put( "qrp", qrp )                             ! heat flux damping 
     112               CALL wrk_dealloc( jpi, jpj, rn_dqdt_s ) 
    105113            ENDIF 
    106114            ! 
    107115            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx)) 
    108                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
     116               CALL wrk_alloc( jpi, jpj, zsrp) 
     117               zsrp = rn_deds 
     118               IF( nn_spp_dedt > 0 ) CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 
    109119!CDIR COLLAPSE 
    110120               DO jj = 1, jpj 
    111121                  DO ji = 1, jpi 
    112                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
     122                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    113123                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )  
    114124                     sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux 
     
    117127               END DO 
    118128               CALL iom_put( "erp", erp )                             ! freshwater flux damping 
     129               CALL wrk_dealloc( jpi,jpj, zsrp ) 
    119130               ! 
    120131            ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns) 
    121                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
     132               CALL wrk_alloc( jpi, jpj, zsrp) 
     133               zsrp = rn_deds 
     134               IF( nn_spp_dedt > 0 ) CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 
    122135               zerp_bnd = rn_sssr_bnd / rday                          !       -              -     
    123136!CDIR COLLAPSE 
    124137               DO jj = 1, jpj 
    125138                  DO ji = 1, jpi                             
    126                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
     139                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    127140                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
    128141                        &        / MAX(  sss_m(ji,jj), 1.e-20   ) 
     
    134147               END DO 
    135148               CALL iom_put( "erp", erp )                             ! freshwater flux damping 
     149               CALL wrk_dealloc( jpi,jpj,zsrp ) 
    136150            ENDIF 
    137151            ! 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r5397 r11384  
    2828   USE wrk_nemo        ! Memory Allocation 
    2929   USE timing          ! Timing 
     30   USE stopack 
    3031 
    3132   IMPLICIT NONE 
     
    4142 
    4243   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd0   ! geothermal heating trend 
     44   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend 
    4345   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read) 
    4446  
     
    8991      ! 
    9092      !                             !  Add the geothermal heat flux trend on temperature 
     93 
     94      IF( nn_spp_geot .GT. 0) THEN 
     95          qgh_trd1 = qgh_trd0 
     96          CALL spp_gen(kt, qgh_trd0, nn_spp_geot, rn_geot_sd, jk_spp_geot) 
     97      ENDIF 
    9198      DO jj = 2, jpjm1 
    9299         DO ji = 2, jpim1 
    93100            ik = mbkt(ji,jj) 
    94             zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik) 
     101            zqgh_trd = qgh_trd1(ji,jj) / fse3t(ji,jj,ik) 
    95102            tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd 
    96103         END DO 
     
    163170         ! 
    164171         ALLOCATE( qgh_trd0(jpi,jpj) )    ! allocation 
     172         ALLOCATE( qgh_trd1(jpi,jpj) )    ! allocation 
    165173         ! 
    166174         SELECT CASE ( nn_geoflx )        ! geothermal heat flux / (rauO * Cp) 
     
    192200            ! 
    193201         END SELECT 
     202         qgh_trd1 = qgh_trd0 
    194203         ! 
    195204      ELSE 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r4990 r11384  
    3939   USE timing         ! Timing 
    4040   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     41   USE stopack 
    4142 
    4243   IMPLICIT NONE 
     
    6768   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM) 
    6869   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points 
     70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_1, ahv_bbl_1   ! diffusive bbl flux coefficients at u and v-points 
    6971   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) 
    7072 
     
    8688         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     & 
    8789         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          & 
     90         &      ahu_bbl_1(jpi,jpj) , ahv_bbl_1(jpi,jpj) ,                                          & 
    8891         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc ) 
    8992         ! 
     
    194197      CALL wrk_alloc( jpi, jpj, zptb ) 
    195198      ! 
     199      ahu_bbl_1 = ahu_bbl 
     200      IF( nn_spp_ahubbl .GT. 0 ) THEN 
     201          CALL spp_gen(1 , ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl ) 
     202      ENDIF 
     203      ahv_bbl_1 = ahv_bbl 
     204      IF( nn_spp_ahvbbl .GT. 0 ) THEN 
     205          CALL spp_gen(1 , ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl ) 
     206      ENDIF 
     207      ! 
    196208      DO jn = 1, kjpt                                     ! tracer loop 
    197209         !                                                ! =========== 
     
    208220               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik) 
    209221               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
    210                   &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   & 
    211                   &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   & 
    212                   &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   & 
    213                   &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr 
     222                  &               + (   ahu_bbl_1(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   & 
     223                  &                   - ahu_bbl_1(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   & 
     224                  &                   + ahv_bbl_1(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   & 
     225                  &                   - ahv_bbl_1(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr 
    214226            END DO 
    215227         END DO 
     
    548560      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    549561 
    550                                         !* sign of grad(H) at u- and v-points 
    551       mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0 
     562                                        !* sign of grad(H) at u- and v-points; zero if grad(H) = 0 
     563      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0 
    552564      DO jj = 1, jpjm1 
    553565         DO ji = 1, jpim1 
    554             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
    555             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     566            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 
     567               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     568            ENDIF 
     569            ! 
     570            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 
     571               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     572            ENDIF 
    556573         END DO 
    557574      END DO 
     
    568585      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1) 
    569586      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1) 
    570  
    571587 
    572588      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r5120 r11384  
    3232   USE wrk_nemo        ! Memory allocation 
    3333   USE timing          ! Timing 
     34   USE stopack 
    3435 
    3536   IMPLICIT NONE 
     
    4344   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile 
    4445   !                                                               !  (key_traldf_ano only) 
     46#if defined key_traldf_c3d 
     47   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: ahtu0,ahtv0,ahtw0,ahtt0 
     48#endif 
     49#if defined key_traldf_c2d 
     50   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:  ) :: ahtu0,ahtv0,ahtw0,ahtt0 
     51#endif 
    4552 
    4653   !! * Substitutions 
     
    6875      ! 
    6976      rldf = 1     ! For active tracers the  
     77      r_fact_lap(:,:,:) = 1.0 
    7078 
    7179      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     
    7482         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    7583      ENDIF 
     84 
     85#if defined key_traldf_c3d 
     86         IF( kt .eq. nit000 .and. (nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt) .gt. 0 ) THEN 
     87             ALLOCATE ( ahtu0(jpi,jpj,jpk), ahtv0(jpi,jpj,jpk) ) 
     88             ALLOCATE ( ahtt0(jpi,jpj,jpk), ahtw0(jpi,jpj,jpk) ) 
     89             ahtu0 = ahtu 
     90             ahtv0 = ahtv 
     91             ahtw0 = ahtw 
     92             ahtt0 = ahtt 
     93         ENDIF 
     94#endif 
     95#if defined key_traldf_c2d 
     96         IF( kt .eq. nit000 .and. (nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt) .gt. 0 ) THEN 
     97             ALLOCATE ( ahtu0(jpi,jpj), ahtv0(jpi,jpj) ) 
     98             ALLOCATE ( ahtt0(jpi,jpj), ahtw0(jpi,jpj) ) 
     99             ahtu0 = ahtu 
     100             ahtv0 = ahtv 
     101             ahtw0 = ahtw 
     102             ahtt0 = ahtt 
     103         ENDIF 
     104#endif 
     105#if defined key_traldf_c3d || defined key_traldf_c2d 
     106         IF( nn_spp_ahtu .GT. 0) THEN 
     107             ahtu = ahtu0 
     108             CALL spp_aht(kt,ahtu,nn_spp_ahtu,rn_ahtu_sd,jk_spp_ahtu) 
     109         ENDIF 
     110         IF( nn_spp_ahtv .GT. 0) THEN 
     111             ahtv = ahtv0 
     112             CALL spp_aht(kt,ahtv,nn_spp_ahtv,rn_ahtv_sd,jk_spp_ahtv) 
     113         ENDIF 
     114         IF( nn_spp_ahtw .GT. 0) THEN 
     115             ahtw = ahtw0 
     116             CALL spp_aht(kt,ahtw,nn_spp_ahtw,rn_ahtw_sd,jk_spp_ahtw) 
     117         ENDIF 
     118         IF( nn_spp_ahtt .GT. 0) THEN 
     119             ahtt = ahtt0 
     120             CALL spp_aht(kt,ahtt,nn_spp_ahtt,rn_ahtt_sd,jk_spp_ahtt) 
     121         ENDIF 
     122#endif 
    76123 
    77124      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 
     
    214261      IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' ) 
    215262      IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' ) 
     263      IF( ln_traldf_hor .AND. ln_traldf_grif )    & 
     264            &   CALL ctl_stop( ' horizontal operator and Griffies triads not available; sitch to isoneutral operator' ) 
     265      IF( ln_traldf_grif .AND. ln_isfcav         )   & 
     266           CALL ctl_stop( ' ice shelf and traldf_grif not tested') 
    216267      IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso )   & 
    217268           CALL ctl_stop( '          eddy induced velocity on tracers',   & 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r5407 r11384  
    3232   USE wrk_nemo       ! Memory Allocation 
    3333   USE timing         ! Timing 
     34   USE stopack 
    3435 
    3536   IMPLICIT NONE 
     
    5152  
    5253   ! Module variables 
    53    REAL(wp) ::   xsi0r                           !: inverse of rn_si0 
     54   REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0 
    5455   REAL(wp) ::   xsi1r                           !: inverse of rn_si1 
    5556   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
     
    179180         !                                        ! ============================================== ! 
    180181         ! 
    181          !                                                ! ------------------------- ! 
     182         !  
     183         IF( nn_spp_qsi0 > 0 ) THEN 
     184             xsi0r = rn_si0 
     185             CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 
     186             xsi0r = 1.e0 / xsi0r 
     187         ENDIF 
     188         !                                               ! ------------------------- ! 
    182189         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration ! 
    183190            !                                             ! ------------------------- ! 
     
    221228!CDIR NOVERRCHK    
    222229                     DO ji = 1, jpi 
    223                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r    ) 
     230                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) ) 
    224231                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 
    225232                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 
     
    237244                  DO jj = 1, jpj 
    238245                     DO ji = 1, jpi 
    239                         zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r    ) 
     246                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r(ji,jj) ) 
    240247                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 
    241248                        zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
     
    267274            !                                             ! ------------------------- ! 
    268275            ! 
    269             IF( lk_vvl ) THEN                                  !* variable volume 
     276            IF( lk_vvl .OR. nn_spp_qsi0 > 0 ) THEN        !* variable volume 
     277 
    270278               zz0   =        rn_abs   * r1_rau0_rcp 
    271279               zz1   = ( 1. - rn_abs ) * r1_rau0_rcp 
     
    273281                  DO jj = 1, jpj 
    274282                     DO ji = 1, jpi 
    275                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    276                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
     283                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
     284                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    277285                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )  
    278286                     END DO 
     
    283291                  DO jj = 1, jpj 
    284292                     DO ji = 1, jpi 
    285                         zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 
    286                         zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 
     293                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 
     294                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 
    287295                        fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 
    288296                     END DO 
     
    445453         !                       ! ===================================== ! 
    446454         ! 
     455         ALLOCATE( xsi0r(jpi,jpj) ) 
    447456         xsi0r = 1.e0 / rn_si0 
    448457         xsi1r = 1.e0 / rn_si1 
     
    499508!CDIR NOVERRCHK    
    500509                        DO ji = 1, jpi 
    501                            zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r    ) 
     510                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) ) 
    502511                           zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 
    503512                           zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 
     
    540549                  DO jj = 1, jpj                              ! top 400 meters 
    541550                     DO ji = 1, jpi 
    542                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    543                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
     551                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
     552                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    544553                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)  
    545554                     END DO 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90

    r5215 r11384  
    2929   USE lib_mpp        ! MPP library 
    3030   USE wrk_nemo       ! Memory allocation 
     31   USE stopack      
    3132 
    3233   IMPLICIT NONE 
     
    6869      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    6970      IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt ) 
     71      IF( ln_dyn_trd .AND. ln_sppt_dyn ) CALL dyn_sppt_collect( putrd, pvtrd, ktrd, kt ) 
    7072          
    7173      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r4990 r11384  
    2929   USE iom            ! I/O manager library 
    3030   USE lib_mpp        ! MPP library 
     31   USE stopack 
    3132   USE wrk_nemo       ! Memory allocation 
    3233 
     
    3839   REAL(wp) ::   r2dt   ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0 
    3940 
    40    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt   ! use to store the temperature trends 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt  ! use to store the temperature trends 
     42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_evd  ! store avt_evd to calculate EVD trend 
    4143 
    4244   !! * Substitutions 
     
    5557      !!                  ***  FUNCTION trd_tra_alloc  *** 
    5658      !!--------------------------------------------------------------------- 
    57       ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc ) 
     59      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc ) 
    5860      ! 
    5961      IF( lk_mpp             )   CALL mpp_sum ( trd_tra_alloc ) 
     
    104106                                 ztrds(:,:,:) = 0._wp 
    105107                                 CALL trd_tra_mng( trdt, ztrds, ktrd, kt ) 
     108         CASE( jptra_evd )   ;   avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 
    106109         CASE DEFAULT                 ! other trends: masked trends 
    107110            trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)              ! mask & store 
     
    128131            zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp 
    129132            DO jk = 2, jpk 
    130                zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     133               zwt(:,:,jk) = avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
    131134               zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
    132135            END DO 
     
    138141            END DO 
    139142            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )   
     143            ! 
     144            !                         ! Also calculate EVD trend at this point.  
     145            zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes 
     146            DO jk = 2, jpk 
     147               zwt(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     148               zws(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     149            END DO 
     150            ! 
     151            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp 
     152            DO jk = 1, jpkm1 
     153               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk) 
     154               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk)  
     155            END DO 
     156            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt )   
    140157            ! 
    141158            CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt ) 
     
    233250      !                   ! 3D output of tracers trends using IOM interface 
    234251      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt ) 
    235  
    236       !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     252      IF( ln_tra_trd .AND. ln_sppt_tra )  CALL tra_sppt_collect( ptrdx, ptrdy, ktrd, kt ) 
     253 
     254      !  Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    237255      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt ) 
    238256 
     
    285303      !! ** Purpose :   output 3D tracer trends using IOM 
    286304      !!---------------------------------------------------------------------- 
    287       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend  
    288       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend 
    289       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index 
    290       INTEGER                   , INTENT(in   ) ::   kt      ! time step 
    291       !! 
    292       INTEGER ::   ji, jj, jk   ! dummy loop indices 
    293       INTEGER ::   ikbu, ikbv   ! local integers 
    294       REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace  
    295       !!---------------------------------------------------------------------- 
    296       ! 
    297 !!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added 
    298       ! 
    299       SELECT CASE( ktrd ) 
    300       CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection 
    301                                CALL iom_put( "strd_xad" , ptrdy ) 
    302       CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection 
    303                                CALL iom_put( "strd_yad" , ptrdy ) 
    304       CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection 
    305                                CALL iom_put( "strd_zad" , ptrdy ) 
    306                                IF( .NOT. lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface 
    307                                   CALL wrk_alloc( jpi, jpj, z2dx, z2dy ) 
    308                                   z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1) 
    309                                   z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1) 
    310                                   CALL iom_put( "ttrd_sad", z2dx ) 
    311                                   CALL iom_put( "strd_sad", z2dy ) 
    312                                   CALL wrk_dealloc( jpi, jpj, z2dx, z2dy ) 
    313                                ENDIF 
    314       CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion 
    315                                CALL iom_put( "strd_ldf" , ptrdy ) 
    316       CASE( jptra_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution) 
    317                                CALL iom_put( "strd_zdf" , ptrdy ) 
    318       CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution) 
    319                                CALL iom_put( "strd_zdfp", ptrdy ) 
    320       CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping) 
    321                                CALL iom_put( "strd_dmp" , ptrdy ) 
    322       CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer 
    323                                CALL iom_put( "strd_bbl" , ptrdy ) 
    324       CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing 
    325                                CALL iom_put( "strd_npc" , ptrdy ) 
    326       CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx )        ! surface forcing + runoff (ln_rnf=T) 
    327                                CALL iom_put( "strd_cdt" , ptrdy ) 
    328       CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature) 
    329       CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature) 
    330       CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter 
    331                                CALL iom_put( "strd_atf" , ptrdy ) 
    332       END SELECT 
    333       ! 
     305     REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend  
     306     REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend 
     307     INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index 
     308     INTEGER                   , INTENT(in   ) ::   kt      ! time step 
     309     !! 
     310     INTEGER ::   ji, jj, jk   ! dummy loop indices 
     311     INTEGER ::   ikbu, ikbv   ! local integers 
     312     REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace  
     313     !!---------------------------------------------------------------------- 
     314     ! 
     315     !!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added 
     316     ! 
     317     ! Trends evaluated every time step that could go to the standard T file and can be output every ts into a 1ts file if 1ts output is selected 
     318     SELECT CASE( ktrd ) 
     319     ! This total trend is done every time step 
     320     CASE( jptra_tot  )   ;   CALL iom_put( "ttrd_tot" , ptrdx )           ! model total trend 
     321        CALL iom_put( "strd_tot" , ptrdy ) 
     322     END SELECT 
     323 
     324     ! These trends are done every second time step. When 1ts output is selected must go different (2ts) file from standard T-file 
     325     IF( MOD( kt, 2 ) == 0 ) THEN 
     326        SELECT CASE( ktrd ) 
     327        CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection 
     328           CALL iom_put( "strd_xad" , ptrdy ) 
     329        CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection 
     330           CALL iom_put( "strd_yad" , ptrdy ) 
     331        CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection 
     332           CALL iom_put( "strd_zad" , ptrdy ) 
     333           IF( .NOT. lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface 
     334              CALL wrk_alloc( jpi, jpj, z2dx, z2dy ) 
     335              z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1) 
     336              z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1) 
     337              CALL iom_put( "ttrd_sad", z2dx ) 
     338              CALL iom_put( "strd_sad", z2dy ) 
     339              CALL wrk_dealloc( jpi, jpj, z2dx, z2dy ) 
     340           ENDIF 
     341        CASE( jptra_totad  ) ;   CALL iom_put( "ttrd_totad" , ptrdx )      ! total   advection 
     342           CALL iom_put( "strd_totad" , ptrdy ) 
     343        CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion 
     344           CALL iom_put( "strd_ldf" , ptrdy ) 
     345        CASE( jptra_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution) 
     346           CALL iom_put( "strd_zdf" , ptrdy ) 
     347        CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution) 
     348           CALL iom_put( "strd_zdfp", ptrdy ) 
     349        CASE( jptra_evd )    ;   CALL iom_put( "ttrd_evd", ptrdx )         ! EVD trend (convection) 
     350           CALL iom_put( "strd_evd", ptrdy ) 
     351        CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping) 
     352           CALL iom_put( "strd_dmp" , ptrdy ) 
     353        CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer 
     354           CALL iom_put( "strd_bbl" , ptrdy ) 
     355        CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing 
     356           CALL iom_put( "strd_npc" , ptrdy ) 
     357        CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature) 
     358        CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx(:,:,1) ) ! surface forcing + runoff (ln_rnf=T) 
     359           CALL iom_put( "strd_cdt" , ptrdy(:,:,1) )        ! output as 2D surface fields 
     360        CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature) 
     361        END SELECT 
     362        ! the Asselin filter trend  is also every other time step but needs to be lagged one time step 
     363        ! Even when 1ts output is selected can go to the same (2ts) file as the trends plotted every even time step. 
     364     ELSE IF( MOD( kt, 2 ) == 1 ) THEN 
     365        SELECT CASE( ktrd ) 
     366        CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter 
     367           CALL iom_put( "strd_atf" , ptrdy ) 
     368        END SELECT 
     369     END IF 
     370     ! 
    334371   END SUBROUTINE trd_tra_iom 
    335372 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r5332 r11384  
    2626   USE wrk_nemo        ! Memory Allocation 
    2727   USE phycst, ONLY: vkarmn 
     28   USE stopack 
    2829 
    2930   IMPLICIT NONE 
     
    5253   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
    5354   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction 
    54    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM) 
     55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d,bfrcoef2d0   ! 2D bottom/top drag coefficient (PUBLIC for TAM) 
    5556 
    5657   !! * Substitutions 
     
    6869      !!                ***  FUNCTION zdf_bfr_alloc  *** 
    6970      !!---------------------------------------------------------------------- 
    70       ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 
     71      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), bfrcoef2d0(jpi,jpj),STAT=zdf_bfr_alloc ) 
    7172      ! 
    7273      IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc ) 
     
    105106         WRITE(numout,*) 'zdf_bfr : Set bottom friction coefficient (non-linear case)' 
    106107         WRITE(numout,*) '~~~~~~~~' 
     108      ENDIF 
     109      !  
     110      IF( nn_spp_bfr > 0 ) THEN 
     111         bfrcoef2d = bfrcoef2d0 
     112         CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr) 
    107113      ENDIF 
    108114      ! 
     
    486492      ENDIF 
    487493      ! 
     494      bfrcoef2d0(:,:) = bfrcoef2d(:,:) 
    488495      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init') 
    489496      ! 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r4990 r11384  
    1919   USE zdf_oce         ! ocean vertical physics variables 
    2020   USE zdfkpp          ! KPP vertical mixing 
     21   USE trd_oce         ! trends: ocean variables 
     22   USE trdtra          ! trends manager: tracers  
    2123   USE in_out_manager  ! I/O manager 
    2224   USE iom             ! for iom_put 
    2325   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2426   USE timing          ! Timing 
     27   USE stopack          
    2528 
    2629   IMPLICIT NONE 
     
    2831 
    2932   PUBLIC   zdf_evd    ! called by step.F90 
     33   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rn_avevd0 
    3034 
    3135   !! * Substitutions 
     
    6771         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    6872         IF(lwp) WRITE(numout,*) 
     73         ALLOCATE ( rn_avevd0(jpi,jpj) ) 
     74         rn_avevd0(:,:) = rn_avevd 
    6975      ENDIF 
    7076 
    7177      zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application 
     78 
     79      IF(nn_spp_aevd.GT.0) THEN 
     80         rn_avevd0(:,:) = rn_avevd 
     81         CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd) 
     82      ENDIF 
    7283 
    7384      SELECT CASE ( nn_evdm ) 
     
    8697                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 
    8798#endif 
    88                      avt (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk) 
    89                      avm (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk) 
    90                      avmu(ji  ,jj  ,jk) = rn_avevd * umask(ji  ,jj  ,jk) 
    91                      avmu(ji-1,jj  ,jk) = rn_avevd * umask(ji-1,jj  ,jk) 
    92                      avmv(ji  ,jj  ,jk) = rn_avevd * vmask(ji  ,jj  ,jk) 
    93                      avmv(ji  ,jj-1,jk) = rn_avevd * vmask(ji  ,jj-1,jk) 
     99                     avt (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk) 
     100                     avm (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk) 
     101                     avmu(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji  ,jj  ,jk) 
     102                     avmu(ji-1,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji-1,jj  ,jk) 
     103                     avmv(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj  ,jk) 
     104                     avmv(ji  ,jj-1,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj-1,jk) 
    94105                  ENDIF 
    95106               END DO 
     
    113124                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   & 
    114125#endif 
    115                      avt(ji,jj,jk) = rn_avevd * tmask(ji,jj,jk) 
     126                     avt(ji,jj,jk) = rn_avevd0(ji,jj) * tmask(ji,jj,jk) 
    116127               END DO 
    117128            END DO 
     
    122133      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
    123134      CALL iom_put( "avt_evd", zavt_evd )              ! output this change 
     135      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 
    124136      ! 
    125137      IF( nn_timing == 1 )  CALL timing_stop('zdf_evd') 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r5109 r11384  
    3232   USE timing         ! Timing 
    3333   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     34   USE stopack 
    3435 
    3536   IMPLICIT NONE 
     
    4243   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfgls = .TRUE.   !: TKE vertical mixing flag 
    4344   ! 
    44    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en      !: now turbulent kinetic energy 
    4545   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length 
    4646   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function 
    47    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k   ! not enhanced Kz 
    48    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k   ! not enhanced Kz 
    49    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k  ! not enhanced Kz 
    50    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmv_k  ! not enhanced Kz 
    5147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
    5248   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points 
     
    120116      !!                ***  FUNCTION zdf_gls_alloc  *** 
    121117      !!---------------------------------------------------------------------- 
    122       ALLOCATE( en(jpi,jpj,jpk),  mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    123          &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                    & 
    124          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk),                    & 
    125          &      ustars2(jpi,jpj), ustarb2(jpi,jpj)                      , STAT= zdf_gls_alloc ) 
     118      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
     119         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
    126120         ! 
    127121      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    329323      !  
    330324      ! One level below 
    331       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
     325      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
     326          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    332327      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    333328      z_elem_a(:,:,2) = 0._wp  
     
    350345      z_elem_a(:,:,2) = 0._wp 
    351346      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
    352       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     347      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
     348           &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
    353349 
    354350      en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
     
    811807            END DO 
    812808         END DO 
     809         IF(nn_spp_avt > 0 ) CALL spp_gen(kt,avt(:,:,jk),nn_spp_avt,rn_avt_sd,jk) 
     810         IF(nn_spp_avm > 0 ) CALL spp_gen(kt,avm(:,:,jk),nn_spp_avm,rn_avm_sd,jk) 
    813811      END DO 
    814812      ! 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r5407 r11384  
    5353   USE timing         ! Timing 
    5454   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     55   USE stopack 
    5556 
    5657   IMPLICIT NONE 
     
    8586   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    8687 
    87    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2] 
     88   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 
     89 
    8890   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    8991   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    90    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz 
     92 
    9293#if defined key_c1d 
    9394   !                                                                        !!** 1D cfg only  **   ('key_c1d') 
     
    115116         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          & 
    116117#endif 
    117          &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     &  
    118          &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          & 
    119          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      ) 
     118         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc      ) 
    120119         ! 
    121120      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc ) 
     
    178177         avmu(:,:,:) = avmu_k(:,:,:)  
    179178         avmv(:,:,:) = avmv_k(:,:,:)  
    180       ENDIF  
     179      ENDIF 
     180      ! 
     181      IF( nn_spp_tkelc > 0 ) THEN 
     182          rn_lc0 = rn_lc 
     183          CALL spp_gen(kt,rn_lc0,nn_spp_tkelc,rn_tkelc_sd,    jk_spp_tkelc ) 
     184      ENDIF 
     185      IF( nn_spp_tkedf > 0 ) THEN 
     186          rn_ediff0 = rn_ediff 
     187          CALL spp_gen(kt,rn_ediff0,nn_spp_tkedf,rn_tkedf_sd, jk_spp_tkedf ) 
     188      ENDIF 
     189      IF( nn_spp_tkeds > 0 ) THEN 
     190          rn_ediss0 = rn_ediss 
     191          CALL spp_gen(kt,rn_ediss0,nn_spp_tkeds,rn_tkeds_sd, jk_spp_tkeds ) 
     192      ENDIF 
     193      IF( nn_spp_tkebb > 0 ) THEN 
     194          rn_ebb0 = rn_ebb 
     195          CALL spp_gen(kt,rn_ebb0,nn_spp_tkebb,rn_tkebb_sd,   jk_spp_tkebb ) 
     196      ENDIF 
     197      IF( nn_spp_tkefr > 0 ) THEN 
     198          rn_efr0 = rn_efr 
     199          CALL spp_gen(kt,rn_efr0,nn_spp_tkefr,rn_tkefr_sd,   jk_spp_tkefr ) 
     200      ENDIF 
    181201      ! 
    182202      CALL tke_tke      ! now tke (en) 
     
    188208      avmu_k(:,:,:) = avmu(:,:,:)  
    189209      avmv_k(:,:,:) = avmv(:,:,:)  
     210      ! 
     211      IF (  kt .eq. nitend ) THEN 
     212        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 
     213      ENDIF 
    190214      ! 
    191215   END SUBROUTINE zdf_tke 
     
    214238      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3 
    215239      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient 
    216       REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars 
    217       REAL(wp) ::   zfact1, zfact2, zfact3          !    -         - 
     240      REAL(wp) ::   zesh2                           ! temporary scalars 
     241      REAL(wp) ::   zfact1                          !    -         - 
    218242      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         - 
    219243      REAL(wp) ::   ztau  , zdif                    !    -         - 
     
    222246!!bfr      REAL(wp) ::   zebot                           !    -         - 
    223247      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc 
    224       REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc 
     248      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc, zbbrau,zfact2,zfact3 
    225249      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 
    226250      !!-------------------------------------------------------------------- 
     
    229253      ! 
    230254      CALL wrk_alloc( jpi,jpj, imlc )    ! integer 
    231       CALL wrk_alloc( jpi,jpj, zhlc )  
     255      CALL wrk_alloc( jpi,jpj, zhlc ) 
     256      CALL wrk_alloc( jpi,jpj, zbbrau ) 
     257      CALL wrk_alloc( jpi,jpj, zfact2 ) 
     258      CALL wrk_alloc( jpi,jpj, zfact3 ) 
    232259      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    233260      ! 
    234       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     261      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation 
    235262      zfact1 = -.5_wp * rdt  
    236       zfact2 = 1.5_wp * rdt * rn_ediss 
    237       zfact3 = 0.5_wp       * rn_ediss 
     263      zfact2 = 1.5_wp * rdt * rn_ediss0 
     264      zfact3 = 0.5_wp       * rn_ediss0 
    238265      ! 
    239266      ! 
     
    250277      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    251278         DO ji = fs_2, fs_jpim1   ! vector opt. 
    252             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     279            en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1) 
    253280         END DO 
    254281      END DO 
     
    315342                  !                                           ! vertical velocity due to LC 
    316343                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 
    317                   zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
     344                  zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    318345                  !                                           ! TKE Langmuir circulation source term 
    319                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     346                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   & 
     347                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     348 
    320349               END DO 
    321350            END DO 
     
    360389               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    361390               zd_lw(ji,jj,jk) = zzd_lw 
    362                zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
     391               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
    363392               ! 
    364393               !                                   ! right hand side in en 
    365394               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    366                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     395                  &                                 + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    367396                  &                                 * wmask(ji,jj,jk) 
    368397            END DO 
     
    420449            DO jj = 2, jpjm1 
    421450               DO ji = fs_2, fs_jpim1   ! vector opt. 
    422                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    423                      &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     451                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     452                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    424453               END DO 
    425454            END DO 
     
    429458            DO ji = fs_2, fs_jpim1   ! vector opt. 
    430459               jk = nmln(ji,jj) 
    431                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    432                   &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     460               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     461                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    433462            END DO 
    434463         END DO 
     
    445474                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    446475                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    447                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    448                      &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     476                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     477                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    449478               END DO 
    450479            END DO 
     
    455484      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
    456485      CALL wrk_dealloc( jpi,jpj, zhlc )  
     486      CALL wrk_dealloc( jpi,jpj, zbbrau )  
     487      CALL wrk_dealloc( jpi,jpj, zfact2 )  
     488      CALL wrk_dealloc( jpi,jpj, zfact3 )  
    457489      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    458490      ! 
     
    637669            DO ji = fs_2, fs_jpim1   ! vector opt. 
    638670               zsqen = SQRT( en(ji,jj,jk) ) 
    639                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     671               zav   = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen 
    640672               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
    641673               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     
    643675            END DO 
    644676         END DO 
     677         IF(nn_spp_avt > 0 ) CALL spp_gen(1 ,avt(:,:,jk),nn_spp_avt,rn_avt_sd, jk_spp_avt, jk) 
     678         IF(nn_spp_avm > 0 ) CALL spp_gen(1 ,avm(:,:,jk),nn_spp_avm,rn_avm_sd, jk_spp_avm, jk) 
    645679      END DO 
    646680      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     
    710744      !!---------------------------------------------------------------------- 
    711745      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    712       INTEGER ::   ios 
     746      INTEGER ::   ios, ierr 
    713747      !! 
    714748      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   & 
     
    767801         rn_mxl0 = rmxl_min 
    768802      ENDIF 
     803 
     804      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc 
     805      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff 
     806      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss 
     807      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb 
     808      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr 
    769809       
    770       IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     810      IF( nn_etau == 2 ) THEN 
     811          ierr = zdf_mxl_alloc() 
     812          nmln(:,:) = nlb10           ! Initialization of nmln 
     813      ENDIF 
    771814 
    772815      !                               !* depth of penetration of surface tke 
     
    836879              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
    837880              ! 
    838               avt_k (:,:,:) = avt (:,:,:) 
    839               avm_k (:,:,:) = avm (:,:,:) 
    840               avmu_k(:,:,:) = avmu(:,:,:) 
    841               avmv_k(:,:,:) = avmv(:,:,:) 
    842               ! 
    843881              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
    844882           ENDIF 
    845883        ELSE                                   !* Start from rest 
    846884           en(:,:,:) = rn_emin * tmask(:,:,:) 
    847            DO jk = 1, jpk                           ! set the Kz to the background value 
    848               avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    849               avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
    850               avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
    851               avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    852            END DO 
    853885        ENDIF 
    854886        ! 
     887        avt_k (:,:,:) = avt (:,:,:) 
     888        avm_k (:,:,:) = avm (:,:,:) 
     889        avmu_k(:,:,:) = avmu(:,:,:) 
     890        avmv_k(:,:,:) = avmv(:,:,:) 
     891         
    855892     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    856893        !                                   ! ------------------- 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5407 r11384  
    8383   USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges  
    8484   USE sbc_oce, ONLY: lk_oasis 
     85   USE stopack 
    8586   USE stopar 
    8687   USE stopts 
     
    453454                            CALL dia_hsb_init   ! heat content, salt content and volume budgets 
    454455                            CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends 
     456                            CALL stopack_init   ! STOPACK scheme 
    455457      IF( lk_diaobs     ) THEN                  ! Observation & model comparison 
    456458                            CALL dia_obs_init            ! Initialize observational data 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5510 r11384  
    105105                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    106106      ENDIF 
     107      IF( ln_stopack )   CALL stopack_pert( kstp ) 
    107108                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
    108109                                                      ! clem: moved here for bdy ice purpose 
     
    110111      ! Update stochastic parameters and random T/S fluctuations 
    111112      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    112                         CALL sto_par( kstp )          ! Stochastic parameters 
     113      IF( ln_sto_eos )   CALL sto_par( kstp )          ! Stochastic parameters 
     114      IF( ln_sto_eos )   CALL sto_pts( tsn  )          ! Random T/S fluctuations 
     115      IF( ln_skeb  )     CALL skeb_comp( kstp ) 
    113116 
    114117      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    134137      ENDIF 
    135138      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths 
    136          DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO 
     139         IF ( nn_spp_arnf .GT. 0 ) THEN 
     140              rn_avt_rnf0 = rn_avt_rnf 
     141              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf ) 
     142         ENDIF 
     143         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO 
    137144      ENDIF 
    138145      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity 
     
    148155      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' ) 
    149156      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' ) 
     157      IF( lrst_oce .AND. ln_stopack)   CALL stopack_rst( kstp, 'WRITE' ) 
    150158      ! 
    151159      !  LATERAL  PHYSICS 
    152160      ! 
    153161      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing 
    154          IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations 
    155162                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density 
    156163         IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
     
    188195          ! Note that the computation of vertical velocity above, hence "after" sea level 
    189196          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    190             IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations 
    191197                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    192198            IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
     
    201207                                  va(:,:,:) = 0.e0 
    202208          IF(  ln_asmiau .AND. & 
    203              & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment 
    204           IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified) 
    205           IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends 
    206                                   CALL dyn_adv      ( kstp )   ! advection (vector or flux form) 
    207                                   CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis 
    208                                   CALL dyn_ldf      ( kstp )   ! lateral mixing 
    209           IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified) 
     209             & ln_dyninc       )  CALL dyn_asm_inc   ( kstp )   ! apply dynamics assimilation increment 
     210          IF( ln_sppt_dyn )       CALL dyn_sppt_apply( kstp, 0 ) 
     211          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! subtract Neptune velocities (simplified) 
     212          IF( lk_bdy           )  CALL bdy_dyn3d_dmp ( kstp )   ! bdy damping trends 
     213                                  CALL dyn_adv       ( kstp )   ! advection (vector or flux form) 
     214                                  CALL dyn_vor       ( kstp )   ! vorticity term including Coriolis 
     215                                  CALL dyn_ldf       ( kstp )   ! lateral mixing 
     216          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! add Neptune velocities (simplified) 
    210217#if defined key_agrif 
    211218          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge 
     
    231238      IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag 
    232239      IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis 
     240                            CALL dia_prod( kstp )        ! ocean model: product diagnostics 
    233241                            CALL dia_wri( kstp )         ! ocean model: outputs 
    234242      ! 
     
    248256                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero 
    249257 
    250       IF(  ln_asmiau .AND. & 
    251          & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment 
     258      IF( ln_sppt_tra )      CALL tra_sppt_apply ( kstp, 0 ) 
     259      IF( ln_asmiau .AND. & 
     260        & ln_trainc      )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment 
    252261                             CALL tra_sbc    ( kstp )       ! surface boundary condition 
    253262      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr 
     
    265274      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge 
    266275#endif 
     276      IF( ln_sppt_tra )      CALL tra_sppt_apply ( kstp, 1 ) 
    267277                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields 
    268278 
     
    270280         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    271281                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    272             IF( ln_sto_eos ) CALL sto_pts( tsn )                 ! Random T/S fluctuations 
     282         IF( ln_sppt_tra )   CALL tra_sppt_apply ( kstp, 2 ) 
    273283                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation 
    274284            IF( ln_zps .AND. .NOT. ln_isfcav)                                & 
     
    281291      ELSE                                                  ! centered hpg  (eos then time stepping) 
    282292         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    283             IF( ln_sto_eos ) CALL sto_pts( tsn )    ! Random T/S fluctuations 
    284293                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
    285294         IF( ln_zps .AND. .NOT. ln_isfcav)                                   & 
     
    293302         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    294303                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
     304         IF( ln_sppt_tra )   THEN 
     305                                CALL tra_sppt_apply ( kstp, 2 ) 
     306                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
     307                                IF( ln_zps .AND. .NOT. ln_isfcav)                & 
     308                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     309                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     310                                IF( ln_zps .AND.       ln_isfcav)                                   &  
     311                                & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     312                                &                   rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     313                                &                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    295314      ENDIF 
    296315 
     
    309328 
    310329                               CALL dyn_bfr( kstp )         ! bottom friction 
     330            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 1 ) 
    311331                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    312332      ELSE 
     
    316336        IF(  ln_asmiau .AND. & 
    317337           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
     338        IF( ln_sppt_dyn )      CALL dyn_sppt_apply ( kstp, 0 ) 
    318339        IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
    319340        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified) 
     
    328349                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
    329350                               CALL dyn_bfr( kstp )         ! bottom friction 
     351            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 1 ) 
    330352                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    331353                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    332354      ENDIF 
    333355                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
     356            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 2 ) 
     357            IF( ln_skeb )      CALL skeb_apply ( kstp ) 
    334358 
    335359                               CALL ssh_swp( kstp )         ! swap of sea surface height 
     
    353377      ENDIF 
    354378      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file 
     379      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters 
    355380 
    356381      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r5501 r11384  
    9595   USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine) 
    9696   USE diaharm 
     97   USE diaprod          ! ocean model: product diagnostics 
    9798   USE flo_oce          ! floats variables 
    9899   USE floats           ! floats computation               (flo_stp routine) 
     
    110111   USE timing           ! Timing 
    111112 
     113   USE stopack          ! Stochastic physics 
     114 
    112115#if defined key_agrif 
    113116   USE agrif_opa_sponge ! Momemtum and tracers sponges 
Note: See TracChangeset for help on using the changeset viewer.