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 6491 for branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Ignore:
Timestamp:
2016-04-21T18:15:17+02:00 (8 years ago)
Author:
davestorkey
Message:

Commit remaining changes to UKMO/r5518_GO6_package branch from following branches:
UKMO/dev_r5021_nn_etau_revision@6238
UKMO/dev_r5107_mld_zint@5534
UKMO/dev_r5107_eorca025_closea@6390
UKMO/restart_datestamp@5539
UKMO/icebergs_latent_heat@5821
UKMO/icebergs_restart_single_file_corrected@6480
UKMO/product_diagnostics@5971
UKMO/antarctic_partial_slip@5961
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5961 cf. r5958 of /branches/UKMO/antarctic_partial_slip/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6349 cf. r5962 of /branches/UKMO/product_diagnostics/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6480 cf. r6479 of /branches/UKMO/icebergs_restart_single_file_corrected/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5986 cf. r5852 of /branches/UKMO/icebergs_restart_single_file/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5821 cf. r5808 of /branches/UKMO/icebergs_latent_heat/NEMOGCM@6490

Location:
branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r6487 r6491  
    1818   USE phycst          ! physical constants 
    1919   USE iom             ! I/O library 
     20   USE eosbn2          ! for zdf_mxl_zint 
    2021   USE lib_mpp         ! MPP library 
    2122   USE wrk_nemo        ! work arrays 
     
    3334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3435   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]  
     37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?  
     38   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3539 
    3640   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3741   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     42 
     43   ! Namelist variables for  namzdf_mldzint 
     44   INTEGER          :: nn_mld_type         ! mixed layer type             
     45   REAL(wp)         :: rn_zref            ! depth of initial T_ref 
     46   REAL(wp)         :: rn_dT_crit          ! Critical temp diff  
     47   REAL(wp)         :: rn_iso_frac         ! Fraction of rn_dT_crit used  
    3848 
    3949   !! * Substitutions 
     
    5262      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5363      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    54          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     64         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       & 
     65        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    5566         ! 
    5667         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
     
    90101      CALL wrk_alloc( jpi,jpj, imld ) 
    91102 
    92       IF( kt == nit000 ) THEN 
     103      IF( kt <= nit000 ) THEN 
    93104         IF(lwp) WRITE(numout,*) 
    94105         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' 
     
    134145      END DO 
    135146      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode 
     147      IF( kt >= nit000 ) THEN               ! workaround for calls before SOMETHING reads the XIOS namelist 
    136148         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth 
    137149         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth 
    138150      ENDIF 
     151      ENDIF 
    139152       
     153      ! Vertically-interpolated mixed-layer depth diagnostic 
     154      IF( iom_use( "mldzint" ) ) THEN 
     155         CALL zdf_mxl_zint( kt ) 
     156         CALL iom_put( "mldzint" , hmld_zint ) 
     157      ENDIF 
     158 
    140159      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    141160      ! 
     
    145164      ! 
    146165   END SUBROUTINE zdf_mxl 
     166 
     167   SUBROUTINE zdf_mxl_zint( kt )  
     168      !!----------------------------------------------------------------------------------  
     169      !!                    ***  ROUTINE zdf_mxl_zint  ***  
     170      !                                                                         
     171      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     172      !             
     173      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     174      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     175      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     176      !   settings set in the namzdf_mldzint namelist.   
     177      !  
     178      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     179      !   density has increased by an amount equivalent to a temperature difference of   
     180      !   0.8C at the surface.  
     181      !  
     182      !   For other values of mld_type the mixed layer is calculated as the depth at   
     183      !   which the temperature differs by 0.8C from the surface temperature.   
     184      !                                                                         
     185      !   David Acreman, Daley Calvert                                       
     186      !  
     187      !!-----------------------------------------------------------------------------------  
     188 
     189      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     190      ! 
     191      ! Local variables 
     192      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels  
     193      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level  
     194      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level  
     195      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density  
     196      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)  
     197      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature  
     198      REAL    :: zT_b                                   ! base temperature  
     199      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT  
     200      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference  
     201      REAL    :: zdz                                    ! depth difference  
     202      REAL    :: zdT                                    ! temperature difference  
     203      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon  
     204      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities  
     205      INTEGER :: ji, jj, jk                             ! loop counter  
     206      INTEGER :: ios 
     207 
     208      NAMELIST/namzdf_mldzint/ nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac 
     209 
     210      !!-------------------------------------------------------------------------------------  
     211      !   
     212      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     213      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     214      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT )  
     215  
     216      IF( kt == nit000 ) THEN 
     217         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     218         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     219901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 
     220 
     221         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     222         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     223902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 
     224         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     225 
     226         WRITE(numout,*) '===== Vertically-interpolated mixed layer =====' 
     227         WRITE(numout,*) 'nn_mld_type = ',nn_mld_type 
     228         WRITE(numout,*) 'rn_zref = ',rn_zref 
     229         WRITE(numout,*) 'rn_dT_crit = ',rn_dT_crit 
     230         WRITE(numout,*) 'rn_iso_frac = ',rn_iso_frac 
     231         WRITE(numout,*) '===============================================' 
     232      ENDIF 
     233  
     234      ! Set the mixed layer depth criterion at each grid point  
     235      IF (nn_mld_type == 1) THEN                                          
     236         ppzdep(:,:)=0.0  
     237         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     238! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     239! [assumes number of tracers less than number of vertical levels]  
     240         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     241         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     242         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     243         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     244         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     245         CALL lbc_lnk( zdelta_T, 'T', 1. )  
     246         zT(:,:,:) = rhop(:,:,:)  
     247      ELSE  
     248         zdelta_T(:,:) = rn_dT_crit                       
     249         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     250      END IF  
     251 
     252      ! Calculate the gradient of zT and absolute difference for use later  
     253      DO jk = 1 ,jpk-2  
     254         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     255         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     256      END DO  
     257 
     258      ! Find density/temperature at the reference level (Kara et al use 10m).           
     259      ! ik_ref is the index of the box centre immediately above or at the reference level  
     260      ! Find rn_zref in the array of model level depths and find the ref     
     261      ! density/temperature by linear interpolation.                                    
     262      DO jk = jpkm1, 2, -1  
     263         WHERE ( fsdept(:,:,jk) > rn_zref )  
     264           ik_ref(:,:) = jk - 1  
     265           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) )  
     266         END WHERE  
     267      END DO  
     268 
     269      ! If the first grid box centre is below the reference level then use the  
     270      ! top model level to get zT_ref  
     271      WHERE ( fsdept(:,:,1) > rn_zref )   
     272         zT_ref = zT(:,:,1)  
     273         ik_ref = 1  
     274      END WHERE  
     275 
     276      ! The number of active tracer levels is 1 less than the number of active w levels  
     277      ikmt(:,:) = mbathy(:,:) - 1  
     278 
     279      ! Search for a uniform density/temperature region where adjacent levels           
     280      ! differ by less than rn_iso_frac * deltaT.                                       
     281      ! ik_iso is the index of the last level in the uniform layer   
     282      ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     283      ik_iso(:,:)   = ik_ref(:,:)  
     284      ll_found(:,:) = .false.  
     285      DO jj = 1, nlcj  
     286         DO ji = 1, nlci  
     287!CDIR NOVECTOR  
     288            DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     289               IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     290                  ik_iso(ji,jj)   = jk  
     291                  ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     292                  EXIT  
     293               END IF  
     294            END DO  
     295         END DO  
     296      END DO  
     297 
     298      ! Use linear interpolation to find depth of mixed layer base where possible  
     299      hmld_zint(:,:) = rn_zref  
     300      DO jj = 1, jpj  
     301         DO ji = 1, jpi  
     302            IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     303               zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     304               hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     305            END IF  
     306         END DO  
     307      END DO  
     308 
     309      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     310      ! from the reference density/temperature  
     311  
     312! Prevent this section from working on land points  
     313      WHERE ( tmask(:,:,1) /= 1.0 )  
     314         ll_found = .true.  
     315      END WHERE  
     316  
     317      DO jk=1, jpk  
     318         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     319      END DO  
     320  
     321! Set default value where interpolation cannot be used (ll_found=false)   
     322      DO jj = 1, jpj  
     323         DO ji = 1, jpi  
     324            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj))  
     325         END DO  
     326      END DO  
     327 
     328      DO jj = 1, jpj  
     329         DO ji = 1, jpi  
     330!CDIR NOVECTOR  
     331            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     332               IF ( ll_found(ji,jj) ) EXIT  
     333               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     334                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     335                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     336                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     337                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     338                  EXIT                                                    
     339               END IF  
     340            END DO  
     341         END DO  
     342      END DO  
     343 
     344      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     345      !   
     346      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     347      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     348      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT )  
     349      !  
     350   END SUBROUTINE zdf_mxl_zint 
    147351 
    148352   !!====================================================================== 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r6487 r6491  
    8383   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1) 
    8484   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
     85   REAL(wp) ::   rn_c      ! fraction of TKE added within the mixed layer by nn_etau 
    8586   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8687   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     
    8889   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
    8990   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m] 
     91   REAL(wp) ::   rhtau                     ! coefficient to relate MLD to htau when nn_htau == 2 
    9092   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3) 
    9193   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    9294 
    9395   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
     96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
     97   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    9498   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    9599#if defined key_c1d 
     
    114118      !!---------------------------------------------------------------------- 
    115119      ALLOCATE(                                                                    & 
     120         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &       
    116121#if defined key_c1d 
    117122         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          & 
     
    420425      END DO 
    421426 
     427      !                                 ! Save TKE prior to nn_etau addition   
     428      e_niw(:,:,:) = en(:,:,:)   
     429      !   
    422430      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    423431      !                            !  TKE due to surface and internal wave breaking 
    424432      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     433      IF( nn_htau == 2 ) THEN           !* mixed-layer depth dependant length scale 
     434         DO jj = 2, jpjm1 
     435            DO ji = fs_2, fs_jpim1   ! vector opt. 
     436               htau(ji,jj) = rhtau * hmlp(ji,jj) 
     437            END DO 
     438         END DO 
     439      ENDIF 
     440#if defined key_iomput 
     441      ! 
     442      CALL iom_put( "htau", htau(:,:) )  ! Check htau (even if constant in time) 
     443#endif 
     444      ! 
    425445      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    426446         DO jk = 2, jpkm1 
     
    457477            END DO 
    458478         END DO 
     479      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up) 
     480         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau  
     481            DO jj = 2, jpjm1 
     482               DO ji = fs_2, fs_jpim1   ! vector opt. 
     483                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 
     484               END DO 
     485            END DO 
     486         ENDIF 
     487         DO jk = 2, jpkm1 
     488            DO jj = 2, jpjm1 
     489               DO ji = fs_2, fs_jpim1   ! vector opt. 
     490                  en(ji,jj,jk) = en(ji,jj,jk) + efr(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     491                     &                                                   * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     492               END DO 
     493            END DO 
     494         END DO 
    459495      ENDIF 
    460496      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     497      ! 
     498      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term   
     499         DO jj = 2, jpjm1   
     500            DO ji = fs_2, fs_jpim1   ! vector opt.   
     501               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)   
     502            END DO   
     503         END DO   
     504      END DO   
     505      !   
     506      CALL lbc_lnk( e_niw, 'W', 1. )   
    461507      ! 
    462508      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
     
    722768         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   & 
    723769         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   & 
    724          &                 nn_etau , nn_htau  , rn_efr    
    725       !!---------------------------------------------------------------------- 
    726       ! 
     770         &                 nn_etau , nn_htau  , rn_efr , rn_c    
     771      !!---------------------------------------------------------------------- 
     772 
    727773      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy 
    728774      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) 
     
    757803         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau 
    758804         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr 
     805         WRITE(numout,*) '      fraction of TKE added within the mixed layer by nn_etau rn_c    = ', rn_c 
    759806         WRITE(numout,*) 
    760807         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri 
     
    767814      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
    768815      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    769       IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     816      IF( nn_htau < 0  .OR.  nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0 to 5    ' ) 
    770817      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    771818 
     
    780827      ENDIF 
    781828 
     829      IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN 
     830          ierr = zdf_mxl_alloc() 
     831          nmln(:,:) = nlb10           ! Initialization of nmln 
     832      ENDIF 
     833 
    782834      !                               !* depth of penetration of surface tke 
    783835      IF( nn_etau /= 0 ) THEN       
     836         htau(:,:) = 0._wp 
    784837         SELECT CASE( nn_htau )             ! Choice of the depth of penetration 
    785838         CASE( 0 )                                 ! constant depth penetration (here 10 meters) 
     
    787840         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees 
    788841            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )             
     842         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer 
     843            rhtau = -1._wp / LOG( 1._wp - rn_c ) 
     844         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees 
     845            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   ) 
     846         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south 
     847            DO jj = 2, jpjm1 
     848               DO ji = fs_2, fs_jpim1   ! vector opt. 
     849                  IF( gphit(ji,jj) <= 0._wp ) THEN 
     850                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     851                  ELSE 
     852                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     853                  ENDIF 
     854               END DO 
     855            END DO 
     856         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south, 
     857            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south 
     858               DO ji = fs_2, fs_jpim1   ! vector opt. 
     859                  IF( gphit(ji,jj) <= -30._wp ) THEN 
     860                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   ) 
     861                  ELSE 
     862                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     863                  ENDIF 
     864               END DO 
     865            END DO 
    789866         END SELECT 
     867         ! 
     868         IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN            ! efr dependant on constant htau 
     869            DO jj = 2, jpjm1 
     870               DO ji = fs_2, fs_jpim1   ! vector opt. 
     871                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 
     872               END DO 
     873            END DO 
     874         ENDIF 
    790875      ENDIF 
    791876      !                               !* set vertical eddy coef. to the background value 
Note: See TracChangeset for help on using the changeset viewer.