Changeset 12585


Ignore:
Timestamp:
2020-03-23T10:58:00+01:00 (9 months ago)
Author:
jcastill
Message:

Update to CO8 version - the previous one crashed after a divide by zero error

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/r12083_mix-lyr_diag/src/OCE/ZDF/zdfmxl.F90

    r12479 r12585  
    2626   PRIVATE 
    2727 
    28    PUBLIC   zdf_mxl_tref  ! called by asminc.F90 
    2928   PUBLIC   zdf_mxl       ! called by zdfphy.F90 
    3029 
     
    8281   END FUNCTION zdf_mxl_alloc 
    8382 
    84    SUBROUTINE zdf_mxl_tref()    
    85       !!----------------------------------------------------------------------    
    86       !!                  ***  ROUTINE zdf_mxl_tref  ***    
    87       !!                       
    88       !! ** Purpose :   Compute the mixed layer depth with temperature criteria.    
    89       !!    
    90       !! ** Method  :   The temperature-defined mixed layer depth is required    
    91       !!                   when assimilating SST in a 2D analysis.     
    92       !!    
    93       !! ** Action  :   hmld_tref    
    94       !!----------------------------------------------------------------------    
    95       !    
    96       INTEGER  ::   ji, jj, jk   ! dummy loop indices    
    97       REAL(wp) ::   t_ref               ! Reference temperature      
    98       REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth      
    99       !!----------------------------------------------------------------------    
    100       !    
    101       ! Initialise array    
    102       IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' )    
    103           
    104       DO jj = 1, jpj      
    105          DO ji = 1, jpi      
    106            hmld_tref(ji,jj)=gdept_n(ji,jj,1  )       
    107            IF(ssmask(ji,jj) > 0.)THEN      
    108              t_ref=tsn(ji,jj,1,jp_tem)     
    109              DO jk=2,jpk      
    110                IF(ssmask(ji,jj)==0.)THEN      
    111                   hmld_tref(ji,jj)=gdept_n(ji,jj,jk )      
    112                   EXIT      
    113                ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN      
    114                   hmld_tref(ji,jj)=gdept_n(ji,jj,jk )      
    115                ELSE      
    116                   EXIT      
    117                ENDIF      
    118              ENDDO      
    119            ENDIF      
    120          ENDDO      
    121       ENDDO    
    122       
    123    END SUBROUTINE zdf_mxl_tref 
    124  
    12583   SUBROUTINE zdf_mxl( kt ) 
    12684      !!---------------------------------------------------------------------- 
     
    210168   END SUBROUTINE zdf_mxl 
    211169 
    212    SUBROUTINE zdf_mxl_zint_mld( sf )   
    213       !!----------------------------------------------------------------------------------   
    214       !!                    ***  ROUTINE zdf_mxl_zint_mld  ***   
    215       !                                                                          
    216       !   Calculate vertically-interpolated mixed layer depth diagnostic.   
    217       !              
    218       !   This routine can calculate the mixed layer depth diagnostic suggested by  
    219       !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate  
    220       !   vertically-interpolated mixed-layer depth diagnostics with other parameter  
    221       !   settings set in the namzdf_mldzint namelist.    
     170   SUBROUTINE zdf_mxl_zint_mld( sf )  
     171      !!----------------------------------------------------------------------------------  
     172      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***  
     173      !                                                                         
     174      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     175      !             
     176      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     177      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     178      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     179      !   settings set in the namzdf_mldzint namelist.   
     180      !  
     181      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     182      !   density has increased by an amount equivalent to a temperature difference of   
     183      !   0.8C at the surface.  
     184      !  
     185      !   For other values of mld_type the mixed layer is calculated as the depth at   
     186      !   which the temperature differs by 0.8C from the surface temperature.   
     187      !                                                                         
     188      !   David Acreman, Daley Calvert                                       
     189      !  
     190      !!-----------------------------------------------------------------------------------  
     191 
     192      TYPE(MXL_ZINT), INTENT(in)  :: sf 
     193 
     194      ! Diagnostic criteria 
     195      INTEGER   :: nn_mld_type   ! mixed layer type      
     196      REAL(wp)  :: rn_zref       ! depth of initial T_ref 
     197      REAL(wp)  :: rn_dT_crit    ! Critical temp diff 
     198      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used 
     199 
     200      ! Local variables 
     201      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value 
     202      INTEGER, DIMENSION(jpi,jpj) :: ikmt          ! number of active tracer levels  
     203      INTEGER, DIMENSION(jpi,jpj) :: ik_ref        ! index of reference level  
     204      INTEGER, DIMENSION(jpi,jpj) :: ik_iso        ! index of last uniform temp level  
     205      REAL, DIMENSION(jpi,jpj,jpk)  :: zT            ! Temperature or density  
     206      REAL, DIMENSION(jpi,jpj)    :: ppzdep        ! depth for use in calculating d(rho)  
     207      REAL, DIMENSION(jpi,jpj)    :: zT_ref        ! reference temperature  
     208      REAL    :: zT_b                                   ! base temperature  
     209      REAL, DIMENSION(jpi,jpj,jpk)  :: zdTdz         ! gradient of zT  
     210      REAL, DIMENSION(jpi,jpj,jpk)  :: zmoddT        ! Absolute temperature difference  
     211      REAL    :: zdz                                    ! depth difference  
     212      REAL    :: zdT                                    ! temperature difference  
     213      REAL, DIMENSION(jpi,jpj)    :: zdelta_T      ! difference critereon  
     214      REAL, DIMENSION(jpi,jpj)    :: zRHO1, zRHO2  ! Densities  
     215      INTEGER :: ji, jj, jk                             ! loop counter  
     216 
     217      !!-------------------------------------------------------------------------------------  
    222218      !   
    223       !   If mld_type=1 the mixed layer depth is calculated as the depth at which the    
    224       !   density has increased by an amount equivalent to a temperature difference of    
    225       !   0.8C at the surface.   
    226       !   
    227       !   For other values of mld_type the mixed layer is calculated as the depth at    
    228       !   which the temperature differs by 0.8C from the surface temperature.    
    229       !                                                                          
    230       !   David Acreman, Daley Calvert                                        
    231       !   
    232       !!-----------------------------------------------------------------------------------   
    233   
    234       TYPE(MXL_ZINT), INTENT(in)  :: sf  
    235   
    236       ! Diagnostic criteria  
    237       INTEGER   :: nn_mld_type   ! mixed layer type       
    238       REAL(wp)  :: rn_zref       ! depth of initial T_ref  
    239       REAL(wp)  :: rn_dT_crit    ! Critical temp diff  
    240       REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used  
    241   
    242       ! Local variables  
    243       REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value  
    244       INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level   
    245       INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level   
    246       REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density   
    247       REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)   
    248       REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature   
    249       REAL    :: zT_b                                   ! base temperature   
    250       REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT   
    251       REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference   
    252       REAL    :: zdz                                    ! depth difference   
    253       REAL    :: zdT                                    ! temperature difference   
    254       REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon   
    255       REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities   
    256       INTEGER :: ji, jj, jk                             ! loop counter   
    257   
    258       !!-------------------------------------------------------------------------------------   
    259       !    
    260       ALLOCATE( ik_ref(jpi,jpj), ik_iso(jpi,jpj), ppzdep(jpi,jpj), zT_ref(jpi,jpj), zdelta_T(jpi,jpj), &  
    261                 zRHO1(jpi,jpj), zRHO2(jpi,jpj), zT(jpi,jpj,jpk), zdTdz(jpi,jpj,jpk),                   & 
    262                 zmoddT(jpi,jpj,jpk), STAT=ji ) 
    263       IF( lk_mpp  )   CALL mpp_sum( 'zdfmxl', ji )  
    264       IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_mld : unable to allocate arrays' )  
    265   
    266       ! Unpack structure  
    267       nn_mld_type = sf%mld_type  
    268       rn_zref     = sf%zref  
    269       rn_dT_crit  = sf%dT_crit  
    270       rn_iso_frac = sf%iso_frac  
    271   
    272       ! Set the mixed layer depth criterion at each grid point   
    273       IF( nn_mld_type == 0 ) THEN  
    274          zdelta_T(:,:) = rn_dT_crit  
     219      ! Unpack structure 
     220      nn_mld_type = sf%mld_type 
     221      rn_zref     = sf%zref 
     222      rn_dT_crit  = sf%dT_crit 
     223      rn_iso_frac = sf%iso_frac 
     224 
     225      ! Set the mixed layer depth criterion at each grid point  
     226      IF( nn_mld_type == 0 ) THEN 
     227         zdelta_T(:,:) = rn_dT_crit 
     228         zT(:,:,:) = rhop(:,:,:) 
     229      ELSE IF( nn_mld_type == 1 ) THEN 
     230         ppzdep(:,:)=0.0  
     231         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     232! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     233! [assumes number of tracers less than number of vertical levels]  
     234         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     235         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     236         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     237         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     238         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     239         CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )  
    275240         zT(:,:,:) = rhop(:,:,:)  
    276       ELSE IF( nn_mld_type == 1 ) THEN  
    277          ppzdep(:,:)=0.0   
    278          call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )   
    279 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST   
    280 ! [assumes number of tracers less than number of vertical levels]   
    281          zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)   
    282          zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit   
    283          CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )   
    284          zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0   
    285          ! RHO from eos (2d version) doesn't calculate north or east halo:   
    286          CALL lbc_lnk( 'zdfmlx', zdelta_T, 'T', 1. )   
    287          zT(:,:,:) = rhop(:,:,:)   
    288       ELSE   
    289          zdelta_T(:,:) = rn_dT_crit                        
    290          zT(:,:,:) = tsn(:,:,:,jp_tem)                             
    291       END IF   
    292   
    293       ! Calculate the gradient of zT and absolute difference for use later   
    294       DO jk = 1 ,jpk-2   
    295          zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)   
    296          zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )   
    297       END DO   
    298   
    299       ! Find density/temperature at the reference level (Kara et al use 10m).            
    300       ! ik_ref is the index of the box centre immediately above or at the reference level   
    301       ! Find rn_zref in the array of model level depths and find the ref      
    302       ! density/temperature by linear interpolation.                                     
    303       DO jk = jpkm1, 2, -1   
    304          WHERE ( gdept_n(:,:,jk) > rn_zref )   
    305            ik_ref(:,:) = jk - 1   
    306            zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )   
    307          END WHERE   
    308       END DO   
    309   
    310       ! If the first grid box centre is below the reference level then use the   
    311       ! top model level to get zT_ref   
    312       WHERE ( gdept_n(:,:,1) > rn_zref )    
    313          zT_ref = zT(:,:,1)   
    314          ik_ref = 1   
    315       END WHERE   
    316   
    317       ! Initialize / reset  
    318       ll_found(:,:) = .false.  
    319   
    320       IF ( rn_iso_frac - zepsilon > 0. ) THEN  
    321          ! Search for a uniform density/temperature region where adjacent levels            
    322          ! differ by less than rn_iso_frac * deltaT.                                        
    323          ! ik_iso is the index of the last level in the uniform layer    
    324          ! ll_found indicates whether the mixed layer depth can be found by interpolation   
    325          ik_iso(:,:)   = ik_ref(:,:)   
    326          DO jj = 1, nlcj   
    327             DO ji = 1, nlci   
    328 !CDIR NOVECTOR   
    329                DO jk = ik_ref(ji,jj), mbkt(ji,jj)-1   
    330                   IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN   
    331                      ik_iso(ji,jj)   = jk   
    332                      ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )   
    333                      EXIT   
    334                   END IF   
    335                END DO   
    336             END DO   
    337          END DO   
    338   
    339          ! Use linear interpolation to find depth of mixed layer base where possible   
    340          hmld_zint(:,:) = rn_zref   
    341          DO jj = 1, jpj   
    342             DO ji = 1, jpi   
    343                IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN   
    344                   zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )   
    345                   hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz   
    346                END IF   
    347             END DO   
    348          END DO   
     241      ELSE  
     242         zdelta_T(:,:) = rn_dT_crit                       
     243         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
    349244      END IF  
    350   
    351       ! If ll_found = .false. then calculate MLD using difference of zdelta_T      
    352       ! from the reference density/temperature   
    353    
    354 ! Prevent this section from working on land points   
    355       WHERE ( tmask(:,:,1) /= 1.0 )   
    356          ll_found = .true.   
    357       END WHERE   
    358    
    359       DO jk=1, jpk   
    360          ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)    
    361       END DO   
    362    
    363 ! Set default value where interpolation cannot be used (ll_found=false)    
    364       DO jj = 1, jpj   
    365          DO ji = 1, jpi   
    366             IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,mbkt(ji,jj))   
    367          END DO   
    368       END DO   
    369   
    370       DO jj = 1, jpj   
    371          DO ji = 1, jpi   
    372 !CDIR NOVECTOR   
    373             DO jk = ik_ref(ji,jj)+1, mbkt(ji,jj)   
    374                IF ( ll_found(ji,jj) ) EXIT   
    375                IF ( ll_belowml(ji,jj,jk) ) THEN                  
    376                   zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )   
    377                   zdT  = zT_b - zT(ji,jj,jk-1)                                        
    378                   zdz  = zdT / zdTdz(ji,jj,jk-1)                                         
    379                   hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz   
    380                   EXIT                                                     
    381                END IF   
    382             END DO   
    383          END DO   
    384       END DO   
    385   
    386       hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)   
    387       !    
    388    END SUBROUTINE zdf_mxl_zint_mld  
    389   
    390    SUBROUTINE zdf_mxl_zint_htc( kt )  
    391       !!----------------------------------------------------------------------  
    392       !!                  ***  ROUTINE zdf_mxl_zint_htc  ***  
    393       !!   
    394       !! ** Purpose :     
    395       !!  
    396       !! ** Method  :     
    397       !!----------------------------------------------------------------------  
    398   
    399       INTEGER, INTENT(in) ::   kt   ! ocean time-step index  
    400   
    401       INTEGER :: ji, jj, jk  
    402       INTEGER :: ikmax  
    403       REAL(wp) :: zc, zcoef  
    404       !  
    405       INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel  
    406       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick  
    407   
    408       !!----------------------------------------------------------------------  
    409   
    410       IF( .NOT. ALLOCATED(ilevel) ) THEN  
    411          ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &  
    412          &         zthick(jpi,jpj), STAT=ji )  
    413          IF( lk_mpp  )   CALL mpp_sum( 'zdfmxl', ji )  
    414          IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )  
    415       ENDIF  
    416   
    417       ! Find last whole model T level above the MLD  
    418       ilevel(:,:)   = 0  
    419       zthick_0(:,:) = 0._wp  
    420   
    421       DO jk = 1, jpkm1    
    422          DO jj = 1, jpj  
    423             DO ji = 1, jpi                      
    424                zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk)  
    425                IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk  
     245 
     246      ! Calculate the gradient of zT and absolute difference for use later  
     247      DO jk = 1 ,jpk-2  
     248         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)  
     249         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     250      END DO  
     251 
     252      ! Find density/temperature at the reference level (Kara et al use 10m).           
     253      ! ik_ref is the index of the box centre immediately above or at the reference level  
     254      ! Find rn_zref in the array of model level depths and find the ref     
     255      ! density/temperature by linear interpolation.                                    
     256      DO jk = jpkm1, 2, -1  
     257         WHERE ( gdept_n(:,:,jk) > rn_zref )  
     258           ik_ref(:,:) = jk - 1  
     259           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )  
     260         END WHERE  
     261      END DO  
     262 
     263      ! If the first grid box centre is below the reference level then use the  
     264      ! top model level to get zT_ref  
     265      WHERE ( gdept_n(:,:,1) > rn_zref )   
     266         zT_ref = zT(:,:,1)  
     267         ik_ref = 1  
     268      END WHERE  
     269 
     270      ! The number of active tracer levels is 1 less than the number of active w levels  
     271      ikmt(:,:) = mbkt(:,:) - 1  
     272 
     273      ! Initialize / reset 
     274      ll_found(:,:) = .false. 
     275 
     276      IF ( rn_iso_frac - zepsilon > 0. ) THEN 
     277         ! Search for a uniform density/temperature region where adjacent levels           
     278         ! differ by less than rn_iso_frac * deltaT.                                       
     279         ! ik_iso is the index of the last level in the uniform layer   
     280         ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     281         ik_iso(:,:)   = ik_ref(:,:)  
     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  
    426292            END DO  
    427293         END DO  
    428          WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)  
    429          WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1)  
     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) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz  
     302               END IF  
     303            END DO  
     304         END DO  
     305      END IF 
     306 
     307      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     308      ! from the reference density/temperature  
     309  
     310! Prevent this section from working on land points  
     311      WHERE ( tmask(:,:,1) /= 1.0 )  
     312         ll_found = .true.  
     313      END WHERE  
     314  
     315      DO jk=1, jpk  
     316         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
    430317      END DO  
    431318  
    432       ! Surface boundary condition  
    433       IF(.NOT.ln_linssh) THEN   ;   zthick(:,:) = 0._wp       ; htc_mld(:,:) = 0._wp                                     
    434       ELSE                      ;   zthick(:,:) = sshn(:,:)   ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)     
    435       ENDIF  
    436   
    437       ! Deepest whole T level above the MLD  
    438       ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )  
    439   
    440       ! Integration down to last whole model T level  
    441       DO jk = 1, ikmax  
    442          DO jj = 1, jpj  
    443             DO ji = 1, jpi  
    444                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel  
    445                zthick(ji,jj) = zthick(ji,jj) + zc  
    446                htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)  
     319! Set default value where interpolation cannot be used (ll_found=false)   
     320      DO jj = 1, jpj  
     321         DO ji = 1, jpi  
     322            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj))  
     323         END DO  
     324      END DO  
     325 
     326      DO jj = 1, jpj  
     327         DO ji = 1, jpi  
     328!CDIR NOVECTOR  
     329            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     330               IF ( ll_found(ji,jj) ) EXIT  
     331               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     332                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     333                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     334                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     335                  hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz  
     336                  EXIT                                                    
     337               END IF  
    447338            END DO  
    448339         END DO  
    449340      END DO  
    450   
    451       ! Subsequent partial T level  
    452       zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD  
    453   
    454       DO jj = 1, jpj  
    455          DO ji = 1, jpi  
    456             htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &   
    457       &                      * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)  
    458          END DO  
    459       END DO  
    460   
    461       WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)  
    462   
    463       ! Convert to heat content  
    464       zcoef = rau0 * rcp  
    465       htc_mld(:,:) = zcoef * htc_mld(:,:)  
    466   
    467    END SUBROUTINE zdf_mxl_zint_htc  
    468   
    469    SUBROUTINE zdf_mxl_zint( kt )  
    470       !!----------------------------------------------------------------------  
    471       !!                  ***  ROUTINE zdf_mxl_zint  ***  
    472       !!   
    473       !! ** Purpose :     
     341 
     342      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     343      !   
     344   END SUBROUTINE zdf_mxl_zint_mld 
     345 
     346   SUBROUTINE zdf_mxl_zint_htc( kt ) 
     347      !!---------------------------------------------------------------------- 
     348      !!                  ***  ROUTINE zdf_mxl_zint_htc  *** 
    474349      !!  
    475       !! ** Method  :     
    476       !!----------------------------------------------------------------------  
    477   
    478       INTEGER, INTENT(in) ::   kt   ! ocean time-step index  
    479   
    480       INTEGER :: ios  
    481       INTEGER :: jn  
    482   
    483       INTEGER :: nn_mld_diag = 0    ! number of diagnostics  
    484  
    485       INTEGER :: i_steps            ! no of timesteps per hour  
    486       INTEGER :: ierror             ! logical error message     
    487  
    488       REAL(wp) :: zdt               ! timestep variable 
    489   
    490       CHARACTER(len=1) :: cmld  
    491   
    492       TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5  
    493       TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags  
    494   
    495       NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5  
    496   
    497       !!----------------------------------------------------------------------  
    498         
    499       IF( kt == nit000 ) THEN  
    500          REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist   
    501          READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)  
    502 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' )  
    503   
    504          REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist   
    505          READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )  
    506 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' )  
    507          IF(lwm) WRITE ( numond, namzdf_mldzint )  
    508   
    509          IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )  
    510   
    511          mld_diags(1) = sn_mld1  
    512          mld_diags(2) = sn_mld2  
    513          mld_diags(3) = sn_mld3  
    514          mld_diags(4) = sn_mld4  
    515          mld_diags(5) = sn_mld5  
    516   
    517          IF( nn_mld_diag > 0 ) THEN  
    518             WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'  
    519             WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'  
    520             DO jn = 1, nn_mld_diag  
    521                WRITE(numout,*) 'MLD criterion',jn,':'  
    522                WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type  
    523                WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref  
    524                WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit  
    525                WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac  
    526             END DO  
    527             WRITE(numout,*) '===================================================================='  
    528          ENDIF  
    529       ENDIF  
    530   
    531       IF( nn_mld_diag > 0 ) THEN  
    532          DO jn = 1, nn_mld_diag  
    533             WRITE(cmld,'(I1)') jn  
    534             IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN  
    535                CALL zdf_mxl_zint_mld( mld_diags(jn) )  
    536   
    537                IF( iom_use( "mldzint_"//cmld ) ) THEN  
    538                   CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )  
    539                ENDIF  
    540   
    541                IF( iom_use( "mldhtc_"//cmld ) )  THEN  
    542                   CALL zdf_mxl_zint_htc( kt )  
    543                   CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) )  
    544                ENDIF  
    545  
    546                IF( iom_use( "mldzint25h_"//cmld ) ) THEN  
    547                   IF( .NOT. mld_25h_write ) mld_25h_write = .TRUE.  
    548                   zdt = rdt  
    549                   IF( MOD( 3600,INT(zdt) ) == 0 ) THEN  
    550                      i_steps = 3600/INT(zdt)  
    551                   ELSE  
    552                      CALL ctl_stop('STOP', 'zdf_mxl_zint 25h: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')  
    553                   ENDIF  
    554                   IF( ( mld_25h_init ) .OR. ( kt == nit000 ) ) THEN  
    555                      i_cnt_25h = 1  
    556                      IF( .NOT. ALLOCATED(hmld_zint_25h) ) THEN  
    557                         ALLOCATE( hmld_zint_25h(jpi,jpj,nn_mld_diag), STAT=ierror )  
    558                         IF( ierror > 0 )  CALL ctl_stop( 'zdf_mxl_zint 25h: unable to allocate hmld_zint_25h' )     
    559                      ENDIF  
    560                      hmld_zint_25h(:,:,jn) = hmld_zint(:,:)  
    561                   ENDIF  
    562                   IF( MOD( kt, i_steps ) == 0 .AND.  kt .NE. nn_it000 ) THEN  
    563                      hmld_zint_25h(:,:,jn) = hmld_zint_25h(:,:,jn) + hmld_zint(:,:)  
    564                   ENDIF  
    565                   IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN  
    566                      CALL iom_put( "mldzint25h_"//cmld , hmld_zint_25h(:,:,jn) / 25._wp   )  
    567                   ENDIF  
     350      !! ** Purpose :    
     351      !! 
     352      !! ** Method  :    
     353      !!---------------------------------------------------------------------- 
     354 
     355      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     356 
     357      INTEGER :: ji, jj, jk 
     358      INTEGER :: ikmax 
     359      REAL(wp) :: zc, zcoef 
     360      ! 
     361      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel 
     362      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick 
     363 
     364      !!---------------------------------------------------------------------- 
     365 
     366      IF( .NOT. ALLOCATED(ilevel) ) THEN 
     367         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 
     368         &         zthick(jpi,jpj), STAT=ji ) 
     369         IF( lk_mpp  )   CALL mpp_sum( 'zdfmxl', ji ) 
     370         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 
     371      ENDIF 
     372 
     373      ! Find last whole model T level above the MLD 
     374      ilevel(:,:)   = 0 
     375      zthick_0(:,:) = 0._wp 
     376 
     377      DO jk = 1, jpkm1   
     378         DO jj = 1, jpj 
     379            DO ji = 1, jpi                     
     380               zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 
     381               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk 
     382            END DO 
     383         END DO 
     384         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
     385         WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1) 
     386      END DO 
     387 
     388      ! Surface boundary condition 
     389      IF( ln_linssh ) THEN  ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)    
     390      ELSE                  ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                    
     391      ENDIF 
     392 
     393      ! Deepest whole T level above the MLD 
     394      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
     395 
     396      ! Integration down to last whole model T level 
     397      DO jk = 1, ikmax 
     398         DO jj = 1, jpj 
     399            DO ji = 1, jpi 
     400               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
     401               zthick(ji,jj) = zthick(ji,jj) + zc 
     402               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     403            END DO 
     404         END DO 
     405      END DO 
     406 
     407      ! Subsequent partial T level 
     408      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD 
     409 
     410      DO jj = 1, jpj 
     411         DO ji = 1, jpi 
     412            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
     413      &                      * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
     414         END DO 
     415      END DO 
     416 
     417      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 
     418 
     419      ! Convert to heat content 
     420      zcoef = rau0 * rcp 
     421      htc_mld(:,:) = zcoef * htc_mld(:,:) 
     422 
     423   END SUBROUTINE zdf_mxl_zint_htc 
     424 
     425   SUBROUTINE zdf_mxl_zint( kt ) 
     426      !!---------------------------------------------------------------------- 
     427      !!                  ***  ROUTINE zdf_mxl_zint  *** 
     428      !!  
     429      !! ** Purpose :    
     430      !! 
     431      !! ** Method  :    
     432      !!---------------------------------------------------------------------- 
     433 
     434      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     435 
     436      INTEGER :: ios 
     437      INTEGER :: jn 
     438 
     439      INTEGER :: nn_mld_diag = 0    ! number of diagnostics 
     440 
     441      CHARACTER(len=1) :: cmld 
     442 
     443      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     444      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags 
     445 
     446      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     447 
     448      !!---------------------------------------------------------------------- 
     449       
     450      IF( kt == nit000 ) THEN 
     451         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     452         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     453901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' ) 
     454 
     455         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     456         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     457902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' ) 
     458         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     459 
     460         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 
     461 
     462         mld_diags(1) = sn_mld1 
     463         mld_diags(2) = sn_mld2 
     464         mld_diags(3) = sn_mld3 
     465         mld_diags(4) = sn_mld4 
     466         mld_diags(5) = sn_mld5 
     467 
     468         IF( lwp .AND. (nn_mld_diag > 0) ) THEN 
     469            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
     470            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     471            DO jn = 1, nn_mld_diag 
     472               WRITE(numout,*) 'MLD criterion',jn,':' 
     473               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
     474               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
     475               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
     476               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     477            END DO 
     478            WRITE(numout,*) '====================================================================' 
     479         ENDIF 
     480      ENDIF 
     481 
     482      IF( nn_mld_diag > 0 ) THEN 
     483         DO jn = 1, nn_mld_diag 
     484            WRITE(cmld,'(I1)') jn 
     485            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 
     486               CALL zdf_mxl_zint_mld( mld_diags(jn) ) 
     487 
     488               IF( iom_use( "mldzint_"//cmld ) ) THEN 
     489                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 
    568490               ENDIF 
    569491 
    570             ENDIF  
    571          END DO  
    572  
    573          IF(  mld_25h_write  ) THEN  
    574             IF( ( MOD( kt, i_steps ) == 0 ) .OR.  mld_25h_init ) THEN  
    575                IF (lwp) THEN  
    576                   WRITE(numout,*) 'zdf_mxl_zint (25h) : Summed the following number of hourly values so far',i_cnt_25h  
    577           ENDIF  
    578                i_cnt_25h = i_cnt_25h + 1  
    579                IF( mld_25h_init ) mld_25h_init = .FALSE.  
    580             ENDIF  
    581             IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN  
    582                i_cnt_25h = 1   
    583                DO jn = 1, nn_mld_diag   
    584                      hmld_zint_25h(:,:,jn) = hmld_zint(:,:)   
    585                ENDDO 
    586             ENDIF  
    587          ENDIF 
    588  
    589       ENDIF  
    590   
     492               IF( iom_use( "mldhtc_"//cmld ) )  THEN 
     493                  CALL zdf_mxl_zint_htc( kt ) 
     494                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   ) 
     495               ENDIF 
     496            ENDIF 
     497         END DO 
     498      ENDIF 
     499 
    591500   END SUBROUTINE zdf_mxl_zint 
    592501 
Note: See TracChangeset for help on using the changeset viewer.