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 7591 – NEMO

Changeset 7591


Ignore:
Timestamp:
2017-01-20T16:53:40+01:00 (7 years ago)
Author:
kingr
Message:

Addedd Kara ML code to avoid conflicts with branches/UKMO/CO6_KD490_amm7_oper

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r6661 r7591  
    2222   USE timing          ! Timing 
    2323   USE trc_oce, ONLY : lk_offline ! offline flag 
     24   USE lbclnk  
     25   USE eosbn2          ! Equation of state 
    2426 
    2527   IMPLICIT NONE 
     
    3537   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
    3638   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m] 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_kara  !: mixed layer depth of Kara et al   [m] 
    3740 
    3841   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3942   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     43 
     44   ! Namelist variables for  namzdf_karaml  
     45  
     46   LOGICAL   :: ln_kara            ! Logical switch for calculating kara mixed 
     47                                     ! layer 
     48   LOGICAL   :: ln_kara_write25h   ! Logical switch for 25 hour outputs 
     49   INTEGER   :: jpmld_type         ! mixed layer type              
     50   REAL(wp)  :: ppz_ref            ! depth of initial T_ref  
     51   REAL(wp)  :: ppdT_crit          ! Critical temp diff   
     52   REAL(wp)  :: ppiso_frac         ! Fraction of ppdT_crit used  
     53    
     54   !Used for 25h mean 
     55   LOGICAL, PRIVATE :: kara_25h_init = .TRUE.   !Logical used to initalise 25h  
     56                                                !outputs. Necissary, because we need to 
     57                                                !initalise the kara_25h on the zeroth 
     58                                                !timestep (i.e in the nemogcm_init call) 
     59   REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h 
    4060 
    4161   !! * Substitutions 
     
    126146      END DO 
    127147      ! depth of the mixing and mixed layers 
     148 
     149      CALL zdf_mxl_kara( kt ) ! kara MLD 
     150 
    128151      DO jj = 1, jpj 
    129152         DO ji = 1, jpi 
     
    193216   END SUBROUTINE zdf_mxl_tref 
    194217 
     218   SUBROUTINE zdf_mxl_kara( kt )  
     219      !!----------------------------------------------------------------------------------  
     220      !!                    ***  ROUTINE zdf_mxl_kara  ***  
     221      !                                                                         
     222      !   Calculate mixed layer depth according to the definition of           
     223      !   Kara et al, 2000, JGR, 105, 16803.   
     224      !  
     225      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     226      !   density has increased by an amount equivalent to a temperature difference of   
     227      !   0.8C at the surface.  
     228      !  
     229      !   For other values of mld_type the mixed layer is calculated as the depth at   
     230      !   which the temperature differs by 0.8C from the surface temperature.   
     231      !                                                                         
     232      !   Original version: David Acreman                                       
     233      !  
     234      !!----------------------------------------------------------------------------------- 
     235      
     236      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index  
     237  
     238      NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, & 
     239      &                       ppiso_frac, ln_kara_write25h  
     240  
     241      ! Local variables                                                         
     242      REAL, DIMENSION(jpi,jpk) :: ppzdep      ! depth for use in calculating d(rho)  
     243      REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d  !Local version of tsn  
     244      LOGICAL :: ll_found(jpi,jpj)              ! Is T_b to be found by interpolation ?  
     245      LOGICAL :: ll_belowml(jpi,jpj,jpk)        ! Flag points below mixed layer when ll_found=F  
     246      INTEGER :: ji, jj, jk                     ! loop counter  
     247      INTEGER :: ik_ref(jpi,jpj)                ! index of reference level  
     248      INTEGER :: ik_iso(jpi,jpj)                ! index of last uniform temp level  
     249      REAL    :: zT(jpi,jpj,jpk)                ! Temperature or denisty  
     250      REAL    :: zT_ref(jpi,jpj)                ! reference temperature  
     251      REAL    :: zT_b                           ! base temperature  
     252      REAL    :: zdTdz(jpi,jpj,jpk-2)           ! gradient of zT  
     253      REAL    :: zmoddT(jpi,jpj,jpk-2)          ! Absolute temperature difference  
     254      REAL    :: zdz                            ! depth difference  
     255      REAL    :: zdT                            ! temperature difference  
     256      REAL    :: zdelta_T(jpi,jpj)              ! difference critereon  
     257      REAL    :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities 
     258      INTEGER, SAVE :: idt, i_steps 
     259      INTEGER, SAVE :: i_cnt_25h     ! Counter for 25 hour means 
     260      INTEGER :: ios                 ! Local integer output status for namelist read 
     261 
     262      
     263  
     264      !!-------------------------------------------------------------------------------------  
     265  
     266      IF( kt == nit000 ) THEN  
     267         ! Default values  
     268         ln_kara      = .FALSE. 
     269         ln_kara_write25h = .FALSE.  
     270         jpmld_type   = 1      
     271         ppz_ref      = 10.0  
     272         ppdT_crit    = 0.2   
     273         ppiso_frac   = 0.1    
     274         ! Read namelist 
     275         REWIND ( numnam_ref ) 
     276         READ   ( numnam_ref, namzdf_karaml, IOSTAT=ios, ERR= 901 ) 
     277901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in reference namelist', lwp ) 
     278 
     279         REWIND( numnam_cfg )              ! Namelist nam_diadiaop in configuration namelist  3D hourly diagnostics 
     280         READ  ( numnam_cfg,  namzdf_karaml, IOSTAT = ios, ERR = 902 ) 
     281902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist', lwp ) 
     282         IF(lwm) WRITE ( numond, namzdf_karaml ) 
     283  
     284 
     285         WRITE(numout,*) '===== Kara mixed layer ====='  
     286         WRITE(numout,*) 'ln_kara = ',    ln_kara 
     287         WRITE(numout,*) 'jpmld_type = ', jpmld_type  
     288         WRITE(numout,*) 'ppz_ref = ',    ppz_ref  
     289         WRITE(numout,*) 'ppdT_crit = ',  ppdT_crit  
     290         WRITE(numout,*) 'ppiso_frac = ', ppiso_frac 
     291         WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h 
     292         WRITE(numout,*) '============================'  
     293       
     294         IF ( .NOT.ln_kara ) THEN 
     295            WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated"  
     296         ELSE 
     297            IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) ) 
     298            IF ( ln_kara_write25h .AND. kara_25h_init ) THEN 
     299               i_cnt_25h = 0 
     300               IF (.NOT. ALLOCATED(hmld_kara_25h) ) & 
     301               &   ALLOCATE( hmld_kara_25h(jpi,jpj) ) 
     302               hmld_kara_25h = 0._wp 
     303               IF( nacc == 1 ) THEN 
     304                  idt = INT(rdtmin) 
     305               ELSE 
     306                  idt = INT(rdt) 
     307               ENDIF 
     308               IF( MOD( 3600,idt) == 0 ) THEN  
     309                  i_steps = 3600 / idt   
     310               ELSE  
     311                  CALL ctl_stop('STOP', & 
     312                  & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// & 
     313                  & ' = 0 otherwise no hourly values are possible')  
     314               ENDIF   
     315            ENDIF 
     316         ENDIF 
     317      ENDIF 
     318       
     319      IF ( ln_kara ) THEN 
     320          
     321         !set critical ln_kara 
     322         ztsn_2d = tsn(:,:,1,:) 
     323         ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit 
     324      
     325         ! Set the mixed layer depth criterion at each grid point  
     326         ppzdep = 0._wp 
     327         IF( jpmld_type == 1 ) THEN                                          
     328            CALL eos ( tsn(:,:,1,:), & 
     329            &          ppzdep(:,:), zRHO1(:,:) )  
     330            CALL eos ( ztsn_2d(:,:,:), & 
     331            &           ppzdep(:,:), zRHO2(:,:) )  
     332            zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     333            ! RHO from eos (2d version) doesn't calculate north or east halo:  
     334            CALL lbc_lnk( zdelta_T, 'T', 1. )  
     335            zT(:,:,:) = rhop(:,:,:)  
     336         ELSE  
     337            zdelta_T(:,:) = ppdT_crit                       
     338            zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     339         ENDIF  
     340      
     341         ! Calculate the gradient of zT and absolute difference for use later  
     342         DO jk = 1 ,jpk-2  
     343            zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     344            zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     345         END DO  
     346      
     347         ! Find density/temperature at the reference level (Kara et al use 10m).           
     348         ! ik_ref is the index of the box centre immediately above or at the reference level  
     349         ! Find ppz_ref in the array of model level depths and find the ref     
     350         ! density/temperature by linear interpolation.                                    
     351         ik_ref = -1 
     352         DO jk = jpkm1, 2, -1  
     353            WHERE( fsdept(:,:,jk) > ppz_ref )  
     354               ik_ref(:,:) = jk - 1  
     355               zT_ref(:,:) = zT(:,:,jk-1) + & 
     356               &             zdTdz(:,:,jk-1) * ( ppz_ref - fsdept(:,:,jk-1) )  
     357            ENDWHERE  
     358         END DO 
     359         IF ( ANY( ik_ref  < 0 ) .OR. ANY( ik_ref  > jpkm1 ) ) THEN 
     360            CALL ctl_stop( "STOP", & 
     361            & "zdf_mxl_kara: unable to find reference level for kara ML" )  
     362         ELSE 
     363            ! If the first grid box centre is below the reference level then use the  
     364            ! top model level to get zT_ref  
     365            WHERE( fsdept(:,:,1) > ppz_ref )   
     366               zT_ref = zT(:,:,1)  
     367               ik_ref = 1  
     368            ENDWHERE  
     369      
     370            ! Search for a uniform density/temperature region where adjacent levels           
     371            ! differ by less than ppiso_frac * deltaT.                                       
     372            ! ik_iso is the index of the last level in the uniform layer   
     373            ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     374            ik_iso(:,:)   = ik_ref(:,:)  
     375            ll_found(:,:) = .false.  
     376            DO jj = 1, nlcj  
     377               DO ji = 1, nlci  
     378                 !CDIR NOVECTOR  
     379                  DO jk = ik_ref(ji,jj), mbathy(ji,jj)-1  
     380                     IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN  
     381                        ik_iso(ji,jj)   = jk  
     382                        ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     383                        EXIT  
     384                     ENDIF  
     385                  END DO  
     386               END DO  
     387            END DO  
     388      
     389            ! Use linear interpolation to find depth of mixed layer base where possible  
     390            hmld_kara(:,:) = ppz_ref  
     391            DO jj = 1, jpj  
     392               DO ji = 1, jpi  
     393                  IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN  
     394                     zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     395                     hmld_kara(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     396                  ENDIF  
     397               END DO  
     398            END DO  
     399      
     400            ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     401            ! from the reference density/temperature  
     402      
     403            ! Prevent this section from working on land points  
     404            WHERE( tmask(:,:,1) /= 1.0 )  
     405               ll_found = .true.  
     406            ENDWHERE  
     407      
     408            DO jk = 1, jpk  
     409               ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= & 
     410               & zdelta_T(:,:) 
     411            END DO  
     412      
     413            ! Set default value where interpolation cannot be used (ll_found=false)   
     414            DO jj = 1, jpj  
     415               DO ji = 1, jpi  
     416                  IF( .NOT. ll_found(ji,jj) )  & 
     417                  &   hmld_kara(ji,jj) = fsdept(ji,jj,mbathy(ji,jj))  
     418               END DO  
     419            END DO  
     420      
     421            DO jj = 1, jpj  
     422               DO ji = 1, jpi  
     423                  !CDIR NOVECTOR  
     424                  DO jk = ik_ref(ji,jj)+1, mbathy(ji,jj)  
     425                     IF( ll_found(ji,jj) ) EXIT  
     426                     IF( ll_belowml(ji,jj,jk) ) THEN                 
     427                        zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * & 
     428                        &      SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     429                        zdT  = zT_b - zT(ji,jj,jk-1)                                       
     430                        zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     431                        hmld_kara(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     432                        EXIT                                                    
     433                     ENDIF  
     434                  END DO  
     435               END DO  
     436            END DO  
     437      
     438            hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1)  
     439  
     440            IF(  ln_kara_write25h  ) THEN 
     441               !Double IF required as i_steps not defined if ln_kara_write25h = 
     442               ! FALSE 
     443               IF ( ( MOD( kt, i_steps ) == 0 ) .OR.  kara_25h_init ) THEN 
     444                  hmld_kara_25h = hmld_kara_25h + hmld_kara 
     445                  i_cnt_25h = i_cnt_25h + 1 
     446                  IF ( kara_25h_init ) kara_25h_init = .FALSE. 
     447               ENDIF 
     448            ENDIF 
     449  
     450#if defined key_iomput  
     451            IF( kt /= nit000 ) THEN  
     452               CALL iom_put( "mldkara"  , hmld_kara )    
     453               IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND.  ln_kara_write25h ) & 
     454                  CALL iom_put( "kara25h"  , ( hmld_kara_25h / 25._wp ) ) 
     455            ENDIF 
     456#endif 
     457  
     458         ENDIF 
     459      ENDIF 
     460        
     461   END SUBROUTINE zdf_mxl_kara  
     462 
    195463   !!====================================================================== 
    196464END MODULE zdfmxl 
Note: See TracChangeset for help on using the changeset viewer.