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 5249 for branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 – NEMO

Ignore:
Timestamp:
2015-04-30T15:21:55+02:00 (10 years ago)
Author:
davestorkey
Message:

UKMO mld_zint branch: implement vertically-interpolated MLD.

File:
1 edited

Legend:

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

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