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

Ignore:
Timestamp:
2018-10-29T13:03:40+01:00 (5 years ago)
Author:
kingr
Message:

Merged AMM15_v3_6_STABLE_package_collate@10237

File:
1 edited

Legend:

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

    r10248 r10249  
    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 
     
    2627   PRIVATE 
    2728 
     29   PUBLIC   zdf_mxl_tref  ! called by asminc.F90 
    2830   PUBLIC   zdf_mxl       ! called by step.F90 
    29    PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init 
    30  
     31 
     32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m] 
    3133   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP) 
    3234   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    3335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3436   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     37   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]  
     38   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint 
     39   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?  
     40   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3541 
    3642   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3743   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     44 
     45   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs 
     46      INTEGER   :: mld_type   ! mixed layer type      
     47      REAL(wp)  :: zref       ! depth of initial T_ref 
     48      REAL(wp)  :: dT_crit    ! Critical temp diff 
     49      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit used 
     50   END TYPE MXL_ZINT 
     51 
     52!Used for 25h mean 
     53   LOGICAL, PRIVATE :: mld_25h_init = .TRUE.    !Logical used to initalise 25h 
     54                                                !outputs. Necassary, because we need to 
     55                                                !initalise the mld_25h on the zeroth 
     56                                                !timestep (i.e in the nemogcm_init call) 
     57   LOGICAL, PRIVATE :: mld_25h_write = .FALSE.  !Logical confirm 25h calculating/processing 
     58 
     59   INTEGER, SAVE :: i_cnt_25h                   ! Counter for 25 hour means 
     60 
     61   REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   hmld_zint_25h 
    3862 
    3963   !! * Substitutions 
     
    5276      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5377      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    54          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     78         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       & 
     79        &          htc_mld(jpi,jpj),                                                                    & 
     80        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    5581         ! 
     82         ALLOCATE(hmld_tref(jpi,jpj)) 
    5683         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
    5784         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.') 
     
    6087   END FUNCTION zdf_mxl_alloc 
    6188 
     89   SUBROUTINE zdf_mxl_tref()  
     90      !!----------------------------------------------------------------------  
     91      !!                  ***  ROUTINE zdf_mxl_tref  ***  
     92      !!                     
     93      !! ** Purpose :   Compute the mixed layer depth with temperature criteria.  
     94      !!  
     95      !! ** Method  :   The temperature-defined mixed layer depth is required  
     96      !!                   when assimilating SST in a 2D analysis.   
     97      !!  
     98      !! ** Action  :   hmld_tref  
     99      !!----------------------------------------------------------------------  
     100      !  
     101      INTEGER  ::   ji, jj, jk   ! dummy loop indices  
     102      REAL(wp) ::   t_ref               ! Reference temperature    
     103      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth    
     104      !!----------------------------------------------------------------------  
     105      !  
     106      ! Initialise array  
     107      IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' )  
     108        
     109      !For the AMM model assimiation uses a temperature based mixed layer depth    
     110      !This is defined here    
     111      DO jj = 1, jpj    
     112         DO ji = 1, jpi    
     113           hmld_tref(ji,jj)=fsdept(ji,jj,1  )     
     114           IF(ssmask(ji,jj) > 0.)THEN    
     115             t_ref=tsn(ji,jj,1,jp_tem)   
     116             DO jk=2,jpk    
     117               IF(ssmask(ji,jj)==0.)THEN    
     118                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )    
     119                  EXIT    
     120               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN    
     121                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )    
     122               ELSE    
     123                  EXIT    
     124               ENDIF    
     125             ENDDO    
     126           ENDIF    
     127         ENDDO    
     128      ENDDO  
     129    
     130   END SUBROUTINE zdf_mxl_tref  
    62131 
    63132   SUBROUTINE zdf_mxl( kt ) 
     
    128197            iikn = nmln(ji,jj) 
    129198            imkt = mikt(ji,jj) 
    130             hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Turbocline depth  
    131             hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Mixed layer depth 
    132             hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
     199            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! Turbocline depth  
     200            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth 
     201            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
    133202         END DO 
    134203      END DO 
     
    138207      ENDIF 
    139208       
     209      ! Vertically-interpolated mixed-layer depth diagnostic 
     210      CALL zdf_mxl_zint( kt ) 
     211 
    140212      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    141213      ! 
     
    145217      ! 
    146218   END SUBROUTINE zdf_mxl 
     219 
     220   SUBROUTINE zdf_mxl_zint_mld( sf )  
     221      !!----------------------------------------------------------------------------------  
     222      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***  
     223      !                                                                         
     224      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     225      !             
     226      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     227      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     228      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     229      !   settings set in the namzdf_mldzint namelist.   
     230      !  
     231      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     232      !   density has increased by an amount equivalent to a temperature difference of   
     233      !   0.8C at the surface.  
     234      !  
     235      !   For other values of mld_type the mixed layer is calculated as the depth at   
     236      !   which the temperature differs by 0.8C from the surface temperature.   
     237      !                                                                         
     238      !   David Acreman, Daley Calvert                                       
     239      !  
     240      !!-----------------------------------------------------------------------------------  
     241 
     242      TYPE(MXL_ZINT), INTENT(in)  :: sf 
     243 
     244      ! Diagnostic criteria 
     245      INTEGER   :: nn_mld_type   ! mixed layer type      
     246      REAL(wp)  :: rn_zref       ! depth of initial T_ref 
     247      REAL(wp)  :: rn_dT_crit    ! Critical temp diff 
     248      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used 
     249 
     250      ! Local variables 
     251      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value 
     252      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels  
     253      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level  
     254      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level  
     255      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density  
     256      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)  
     257      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature  
     258      REAL    :: zT_b                                   ! base temperature  
     259      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT  
     260      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference  
     261      REAL    :: zdz                                    ! depth difference  
     262      REAL    :: zdT                                    ! temperature difference  
     263      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon  
     264      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities  
     265      INTEGER :: ji, jj, jk                             ! loop counter  
     266 
     267      !!-------------------------------------------------------------------------------------  
     268      !   
     269      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     270      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     271      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT )  
     272 
     273      ! Unpack structure 
     274      nn_mld_type = sf%mld_type 
     275      rn_zref     = sf%zref 
     276      rn_dT_crit  = sf%dT_crit 
     277      rn_iso_frac = sf%iso_frac 
     278 
     279      ! Set the mixed layer depth criterion at each grid point  
     280      IF( nn_mld_type == 0 ) THEN 
     281         zdelta_T(:,:) = rn_dT_crit 
     282         zT(:,:,:) = rhop(:,:,:) 
     283      ELSE IF( nn_mld_type == 1 ) THEN 
     284         ppzdep(:,:)=0.0  
     285         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     286! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     287! [assumes number of tracers less than number of vertical levels]  
     288         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     289         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     290         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     291         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     292         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     293         CALL lbc_lnk( zdelta_T, 'T', 1. )  
     294         zT(:,:,:) = rhop(:,:,:)  
     295      ELSE  
     296         zdelta_T(:,:) = rn_dT_crit                       
     297         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     298      END IF  
     299 
     300      ! Calculate the gradient of zT and absolute difference for use later  
     301      DO jk = 1 ,jpk-2  
     302         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     303         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     304      END DO  
     305 
     306      ! Find density/temperature at the reference level (Kara et al use 10m).           
     307      ! ik_ref is the index of the box centre immediately above or at the reference level  
     308      ! Find rn_zref in the array of model level depths and find the ref     
     309      ! density/temperature by linear interpolation.                                    
     310      DO jk = jpkm1, 2, -1  
     311         WHERE ( fsdept(:,:,jk) > rn_zref )  
     312           ik_ref(:,:) = jk - 1  
     313           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) )  
     314         END WHERE  
     315      END DO  
     316 
     317      ! If the first grid box centre is below the reference level then use the  
     318      ! top model level to get zT_ref  
     319      WHERE ( fsdept(:,:,1) > rn_zref )   
     320         zT_ref = zT(:,:,1)  
     321         ik_ref = 1  
     322      END WHERE  
     323 
     324      ! The number of active tracer levels is 1 less than the number of active w levels  
     325      ikmt(:,:) = mbathy(:,:) - 1  
     326 
     327      ! Initialize / reset 
     328      ll_found(:,:) = .false. 
     329 
     330      IF ( rn_iso_frac - zepsilon > 0. ) THEN 
     331         ! Search for a uniform density/temperature region where adjacent levels           
     332         ! differ by less than rn_iso_frac * deltaT.                                       
     333         ! ik_iso is the index of the last level in the uniform layer   
     334         ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     335         ik_iso(:,:)   = ik_ref(:,:)  
     336         DO jj = 1, nlcj  
     337            DO ji = 1, nlci  
     338!CDIR NOVECTOR  
     339               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     340                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     341                     ik_iso(ji,jj)   = jk  
     342                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     343                     EXIT  
     344                  END IF  
     345               END DO  
     346            END DO  
     347         END DO  
     348 
     349         ! Use linear interpolation to find depth of mixed layer base where possible  
     350         hmld_zint(:,:) = rn_zref  
     351         DO jj = 1, jpj  
     352            DO ji = 1, jpi  
     353               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     354                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     355                  hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     356               END IF  
     357            END DO  
     358         END DO  
     359      END IF 
     360 
     361      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     362      ! from the reference density/temperature  
     363  
     364! Prevent this section from working on land points  
     365      WHERE ( tmask(:,:,1) /= 1.0 )  
     366         ll_found = .true.  
     367      END WHERE  
     368  
     369      DO jk=1, jpk  
     370         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     371      END DO  
     372  
     373! Set default value where interpolation cannot be used (ll_found=false)   
     374      DO jj = 1, jpj  
     375         DO ji = 1, jpi  
     376            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj))  
     377         END DO  
     378      END DO  
     379 
     380      DO jj = 1, jpj  
     381         DO ji = 1, jpi  
     382!CDIR NOVECTOR  
     383            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     384               IF ( ll_found(ji,jj) ) EXIT  
     385               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     386                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     387                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     388                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     389                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     390                  EXIT                                                    
     391               END IF  
     392            END DO  
     393         END DO  
     394      END DO  
     395 
     396      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     397      !   
     398      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     399      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     400      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT )  
     401      !  
     402   END SUBROUTINE zdf_mxl_zint_mld 
     403 
     404   SUBROUTINE zdf_mxl_zint_htc( kt ) 
     405      !!---------------------------------------------------------------------- 
     406      !!                  ***  ROUTINE zdf_mxl_zint_htc  *** 
     407      !!  
     408      !! ** Purpose :    
     409      !! 
     410      !! ** Method  :    
     411      !!---------------------------------------------------------------------- 
     412 
     413      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     414 
     415      INTEGER :: ji, jj, jk 
     416      INTEGER :: ikmax 
     417      REAL(wp) :: zc, zcoef 
     418      ! 
     419      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel 
     420      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick 
     421 
     422      !!---------------------------------------------------------------------- 
     423 
     424      IF( .NOT. ALLOCATED(ilevel) ) THEN 
     425         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 
     426         &         zthick(jpi,jpj), STAT=ji ) 
     427         IF( lk_mpp  )   CALL mpp_sum(ji) 
     428         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 
     429      ENDIF 
     430 
     431      ! Find last whole model T level above the MLD 
     432      ilevel(:,:)   = 0 
     433      zthick_0(:,:) = 0._wp 
     434 
     435      DO jk = 1, jpkm1   
     436         DO jj = 1, jpj 
     437            DO ji = 1, jpi                     
     438               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk) 
     439               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk 
     440            END DO 
     441         END DO 
     442         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
     443         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1) 
     444      END DO 
     445 
     446      ! Surface boundary condition 
     447      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                    
     448      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)    
     449      ENDIF 
     450 
     451      ! Deepest whole T level above the MLD 
     452      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
     453 
     454      ! Integration down to last whole model T level 
     455      DO jk = 1, ikmax 
     456         DO jj = 1, jpj 
     457            DO ji = 1, jpi 
     458               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
     459               zthick(ji,jj) = zthick(ji,jj) + zc 
     460               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     461            END DO 
     462         END DO 
     463      END DO 
     464 
     465      ! Subsequent partial T level 
     466      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD 
     467 
     468      DO jj = 1, jpj 
     469         DO ji = 1, jpi 
     470            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
     471      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
     472         END DO 
     473      END DO 
     474 
     475      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 
     476 
     477      ! Convert to heat content 
     478      zcoef = rau0 * rcp 
     479      htc_mld(:,:) = zcoef * htc_mld(:,:) 
     480 
     481   END SUBROUTINE zdf_mxl_zint_htc 
     482 
     483   SUBROUTINE zdf_mxl_zint( kt ) 
     484      !!---------------------------------------------------------------------- 
     485      !!                  ***  ROUTINE zdf_mxl_zint  *** 
     486      !!  
     487      !! ** Purpose :    
     488      !! 
     489      !! ** Method  :    
     490      !!---------------------------------------------------------------------- 
     491 
     492      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     493 
     494      INTEGER :: ios 
     495      INTEGER :: jn 
     496 
     497      INTEGER :: nn_mld_diag = 0    ! number of diagnostics 
     498 
     499      INTEGER :: i_steps          ! no of timesteps per hour 
     500      INTEGER :: ierror             ! logical error message    
     501 
     502      REAL(wp) :: zdt             ! timestep variable   
     503 
     504      CHARACTER(len=1) :: cmld 
     505 
     506      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     507      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags 
     508 
     509      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     510 
     511      !!---------------------------------------------------------------------- 
     512       
     513      IF( kt == nit000 ) THEN 
     514         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     515         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     516901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 
     517 
     518         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     519         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     520902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 
     521         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     522 
     523         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 
     524 
     525         mld_diags(1) = sn_mld1 
     526         mld_diags(2) = sn_mld2 
     527         mld_diags(3) = sn_mld3 
     528         mld_diags(4) = sn_mld4 
     529         mld_diags(5) = sn_mld5 
     530 
     531         IF( nn_mld_diag > 0 ) THEN 
     532            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
     533            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     534            DO jn = 1, nn_mld_diag 
     535               WRITE(numout,*) 'MLD criterion',jn,':' 
     536               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
     537               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
     538               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
     539               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     540            END DO 
     541            WRITE(numout,*) '====================================================================' 
     542         ENDIF 
     543      ENDIF 
     544 
     545      IF( nn_mld_diag > 0 ) THEN 
     546         DO jn = 1, nn_mld_diag 
     547            WRITE(cmld,'(I1)') jn 
     548            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 
     549               CALL zdf_mxl_zint_mld( mld_diags(jn) ) 
     550 
     551               IF( iom_use( "mldzint_"//cmld ) ) THEN 
     552                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 
     553               ENDIF 
     554 
     555               IF( iom_use( "mldhtc_"//cmld ) )  THEN 
     556                  CALL zdf_mxl_zint_htc( kt ) 
     557                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   ) 
     558               ENDIF 
     559 
     560               IF( iom_use( "mldzint25h_"//cmld ) ) THEN 
     561                  IF( .NOT. mld_25h_write ) mld_25h_write = .TRUE. 
     562                  zdt = rdt 
     563                  IF( nacc == 1 ) zdt = rdtmin 
     564                  IF( MOD( 3600,INT(zdt) ) == 0 ) THEN 
     565                     i_steps = 3600/INT(zdt) 
     566                  ELSE 
     567                     CALL ctl_stop('STOP', 'zdf_mxl_zint 25h: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible') 
     568                  ENDIF 
     569                  IF( ( mld_25h_init ) .OR. ( kt == nit000 ) ) THEN 
     570                     i_cnt_25h = 1 
     571                     IF( .NOT. ALLOCATED(hmld_zint_25h) ) THEN 
     572                        ALLOCATE( hmld_zint_25h(jpi,jpj,nn_mld_diag), STAT=ierror ) 
     573                        IF( ierror > 0 )  CALL ctl_stop( 'zdf_mxl_zint 25h: unable to allocate hmld_zint_25h' )    
     574                     ENDIF 
     575                     hmld_zint_25h(:,:,jn) = hmld_zint(:,:) 
     576                  ENDIF 
     577                  IF( MOD( kt, i_steps ) == 0 .AND.  kt .NE. nn_it000 ) THEN 
     578                     hmld_zint_25h(:,:,jn) = hmld_zint_25h(:,:,jn) + hmld_zint(:,:) 
     579                  ENDIF 
     580                  IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
     581                     CALL iom_put( "mldzint25h_"//cmld , hmld_zint_25h(:,:,jn) / 25._wp   ) 
     582                  ENDIF 
     583               ENDIF 
     584 
     585            ENDIF 
     586         END DO 
     587                   
     588         IF(  mld_25h_write  ) THEN 
     589            IF( ( MOD( kt, i_steps ) == 0 ) .OR.  mld_25h_init ) THEN 
     590               IF (lwp) THEN 
     591                  WRITE(numout,*) 'zdf_mxl_zint (25h) : Summed the following number of hourly values so far',i_cnt_25h 
     592          ENDIF 
     593               i_cnt_25h = i_cnt_25h + 1 
     594               IF( mld_25h_init ) mld_25h_init = .FALSE. 
     595            ENDIF 
     596            IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
     597               i_cnt_25h = 1  
     598               DO jn = 1, nn_mld_diag  
     599                     hmld_zint_25h(:,:,jn) = hmld_zint(:,:)  
     600               ENDDO  
     601            ENDIF 
     602         ENDIF 
     603                   
     604      ENDIF 
     605 
     606   END SUBROUTINE zdf_mxl_zint 
    147607 
    148608   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.