New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7367 for branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Ignore:
Timestamp:
2016-11-29T11:52:31+01:00 (8 years ago)
Author:
deazer
Message:

Contains merged code for CO5 reference.

Location:
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r7363 r7367  
    88   !!            3.2  ! 2009-09  (A.C.Coward)  Correction to include barotropic contribution 
    99   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     10   !!            3.4  ! 2011-11  (H. Liu) implementation of semi-implicit bottom friction option 
     11   !!                 ! 2012-06  (H. Liu) implementation of Log Layer bottom friction option 
    1012   !!---------------------------------------------------------------------- 
    1113 
     
    3032   PUBLIC   zdf_bfr_init    ! called by opa.F90 
    3133 
     34   REAL(wp), PARAMETER :: karman = 0.41_wp ! von Karman constant   
    3235   !                                    !!* Namelist nambfr: bottom friction namelist * 
    33    INTEGER  ::   nn_bfr    = 0           ! = 0/1/2/3 type of bottom friction  
    34    REAL(wp) ::   rn_bfri1  = 4.0e-4_wp   ! bottom drag coefficient (linear case)  
    35    REAL(wp) ::   rn_bfri2  = 1.0e-3_wp   ! bottom drag coefficient (non linear case) 
    36    REAL(wp) ::   rn_bfeb2  = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2] 
    37    REAL(wp) ::   rn_bfrien = 30._wp      ! local factor to enhance coefficient bfri 
    38    LOGICAL  ::   ln_bfr2d  = .false.     ! logical switch for 2D enhancement 
    39    LOGICAL , PUBLIC                            ::  ln_bfrimp = .false.  ! logical switch for implicit bottom friction 
     36   INTEGER  ::   nn_bfr      = 0           ! = 0/1/2/3 type of bottom friction  
     37   REAL(wp) ::   rn_bfri1    = 4.0e-4_wp   ! bottom drag coefficient (linear case)  
     38   REAL(wp) ::   rn_bfri2    = 1.0e-3_wp   ! bottom drag coefficient (non linear case) 
     39   REAL(wp) ::   rn_bfeb2    = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2] 
     40   REAL(wp) ::   rn_bfrien   = 30._wp      ! local factor to enhance coefficient bfri 
     41   REAL(wp) ::   rn_bfrz0    = 0.003_wp    ! bottom roughness for loglayer bfr coeff 
     42   LOGICAL  ::   ln_bfr2d    = .false.     ! logical switch for 2D enhancement 
     43   LOGICAL  ::   ln_loglayer = .false.     ! switch for log layer bfr coeff. 
     44   LOGICAL , PUBLIC                            ::  ln_bfrimp   = .false.  ! switch for implicit bottom friction 
    4045   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  bfrcoef2d            ! 2D bottom drag coefficient 
    4146 
     
    8287      INTEGER  ::   ikbu, ikbv   ! local integers 
    8388      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars 
     89      REAL(wp) ::   ztmp         ! temporary scalars 
    8490      !!---------------------------------------------------------------------- 
    8591      ! 
     
    9298         ! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]} 
    9399         ! 
     100 
     101         IF(ln_loglayer) THEN       ! "log layer" bottom friction coefficient 
     102#  if defined key_vectopt_loop 
     103           DO jj = 1, 1 
     104             DO ji = 1, jpij   ! vector opt. (forced unrolling) 
     105#  else 
     106           DO jj = 1, jpj 
     107             DO ji = 1, jpi 
     108#  endif 
     109                ztmp = 0.5_wp * fse3t(ji,jj,mbkt(ji,jj)) 
     110                ztmp = max(ztmp, rn_bfrz0) 
     111                bfrcoef2d(ji,jj) = ( log( ztmp / rn_bfrz0 ) / karman ) ** (-2) 
     112             END DO 
     113           END DO 
     114         ENDIF 
     115 
    94116# if defined key_vectopt_loop 
    95117         DO jj = 1, 1 
     
    117139            END DO 
    118140         END DO 
     141 
     142 
    119143         ! 
    120144         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
     
    141165      USE iom   ! I/O module for ehanced bottom friction file 
    142166      !! 
    143       INTEGER ::   inum         ! logical unit for enhanced bottom friction file 
    144       INTEGER ::   ji, jj       ! dummy loop indexes 
    145       INTEGER ::   ikbu, ikbv   ! temporary integers 
    146       INTEGER ::   ictu, ictv   !    -          - 
    147       REAL(wp) ::  zminbfr, zmaxbfr   ! temporary scalars 
    148       REAL(wp) ::  zfru, zfrv         !    -         - 
    149       !! 
    150       NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien, ln_bfrimp 
     167      INTEGER   ::   inum         ! logical unit for enhanced bottom friction file 
     168      INTEGER   ::   ji, jj       ! dummy loop indexes 
     169      INTEGER   ::   ikbu, ikbv   ! temporary integers 
     170      INTEGER   ::   ictu, ictv   !    -          - 
     171      REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars 
     172      REAL(wp)  ::   zfru, zfrv         !    -         - 
     173      !! 
     174      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 
     175                    &  rn_bfrien, ln_bfrimp, ln_loglayer 
    151176      !!---------------------------------------------------------------------- 
    152177      ! 
     
    212237         ENDIF 
    213238         bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
     239 
    214240         ! 
    215241         IF(ln_bfr2d) THEN  
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r7363 r7367  
    227227      ENDIF 
    228228      ! 
    229       !                              ! allocate zdfddm arrays 
     229      !                               ! allocate zdfddm arrays 
    230230      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' ) 
     231      !                               ! initialization to masked Kz 
     232      avs(:,:,:) = rn_avt0 * tmask(:,:,:)  
    231233      ! 
    232234   END SUBROUTINE zdf_ddm_init 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r7363 r7367  
    1515   USE prtctl          ! Print control 
    1616   USE iom             ! I/O library 
     17   USE eosbn2          ! Equation of state  
     18   USE phycst, ONLY : rau0 ! reference density  
     19   USE lbclnk  
    1720   USE lib_mpp         ! MPP library 
    1821   USE wrk_nemo        ! work arrays 
     
    2427 
    2528   PUBLIC   zdf_mxl       ! called by step.F90 
    26  
     29    
     30   ! Namelist variables for  namzdf_karaml  
     31  
     32   LOGICAL   :: ln_kara            ! Logical switch for calculating kara mixed 
     33                                     ! layer 
     34   LOGICAL   :: ln_kara_write25h   ! Logical switch for 25 hour outputs 
     35   INTEGER   :: jpmld_type         ! mixed layer type              
     36   REAL(wp)  :: ppz_ref            ! depth of initial T_ref  
     37   REAL(wp)  :: ppdT_crit          ! Critical temp diff   
     38   REAL(wp)  :: ppiso_frac         ! Fraction of ppdT_crit used  
     39    
     40   !Used for 25h mean 
     41   LOGICAL, PRIVATE :: kara_25h_init = .TRUE.   !Logical used to initalise 25h  
     42                                                !outputs. Necissary, because we need to 
     43                                                !initalise the kara_25h on the zeroth 
     44                                                !timestep (i.e in the nemogcm_init call) 
     45   REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h 
     46    
    2747   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP) 
     48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_kara  !: mixed layer depth of Kara et al   [m] 
    2849   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    2950   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3051   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
    31  
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m] 
     53    
    3254   !! * Substitutions 
    3355#  include "domzgr_substitute.h90" 
     
    4567      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    4668      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    47          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     69         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), &  
     70                                            hmld_tref(jpi,jpj), STAT= zdf_mxl_alloc )  
    4871         ! 
    4972         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
     
    5982      !!                    
    6083      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth 
    61       !!              with density criteria. 
     84      !!      with a simple density criteria. Also calculates the mixed layer  
     85      !!      depth of Kara et al by calling zdf_mxl_kara.  
    6286      !! 
    6387      !! ** Method  :   The mixed layer depth is the shallowest W depth with  
     
    78102      REAL(wp) ::   zrho_c = 0.01_wp    ! density criterion for mixed layer depth 
    79103      REAL(wp) ::   zavt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     104      REAL(wp) ::   t_ref               ! Reference temperature   
     105      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth   
    80106      !!---------------------------------------------------------------------- 
    81107      ! 
     
    104130      END DO 
    105131      ! depth of the mixing and mixed layers 
     132       
     133      CALL zdf_mxl_kara( kt )  
     134       
    106135      DO jj = 1, jpj 
    107136         DO ji = 1, jpi 
     
    113142         END DO 
    114143      END DO 
     144#if defined key_iomput  
    115145      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode 
    116146         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth 
    117147         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth 
    118148      ENDIF 
    119        
     149#endif 
     150 
     151      !For the AMM model assimiation uses a temperature based mixed layer depth   
     152      !This is defined here   
     153      DO jj = 1, jpj   
     154         DO ji = 1, jpi   
     155           hmld_tref(ji,jj)=fsdept(ji,jj,1  )    
     156           IF(tmask(ji,jj,1) > 0.)THEN   
     157             t_ref=tsn(ji,jj,1,jp_tem)  
     158             DO jk=2,jpk   
     159               IF(tmask(ji,jj,jk)==0.)THEN   
     160                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )   
     161                  EXIT   
     162               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN   
     163                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )   
     164               ELSE   
     165                  EXIT   
     166               ENDIF   
     167             ENDDO   
     168           ENDIF   
     169         ENDDO   
     170      ENDDO 
     171             
    120172      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    121173      ! 
     
    125177      ! 
    126178   END SUBROUTINE zdf_mxl 
    127  
     179    
     180    
     181   SUBROUTINE zdf_mxl_kara( kt )  
     182      !!----------------------------------------------------------------------------------  
     183      !!                    ***  ROUTINE zdf_mxl_kara  ***  
     184      !                                                                         
     185      !   Calculate mixed layer depth according to the definition of           
     186      !   Kara et al, 2000, JGR, 105, 16803.   
     187      !  
     188      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     189      !   density has increased by an amount equivalent to a temperature difference of   
     190      !   0.8C at the surface.  
     191      !  
     192      !   For other values of mld_type the mixed layer is calculated as the depth at   
     193      !   which the temperature differs by 0.8C from the surface temperature.   
     194      !                                                                         
     195      !   Original version: David Acreman                                       
     196      !  
     197      !!----------------------------------------------------------------------------------- 
     198      
     199      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index  
     200  
     201      NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, & 
     202      &                       ppiso_frac, ln_kara_write25h  
     203  
     204      ! Local variables                                                         
     205      REAL, DIMENSION(jpi,jpk) :: ppzdep      ! depth for use in calculating d(rho)  
     206      REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d  !Local version of tsn  
     207      LOGICAL :: ll_found(jpi,jpj)              ! Is T_b to be found by interpolation ?  
     208      LOGICAL :: ll_belowml(jpi,jpj,jpk)        ! Flag points below mixed layer when ll_found=F  
     209      INTEGER :: ji, jj, jk                     ! loop counter  
     210      INTEGER :: ik_ref(jpi,jpj)                ! index of reference level  
     211      INTEGER :: ik_iso(jpi,jpj)                ! index of last uniform temp level  
     212      REAL    :: zT(jpi,jpj,jpk)                ! Temperature or denisty  
     213      REAL    :: zT_ref(jpi,jpj)                ! reference temperature  
     214      REAL    :: zT_b                           ! base temperature  
     215      REAL    :: zdTdz(jpi,jpj,jpk-2)           ! gradient of zT  
     216      REAL    :: zmoddT(jpi,jpj,jpk-2)          ! Absolute temperature difference  
     217      REAL    :: zdz                            ! depth difference  
     218      REAL    :: zdT                            ! temperature difference  
     219      REAL    :: zdelta_T(jpi,jpj)              ! difference critereon  
     220      REAL    :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities 
     221      INTEGER, SAVE :: idt, i_steps 
     222      INTEGER, SAVE :: i_cnt_25h     ! Counter for 25 hour means 
     223      
     224  
     225      !!-------------------------------------------------------------------------------------  
     226  
     227      IF( kt == nit000 ) THEN  
     228         ! Default values  
     229         ln_kara      = .FALSE. 
     230         ln_kara_write25h = .FALSE.  
     231         jpmld_type   = 1      
     232         ppz_ref      = 10.0  
     233         ppdT_crit    = 0.2   
     234         ppiso_frac   = 0.1    
     235         ! Read namelist  
     236         REWIND ( numnam )     
     237         READ   ( numnam, namzdf_karaml )  
     238         WRITE(numout,*) '===== Kara mixed layer ====='  
     239         WRITE(numout,*) 'ln_kara = ',    ln_kara 
     240         WRITE(numout,*) 'jpmld_type = ', jpmld_type  
     241         WRITE(numout,*) 'ppz_ref = ',    ppz_ref  
     242         WRITE(numout,*) 'ppdT_crit = ',  ppdT_crit  
     243         WRITE(numout,*) 'ppiso_frac = ', ppiso_frac 
     244         WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h 
     245         WRITE(numout,*) '============================'  
     246       
     247         IF ( .NOT.ln_kara ) THEN 
     248            WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated"  
     249         ELSE 
     250            IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) ) 
     251            IF ( ln_kara_write25h .AND. kara_25h_init ) THEN 
     252               i_cnt_25h = 0 
     253               IF (.NOT. ALLOCATED(hmld_kara_25h) ) & 
     254               &   ALLOCATE( hmld_kara_25h(jpi,jpj) ) 
     255               hmld_kara_25h = 0._wp 
     256               IF( nacc == 1 ) THEN 
     257                  idt = INT(rdtmin) 
     258               ELSE 
     259                  idt = INT(rdt) 
     260               ENDIF 
     261               IF( MOD( 3600,idt) == 0 ) THEN  
     262                  i_steps = 3600 / idt   
     263               ELSE  
     264                  CALL ctl_stop('STOP', & 
     265                  & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// & 
     266                  & ' = 0 otherwise no hourly values are possible')  
     267               ENDIF   
     268            ENDIF 
     269         ENDIF 
     270      ENDIF 
     271       
     272      IF ( ln_kara ) THEN 
     273          
     274         !set critical ln_kara 
     275         ztsn_2d = tsn(:,:,1,:) 
     276         ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit 
     277      
     278         ! Set the mixed layer depth criterion at each grid point  
     279         ppzdep = 0._wp 
     280         IF( jpmld_type == 1 ) THEN                                          
     281            CALL eos ( tsn(:,:,1,:), & 
     282            &          ppzdep(:,:), zRHO1(:,:) )  
     283            CALL eos ( ztsn_2d(:,:,:), & 
     284            &           ppzdep(:,:), zRHO2(:,:) )  
     285            zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     286            ! RHO from eos (2d version) doesn't calculate north or east halo:  
     287            CALL lbc_lnk( zdelta_T, 'T', 1. )  
     288            zT(:,:,:) = rhop(:,:,:)  
     289         ELSE  
     290            zdelta_T(:,:) = ppdT_crit                       
     291            zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     292         ENDIF  
     293      
     294         ! Calculate the gradient of zT and absolute difference for use later  
     295         DO jk = 1 ,jpk-2  
     296            zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     297            zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     298         END DO  
     299      
     300         ! Find density/temperature at the reference level (Kara et al use 10m).           
     301         ! ik_ref is the index of the box centre immediately above or at the reference level  
     302         ! Find ppz_ref in the array of model level depths and find the ref     
     303         ! density/temperature by linear interpolation.                                    
     304         ik_ref = -1 
     305         DO jk = jpkm1, 2, -1  
     306            WHERE( fsdept(:,:,jk) > ppz_ref )  
     307               ik_ref(:,:) = jk - 1  
     308               zT_ref(:,:) = zT(:,:,jk-1) + & 
     309               &             zdTdz(:,:,jk-1) * ( ppz_ref - fsdept(:,:,jk-1) )  
     310            ENDWHERE  
     311         END DO 
     312         IF ( ANY( ik_ref  < 0 ) .OR. ANY( ik_ref  > jpkm1 ) ) THEN 
     313            CALL ctl_stop( "STOP", & 
     314            & "zdf_mxl_kara: unable to find reference level for kara ML" )  
     315         ELSE 
     316            ! If the first grid box centre is below the reference level then use the  
     317            ! top model level to get zT_ref  
     318            WHERE( fsdept(:,:,1) > ppz_ref )   
     319               zT_ref = zT(:,:,1)  
     320               ik_ref = 1  
     321            ENDWHERE  
     322      
     323            ! Search for a uniform density/temperature region where adjacent levels           
     324            ! differ by less than ppiso_frac * deltaT.                                       
     325            ! ik_iso is the index of the last level in the uniform layer   
     326            ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     327            ik_iso(:,:)   = ik_ref(:,:)  
     328            ll_found(:,:) = .false.  
     329            DO jj = 1, nlcj  
     330               DO ji = 1, nlci  
     331                 !CDIR NOVECTOR  
     332                  DO jk = ik_ref(ji,jj), mbathy(ji,jj)-1  
     333                     IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN  
     334                        ik_iso(ji,jj)   = jk  
     335                        ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     336                        EXIT  
     337                     ENDIF  
     338                  END DO  
     339               END DO  
     340            END DO  
     341      
     342            ! Use linear interpolation to find depth of mixed layer base where possible  
     343            hmld_kara(:,:) = ppz_ref  
     344            DO jj = 1, jpj  
     345               DO ji = 1, jpi  
     346                  IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN  
     347                     zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     348                     hmld_kara(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     349                  ENDIF  
     350               END DO  
     351            END DO  
     352      
     353            ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     354            ! from the reference density/temperature  
     355      
     356            ! Prevent this section from working on land points  
     357            WHERE( tmask(:,:,1) /= 1.0 )  
     358               ll_found = .true.  
     359            ENDWHERE  
     360      
     361            DO jk = 1, jpk  
     362               ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= & 
     363               & zdelta_T(:,:) 
     364            END DO  
     365      
     366            ! Set default value where interpolation cannot be used (ll_found=false)   
     367            DO jj = 1, jpj  
     368               DO ji = 1, jpi  
     369                  IF( .NOT. ll_found(ji,jj) )  & 
     370                  &   hmld_kara(ji,jj) = fsdept(ji,jj,mbathy(ji,jj))  
     371               END DO  
     372            END DO  
     373      
     374            DO jj = 1, jpj  
     375               DO ji = 1, jpi  
     376                  !CDIR NOVECTOR  
     377                  DO jk = ik_ref(ji,jj)+1, mbathy(ji,jj)  
     378                     IF( ll_found(ji,jj) ) EXIT  
     379                     IF( ll_belowml(ji,jj,jk) ) THEN                 
     380                        zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * & 
     381                        &      SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     382                        zdT  = zT_b - zT(ji,jj,jk-1)                                       
     383                        zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     384                        hmld_kara(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     385                        EXIT                                                    
     386                     ENDIF  
     387                  END DO  
     388               END DO  
     389            END DO  
     390      
     391            hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1)  
     392  
     393            IF(  ln_kara_write25h  ) THEN 
     394               !Double IF required as i_steps not defined if ln_kara_write25h = 
     395               ! FALSE 
     396               IF ( ( MOD( kt, i_steps ) == 0 ) .OR.  kara_25h_init ) THEN 
     397                  hmld_kara_25h = hmld_kara_25h + hmld_kara 
     398                  i_cnt_25h = i_cnt_25h + 1 
     399                  IF ( kara_25h_init ) kara_25h_init = .FALSE. 
     400               ENDIF 
     401            ENDIF 
     402  
     403#if defined key_iomput  
     404            IF( kt /= nit000 ) THEN  
     405               CALL iom_put( "mldkara"  , hmld_kara )    
     406               IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND.  ln_kara_write25h ) & 
     407                  CALL iom_put( "kara25h"  , ( hmld_kara_25h / 25._wp ) ) 
     408            ENDIF 
     409#endif 
     410  
     411         ENDIF 
     412      ENDIF 
     413        
     414   END SUBROUTINE zdf_mxl_kara  
     415    
    128416   !!====================================================================== 
    129417END MODULE zdfmxl 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r7363 r7367  
    8787   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    8888   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
     89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz 
     90   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz 
    8991#if defined key_c1d 
    9092   !                                                                        !!** 1D cfg only  **   ('key_c1d') 
     
    112114         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          & 
    113115#endif 
    114          &      en   (jpi,jpj,jpk) , htau (jpi,jpj)     , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 
     116         &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     &  
     117         &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          & 
     118         &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      ) 
    115119         ! 
    116120      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc ) 
     
    168172      !!---------------------------------------------------------------------- 
    169173      ! 
     174      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
     175         avt (:,:,:) = avt_k (:,:,:)  
     176         avm (:,:,:) = avm_k (:,:,:)  
     177         avmu(:,:,:) = avmu_k(:,:,:)  
     178         avmv(:,:,:) = avmv_k(:,:,:)  
     179      ENDIF  
     180      ! 
    170181      CALL tke_tke      ! now tke (en) 
    171182      ! 
    172183      CALL tke_avn      ! now avt, avm, avmu, avmv 
     184      ! 
     185      avt_k (:,:,:) = avt (:,:,:)  
     186      avm_k (:,:,:) = avm (:,:,:)  
     187      avmu_k(:,:,:) = avmu(:,:,:)  
     188      avmv_k(:,:,:) = avmv(:,:,:)  
    173189      ! 
    174190   END SUBROUTINE zdf_tke 
     
    811827        !                                   ! ------------------- 
    812828        IF(lwp) WRITE(numout,*) '---- tke-rst ----' 
    813         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    ) 
    814         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   ) 
    815         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   ) 
    816         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  ) 
    817         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  ) 
    818         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 
     829        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
     830        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  ) 
     831        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  ) 
     832        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
     833        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 
     834        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  ) 
    819835        ! 
    820836     ENDIF 
Note: See TracChangeset for help on using the changeset viewer.