MODULE trdmld !!====================================================================== !! *** MODULE trdmld *** !! Ocean diagnostics: mixed layer T-S trends !!===================================================================== #if defined key_trdmld || defined key_esopa !!---------------------------------------------------------------------- !! 'key_trdmld' mixed layer trend diagnostics !!---------------------------------------------------------------------- !! trd_mld : T and S trends averaged over the mixed layer !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain variables USE ldftra_oce ! ocean active tracers: lateral physics USE trdtra_oce ! ocean active tracer trend variables USE zdf_oce ! ocean vertical physics USE in_out_manager ! I/O manager USE phycst ! Define parameters for the routines USE daymod ! calendar USE dianam ! build the name of file (routine) USE ldfslp ! iso-neutral slopes USE zdfmxl USE zdfddm USE ioipsl USE lbclnk #if defined key_dimgout USE diawri, ONLY : dia_wri_dimg #endif IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC trd_mld ! routine called by step.F90 !! * Shared module variables LOGICAL, PUBLIC, PARAMETER :: lk_trdmld = .TRUE. ! momentum trend flag !! * Module variables INTEGER, DIMENSION(jpi,jpj) :: & nmld, & ! mixed layer depth nbol INTEGER :: & nh_t, nmoymltrd, & ! ??? nidtrd,nhoridtrd, & ndextrd1(jpi*jpj),ndimtrd1 REAL(wp), DIMENSION(jpi,jpj) :: & rmld , & ! mld depth (m) corresponding to nmld tml , sml , & ! average T and S over mixed layer tmlb , smlb , & ! before tml and sml (kt-1) tmlbb , smlbb, & ! tml and sml at begining of the nwrite-1 ! ! timestep averaging period tmlbn , smlbn, & ! after tml and sml at time step after the ! ! begining of the NWRITE-1 timesteps tmltrdm, smltrdm ! ! REAL(wp), DIMENSION(jpi,jpj,jpltrd) :: & ! Must be jpk for mpp lbc_lnk ! TO BE FIXED ??? REAL(wp), DIMENSION(jpi,jpj,jpk) :: & tmltrd , & ! total cumulative trends of temperature and smltrd ! salinity over nwrite-1 time steps INTEGER :: iyear,imon,iday CHARACTER(LEN=80) :: clname, cltext, clmode !! * Substitutions # include "domzgr_substitute.h90" # include "ldftra_substitute.h90" # include "zdfddm_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2003) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trd_mld( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE trd_mld *** !! !! ** Purpose : computation of vertically integrated T and S budgets !! from ocean surface down to control surface (NetCDF output) !! !! ** Method/usage : !! integration done over nwrite-1 time steps !! Control surface can be either a mixed layer depth (time varying) !! or a fixed surface (jk level or bowl). !! Choose control surface with nctls in namelist NAMDIA. !! nctls = 0 : use mixed layer with density criterion !! nctls = 1 : read index from file 'ctlsurf_idx' !! nctls > 1 : use fixed level surface jk = nctls !! Note: in the remainder of the routine, the volume between the !! surface and the control surface is called "mixed-layer" !! Method check : if the control surface is fixed, the residual dh/dt !! entrainment should be zero !! !! ** Action : !! /commld/ : rmld mld depth corresponding to nmld !! tml average T over mixed layer !! tmlb tml at kt-1 !! tmlbb tml at begining of the NWRITE-1 !! time steps averaging period !! tmlbn tml at time step after the !! begining of the NWRITE-1 time !! steps averaging period !! !! mixed layer trends : !! !! tmltrd (,,1) = zonal advection !! tmltrd (,,2) = meridional advection !! tmltrd (,,3) = vertical advection !! tmltrd (,,4) = lateral diffusion (horiz. component+Beckman) !! tmltrd (,,5) = forcing !! tmltrd (,,6) = entrainment due to vertical diffusion (TKE) !! if iso tmltrd (,,7) = lateral diffusion (vertical component) !! tmltrd (,,8) = eddy induced zonal advection !! tmltrd (,,9) = eddy induced meridional advection !! tmltrd (,,10) = eddy induced vertical advection !! !! tmltrdm(,) : total cumulative trends over nwrite-1 time steps !! ztmltot(,) : dT/dt over the NWRITE-1 time steps !! averaging period (including Asselin !! terms) !! ztmlres(,) : residual = dh/dt entrainment !! !! trends output in netCDF format using ioipsl !! !! History : !! ! 95-04 (J. Vialard) Original code !! ! 97-02 (E. Guilyardi) Adaptation global + base cmo !! ! 99-09 (E. Guilyardi) Re-writing + netCDF output !! 8.5 ! 02-06 (G. Madec) F90: Free form and module !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ilseq INTEGER :: ji, jj, jk, jl, ik, ikb, idebug, isum, it INTEGER, DIMENSION(jpi,jpj) :: zvlmsk REAL(wp) :: zmean, zavt, zjulian, zsto, zout REAL(wp) ,DIMENSION(jpi,jpj,jpktrd) :: zwkx REAL(wp) ,DIMENSION(jpi,jpj) :: & & ztmltot, ztmlres, z2d, zsmltot, zsmlres CHARACTER (len=21) :: & clold ='OLD' , & ! open specifier (direct access files) clunf ='UNFORMATTED', & ! open specifier (direct access files) clseq ='SEQUENTIAL' ! open specifier (direct access files) CHARACTER (len=80) :: clname CHARACTER (len=40) :: clhstnam CHARACTER (len=40) :: clop CHARACTER (len=12) :: clmxl NAMELIST/namtrd/ ntrd, nctls !!---------------------------------------------------------------------- ! =================== ! 0. initialization ! =================== ! Open specifier ilseq = 1 idebug = 0 ! set it to 1 in case of problem to have more print IF( kt == nit000 ) THEN ! namelist namtrd : trend diagnostic REWIND( numnam ) READ ( numnam, namtrd ) IF(lwp) THEN WRITE(numout,*) 'namtrd' WRITE(numout,*) ' ' WRITE(numout,*) ' time step frequency trend ntrd = ',ntrd WRITE(numout,*) ' control surface for trends nctls = ',nctls WRITE(numout,*) ' ' ENDIF ! cumulated trends array init nmoymltrd = 0 tmltrdm(:,:) = 0. smltrdm(:,:) = 0. ENDIF ! set before values of vertically average T and S IF( kt > nit000 ) THEN tmlb(:,:) = tml(:,:) smlb(:,:) = sml(:,:) ENDIF ! read control surface from file ctlsurf_idx IF( kt == nit000 .and. nctls == - 1 ) THEN clname ='ctlsurf_idx' CALL ctlopn(numbol,clname,clold,clunf,clseq, & ilseq,numout,lwp,1) REWIND (numbol) READ(numbol) nbol ENDIF IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: 0. done ' CALL FLUSH(numout) ENDIF ! ======================================================== ! I. definition of control surface and associated fields ! ======================================================== ! I.1 set nmld(ji,jj) = index of first T point below control surface ! ------------------- or outside mixed-layer ! clmxl = legend root for netCDF output IF( nctls == 0 ) THEN ! control surface = mixed-layer with density criterion ! (array nmln computed in zdfmxl.F90) nmld(:,:) = nmln(:,:) clmxl = 'Mixed Layer ' ELSE IF( nctls == 1 ) THEN ! control surface = read index from file nmld(:,:) = nbol(:,:) clmxl = ' Bowl ' ELSE IF( nctls >= 2 ) THEN ! control surface = model level nctls = MIN( nctls, jpktrd - 1 ) nmld(:,:) = nctls + 1 WRITE(clmxl,'(A9,I2,1X)') 'Levels 1-', nctls ENDIF ! Check of validity : nmld(ji,jj) =< jpktrd isum = 0 IF( jpktrd < jpk ) THEN DO jj = 1, jpj DO ji = 1, jpi IF( nmld(ji,jj) <= jpktrd ) THEN zvlmsk(ji,jj) = tmask(ji,jj,1) ELSE isum = isum + 1 zvlmsk(ji,jj) = 0. ENDIF END DO END DO ENDIF IF( idebug /= 0 ) THEN ! CALL prihre (zvlmsk,jpi,jpj,1,jpi,2,1,jpj,2,3,numout) WRITE(numout,*) ' debuging trd_mld: I.1 done ' CALL FLUSH(numout) ENDIF ! I.2 probability density function of presence in mixed-layer ! -------------------------------- ! (i.e. weight of each grid point in vertical integration : zwkx(ji,jj,jk) ! initialize zwkx with vertical scale factor in mixed-layer zwkx(:,:,:) = 0.e0 DO jk = 1, jpktrd DO jj = 1,jpj DO ji = 1,jpi IF( jk - nmld(ji,jj) < 0. ) zwkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk) END DO END DO END DO ! compute mixed-layer depth : rmld rmld(:,:) = 0. DO jk = 1, jpktrd rmld(:,:) = rmld(:,:) + zwkx(:,:,jk) END DO ! compute PDF DO jk = 1, jpktrd zwkx(:,:,jk) = zwkx(:,:,jk) / MAX( 1., rmld(:,:) ) END DO IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: I.2 done ' CALL FLUSH(numout) ENDIF ! I.3 vertically integrated T and S ! --------------------------------- tml(:,:) = 0. sml(:,:) = 0. DO jk = 1, jpktrd - 1 tml(:,:) = tml(:,:) + zwkx(:,:,jk) * tn(:,:,jk) sml(:,:) = sml(:,:) + zwkx(:,:,jk) * sn(:,:,jk) END DO IF(idebug /= 0) THEN WRITE(numout,*) ' debuging trd_mld: I.3 done' CALL FLUSH(numout) ENDIF ! =================================== ! II. netCDF output initialization ! =================================== #if defined key_dimgout #else #include "trdmld_ncinit.h90" #endif IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: II. done' CALL FLUSH(numout) ENDIF ! ==================================================== ! III. vertical integration of trends in mixed-layer ! ==================================================== ! III.0 initializations ! --------------------- tmltrd(:,:,:) = 0.e0 smltrd(:,:,:) = 0.e0 IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: III.0 done' CALL FLUSH(numout) ENDIF ! III.1 vertical integration of 3D trends ! --------------------------------------- DO jk = 1,jpktrd ! Temperature tmltrd(:,:,1) = tmltrd(:,:,1) + ttrdh(:,:,jk,1) * zwkx(:,:,jk) ! zonal advection tmltrd(:,:,2) = tmltrd(:,:,2) + ttrdh(:,:,jk,2) * zwkx(:,:,jk) ! meridional advection tmltrd(:,:,3) = tmltrd(:,:,3) + ttrd (:,:,jk,2) * zwkx(:,:,jk) ! vertical advection tmltrd(:,:,4) = tmltrd(:,:,4) + ttrd (:,:,jk,3) * zwkx(:,:,jk) ! lateral diffusion (hor. part) tmltrd(:,:,5) = tmltrd(:,:,5) + ttrd (:,:,jk,7) * zwkx(:,:,jk) ! forcing (penetrative) IF( l_traldf_iso ) THEN tmltrd(:,:,7) = tmltrd(:,:,7) + ttrd (:,:,jk,4) * zwkx(:,:,jk) ! lateral diffusion (explicit ! ! vert. part (isopycnal diff.) ENDIF !#if defined key_traldf_eiv ! tmltrd(:,:,8 ) = tmltrdg(:,:,8) + ttrdh(:,:,jk,3) * zwkx(:,:,jk) ! eddy induced zonal advection ! tmltrd(:,:,9 ) = tmltrdg(:,:,9) + ttrdh(:,:,jk,4) * zwkx(:,:,jk) ! eddy induced merid. advection ! tmltrd(:,:,10) = tmltrdg(:,:,10) + ttrd(:,:,jk,6) * zwkx(:,:,jk) ! eddy induced vert. advection !#endif ! Salinity smltrd(:,:,1) = smltrd(:,:,1) + strdh(:,:,jk,1) * zwkx(:,:,jk) ! zonal advection smltrd(:,:,2) = smltrd(:,:,2) + strdh(:,:,jk,2) * zwkx(:,:,jk) ! meridional advection smltrd(:,:,3) = smltrd(:,:,3) + strd (:,:,jk,2) * zwkx(:,:,jk) ! vertical advection smltrd(:,:,4) = smltrd(:,:,4) + strd (:,:,jk,3) * zwkx(:,:,jk) ! lateral diffusion (hor. part) IF( l_traldf_iso ) THEN smltrd(:,:,7) = smltrd(:,:,7) + strd (:,:,jk,4) * zwkx(:,:,jk) ! lateral diffusion (explicit ! ! vert. part (isopycnal diff.) ENDIF !#if defined key_traldf_eiv ! smltrd(:,:,8) = smltrdg(:,:,8) + strdh(:,:,jk,3) * zwkx(:,:,jk) ! eddy induced zonal advection ! smltrd(:,:,9) = smltrdg(:,:,9) + strdh(:,:,jk,4) * zwkx(:,:,jk) ! eddy induced merid. advection ! smltrd(:,:,10) = smltrdg(:,:,10) + strd(:,:,jk,6) * zwkx(:,:,jk) ! eddy induced vert. advection !#endif END DO IF( idebug /= 0 ) THEN IF(lwp) WRITE(numout,*) ' debuging trd_mld: III.1 done' CALL FLUSH(numout) ENDIF ! III.2 trends terms at upper and lower boundaries of mixed-layer ! --------------------------------------------------------------- DO jj = 1,jpj DO ji = 1,jpi ik = nmld(ji,jj) ! Temperature ! forcing (non penetrative) tmltrd(ji,jj,5) = tmltrd(ji,jj,5) + flxtrd(ji,jj,1) * zwkx(ji,jj,1) ! entrainment due to vertical diffusion ! - due to vertical mixing scheme (TKE) zavt = avt(ji,jj,ik) tmltrd(ji,jj,6) = - 1. * zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) ) & / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1) ! Salinity ! forcing smltrd(ji,jj,5) = flxtrd(ji,jj,2) * zwkx(ji,jj,1) ! entrainment due to vertical diffusion ! - due to vertical mixing scheme (TKE) zavt = fsavs(ji,jj,ik) smltrd(ji,jj,6) = -1. * zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & & * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) ) & & / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1) END DO END DO IF( l_traldf_iso ) THEN !!Clem On retire de la diffusion verticale TOTALE calculee par tmltrd(:,:,7) !!Clem ds trazdf.isopycnal et implicit, la partie verticale due au Kz afin de ne garder !!Clem effectivement que la diffusion verticale isopycnale (ie composante de la !!Clem diff isopycnale sur la verticale) : tmltrd(:,:,7) = tmltrd(:,:,7) - tmltrd(:,:,6) ! - due to isopycnal mixing scheme (implicit part) smltrd(:,:,7) = smltrd(:,:,7) - smltrd(:,:,6) ! - due to isopycnal mixing scheme (implicit part) ENDIF ! Boundary conditions DO jk = 1, jpltrd CALL lbc_lnk( tmltrd, 'T', 1. ) CALL lbc_lnk( smltrd, 'T', 1. ) END DO IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: III.2 done' CALL FLUSH(numout) ENDIF #if defined key_trabbl_dif ! III.3 trends terms from beckman over-flow parameterization ! ---------------------------------------------------------- DO jj = 1,jpj DO ji = 1,jpi ikb = MAX( mbathy(ji,jj)-1, 1 ) ! beckmann component -> horiz. part of lateral diffusion tmltrd(ji,jj,4) = tmltrd(ji,jj,4) + bbltrd(ji,jj,1) * zwkx(ji,jj,ikb) smltrd(ji,jj,4) = smltrd(ji,jj,4) + bbltrd(ji,jj,2) * zwkx(ji,jj,ikb) END DO END DO #endif IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: III.3 done' CALL FLUSH(numout) ENDIF ! ================================= ! IV. Cumulated trends ! ================================= ! IV.1 set `before' mixed layer values for kt = nit000+1 ! -------------------------------------------------------- IF( kt == nit000+1 ) THEN tmlbb(:,:) = tmlb(:,:) tmlbn(:,:) = tml (:,:) smlbb(:,:) = smlb(:,:) smlbn(:,:) = sml (:,:) ENDIF IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: IV.1 done' CALL FLUSH(numout) ENDIF ! IV.2 cumulated trends over analysis period (kt=2 to nwrite) ! ---------------------- ! trends cumulated over nwrite-2 time steps IF( kt >= nit000+2 ) THEN nmoymltrd = nmoymltrd + 1 DO jl = 1, jpltrd tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl) smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl) END DO ENDIF IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: IV.2 done' CALL FLUSH(numout) ENDIF ! ============================================= ! V. Output in netCDF + residual computation ! ============================================= IF( MOD( kt - nit000+1, nwrite ) == 0 ) THEN ! V.1 compute total trend ! ------------------------ zmean = float(nmoymltrd) ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / (zmean * 2. * rdt) zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / (zmean * 2. * rdt) IF(idebug /= 0) THEN WRITE(numout,*) ' zmean = ',zmean WRITE(numout,*) ' debuging trd_mld: V.1 done' CALL FLUSH(numout) ENDIF ! V.2 compute residual ! --------------------- ztmlres(:,:) = ztmltot(:,:) - tmltrdm(:,:) / zmean zsmlres(:,:) = zsmltot(:,:) - smltrdm(:,:) / zmean ! Boundary conditions CALL lbc_lnk( ztmltot, 'T', 1. ) CALL lbc_lnk( ztmlres, 'T', 1. ) CALL lbc_lnk( zsmltot, 'T', 1. ) CALL lbc_lnk( zsmlres, 'T', 1. ) IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: V.2 done' CALL FLUSH(numout) ENDIF ! V.3 time evolution array swap ! ------------------------------ tmlbb(:,:) = tmlb(:,:) tmlbn(:,:) = tml (:,:) smlbb(:,:) = smlb(:,:) smlbn(:,:) = sml (:,:) IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: V.3 done' CALL FLUSH(numout) ENDIF ! V.4 zero cumulative array ! --------------------------- nmoymltrd = 0 tmltrdm(:,:) = 0. smltrdm(:,:) = 0. IF(idebug /= 0) THEN WRITE(numout,*) ' debuging trd_mld: IV.4 done' CALL FLUSH(numout) ENDIF ENDIF ! IV.5 write trends to output ! --------------------------- #if defined key_dimgout ! code for dimg mpp output IF ( MOD(kt,nwrite) == 0 ) THEN WRITE(clmode,'(f5.1,a)' ) nwrite*rdt/86400.,' days average' iyear = ndastp/10000 imon = (ndastp-iyear*10000)/100 iday = ndastp - imon*100 - iyear*10000 WRITE(clname,9000) TRIM(cexper),'MLDiags',iyear,imon,iday cltext=TRIM(cexper)//' mld diags'//TRIM(clmode) CALL dia_wri_dimg (clname, cltext, smltrd, jpltrd, '2') 9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc") END IF #else IF( kt >= nit000+1 ) THEN #include "trdmld_ncwrite.h90" IF( idebug /= 0 ) THEN WRITE(numout,*) ' debuging trd_mld: IV.5 done' CALL FLUSH(numout) ENDIF ENDIF IF( kt == nitend ) CALL histclo( nidtrd ) #endif END SUBROUTINE trd_mld #else !!---------------------------------------------------------------------- !! Default option : Empty module !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_trdmld = .FALSE. ! momentum trend flag CONTAINS SUBROUTINE trd_mld( kt ) ! Empty routine WRITE(*,*) kt END SUBROUTINE trd_mld #endif !!====================================================================== END MODULE trdmld