MODULE trdmld_trc !!====================================================================== !! *** MODULE trdmld_trc *** !! Ocean diagnostics: mixed layer passive tracer trends !!====================================================================== !! History : 9.0 ! 06-08 (C. Deltel) Original code (from trdmld.F90) !! ! 07-04 (C. Deltel) Bug fix : add trcrad trends !! ! 07-06 (C. Deltel) key_gyre : do not call lbc_lnk !!---------------------------------------------------------------------- #if defined key_top && ( defined key_trdmld_trc || defined key_esopa ) !!---------------------------------------------------------------------- !! 'key_trdmld_trc' mixed layer trend diagnostics !!---------------------------------------------------------------------- !! trd_mld_trc : passive tracer cumulated trends averaged over ML !! trd_mld_trc_zint : passive tracer trends vertical integration !! trd_mld_trc_init : initialization step !!---------------------------------------------------------------------- USE trp_trc ! tracer definitions (trn, trb, tra, etc.) USE oce_trc ! needed for namelist logicals, and euphotic layer arrays USE trctrp_lec USE trdmld_trc_oce ! definition of main arrays used for trends computations USE in_out_manager ! I/O manager USE dianam ! build the name of file (routine) USE ldfslp ! iso-neutral slopes USE ioipsl ! NetCDF library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE trdmld_trc_rst ! restart for diagnosing the ML trends USE prtctl ! print control USE sms_pisces USE sms_lobster USE trcsms_cfc USE trc USE trcrst ! for lrst_trc -> circ. dep. ??? we put lrst_trc in trc_oce IMPLICIT NONE PRIVATE PUBLIC trd_mod_trc ! routine called by step.F90 PUBLIC trd_mld_trc PUBLIC trd_mld_trc_init CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file INTEGER :: nmoymltrd INTEGER :: ndextrd1(jpi*jpj) INTEGER, DIMENSION(jptra) :: nidtrd, nh_t INTEGER :: ndimtrd1 INTEGER, SAVE :: ionce, icount LOGICAL :: llwarn = .TRUE. ! this should always be .TRUE. LOGICAL :: lldebug = .TRUE. !! * Substitutions # include "top_substitute.h90" !!---------------------------------------------------------------------- !! TOP 1.0 , LOCEAN-IPSL (2007) !! $Header: $ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trd_mod_trc( ptrtrd, kjn, ktrd, kt ) !!---------------------------------------------------------------------- !! *** ROUTINE trd_mod_trc *** !!---------------------------------------------------------------------- #if defined key_trcbbl_adv REAL(wp), DIMENSION(jpi,jpj,jpk) :: zun, zvn ! temporary arrays #else USE oce_trc, zun => un ! When no bbl, zun == un USE oce_trc, zvn => vn ! When no bbl, zvn == vn #endif INTEGER, INTENT( in ) :: kt ! time step INTEGER, INTENT( in ) :: kjn ! tracer index INTEGER, INTENT( in ) :: ktrd ! tracer trend index REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: ptrtrd ! Temperature or U trend !!---------------------------------------------------------------------- IF( kt == nittrc000 ) THEN ! IF(lwp)WRITE(numout,*) ! IF(lwp)WRITE(numout,*) 'trd_mod_trc:' ! IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~' ENDIF !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Mixed layer trends for passive tracers !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< SELECT CASE ( ktrd ) CASE ( jptrc_trd_xad ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xad , '3D', kjn ) CASE ( jptrc_trd_yad ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yad , '3D', kjn ) CASE ( jptrc_trd_zad ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zad , '3D', kjn ) CASE ( jptrc_trd_ldf ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf , '3D', kjn ) CASE ( jptrc_trd_xei ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xei , '3D', kjn ) CASE ( jptrc_trd_yei ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yei , '3D', kjn ) CASE ( jptrc_trd_bbl ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbl , '3D', kjn ) CASE ( jptrc_trd_zdf ) IF( ln_trcldf_iso ) THEN CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf, '3D', kjn ) ELSE CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zdf, '3D', kjn ) ENDIF CASE ( jptrc_trd_zei ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zei , '3D', kjn ) CASE ( jptrc_trd_dmp ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_dmp , '3D', kjn ) CASE ( jptrc_trd_sbc ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sbc , '2D', kjn ) #if defined key_lobster CASE ( jptrc_trd_sms_sed ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_sed, '3D', kjn ) CASE ( jptrc_trd_sms_bio ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_bio, '3D', kjn ) CASE ( jptrc_trd_sms_exp ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_exp, '3D', kjn ) #else CASE ( jptrc_trd_sms ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms , '3D', kjn ) #endif CASE ( jptrc_trd_bbc ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbc , '3D', kjn ) CASE ( jptrc_trd_radb ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radb , '3D', kjn ) CASE ( jptrc_trd_radn ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radn , '3D', kjn ) CASE ( jptrc_trd_atf ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_atf , '3D', kjn ) END SELECT END SUBROUTINE trd_mod_trc SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn ) !!---------------------------------------------------------------------- !! *** ROUTINE trd_mld_trc_zint *** !! !! ** Purpose : Compute the vertical average of the 3D fields given as arguments !! to the subroutine. This vertical average is performed from ocean !! surface down to a chosen control surface. !! !! ** Method/usage : !! The control surface can be either a mixed layer depth (time varying) !! or a fixed surface (jk level or bowl). !! Choose control surface with nctls_trc in namelist NAMTRD : !! nctls_trc = -2 : use isopycnal surface !! nctls_trc = -1 : use euphotic layer with light criterion !! nctls_trc = 0 : use mixed layer with density criterion !! nctls_trc = 1 : read index from file 'ctlsurf_idx' !! nctls_trc > 1 : use fixed level surface jk = nctls_trc !! Note: in the remainder of the routine, the volume between the !! surface and the control surface is called "mixed-layer" !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: ptrc_trdmld ! passive tracer trend INTEGER :: ji, jj, jk, isum REAL(wp), DIMENSION(jpi,jpj) :: zvlmsk !!---------------------------------------------------------------------- ! I. Definition of control surface and integration weights ! -------------------------------------------------------- ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN ! tmltrd_trc(:,:,:,:) = 0.e0 ! <<< reset trend arrays to zero ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer SELECT CASE ( nctls_trc ) ! choice of the control surface CASE ( -2 ) ; STOP 'trdmld_trc : not ready ' ! -> isopycnal surface (see ???) #if defined key_pisces || defined key_lobster CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion #endif CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl) CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file CASE ( 2: ) ; nctls_trc = MIN( nctls_trc, jpktrd_trc - 1 ) nmld_trc(:,:) = nctls_trc + 1 ! -> model level END SELECT ! ... Compute ndextrd1 and ndimtrd1 ??? role de jpktrd_trc IF( ionce == 1 ) THEN ! isum = 0 ; zvlmsk(:,:) = 0.e0 IF( jpktrd_trc < jpk ) THEN ! description ??? DO jj = 1, jpj DO ji = 1, jpi IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN zvlmsk(ji,jj) = tmask(ji,jj,1) ELSE isum = isum + 1 zvlmsk(ji,jj) = 0.e0 ENDIF END DO END DO ENDIF IF( isum > 0 ) THEN ! index of ocean points (2D only) WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 ) ELSE CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 ) ENDIF ionce = 0 ! no more pass here ! ENDIF ! ionce == 1 ! ... Weights for vertical averaging wkx_trc(:,:,:) = 0.e0 DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer DO jj = 1, jpj DO ji = 1, jpi IF( jk - nmld_trc(ji,jj) < 0 ) wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk) END DO END DO END DO rmld_trc(:,:) = 0.e0 DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk) END DO DO jk = 1, jpktrd_trc ! compute integration weights wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) ) END DO icount = 0 ! <<< flag = off : control surface & integr. weights ! ! computed only once per time step ENDIF ONCE_PER_TIME_STEP ! II. Vertical integration of trends in the mixed-layer ! ----------------------------------------------------- SELECT CASE ( ctype ) CASE ( '3D' ) ! mean passive tracer trends in the mixed-layer DO jk = 1, jpktrd_trc tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,jk) * wkx_trc(:,:,jk) END DO CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,1) * wkx_trc(:,:,1) ! non penetrative END SELECT END SUBROUTINE trd_mld_trc_zint SUBROUTINE trd_mld_trc( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE trd_mld_trc *** !! !! ** Purpose : Compute and cumulate the mixed layer trends over an analysis !! period, and write NetCDF (or dimg) outputs. !! !! ** Method/usage : !! The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant !! logical namelist variable) : !! 1) to explain the difference between initial and final !! mixed-layer T & S (where initial and final relate to the !! current analysis window, defined by ntrc_trc in the namelist) !! 2) to explain the difference between the current and previous !! TIME-AVERAGED mixed-layer T & S (where time-averaging is !! performed over each analysis window). !! !! ** Consistency check : !! If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt !! entrainment) should be zero, at machine accuracy. Note that in the case !! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO !! over the first two analysis windows (except if restart). !! N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, ucf_trc=1., nctls_trc=8 !! for checking residuals. !! On a NEC-SX5 computer, this typically leads to: !! O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false. !! O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true. !! !! ** Action : !! At each time step, mixed-layer averaged trends are stored in the !! tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx). !! This array is known when trd_mld is called, at the end of the stp subroutine, !! except for the purely vertical K_z diffusion term, which is embedded in the !! lateral diffusion trend. !! !! In I), this K_z term is diagnosed and stored, thus its contribution is removed !! from the lateral diffusion trend. !! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative !! arrays are updated. !! In III), called only once per analysis window, we compute the total trends, !! along with the residuals and the Asselin correction terms. !! In IV), the appropriate trends are written in the trends NetCDF file. !! !! References : !! - Vialard & al. !! - See NEMO documentation (in preparation) !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! ocean time-step index INTEGER :: ji, jj, jk, jl, ik, it, jn REAL(wp) :: zavt, zfn, zfn2 !! REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmltot ! d(trc)/dt over the anlysis window (incl. Asselin) REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlres ! residual = dh/dt entrainment term REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlatf ! for storage only REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlrad ! for storage only (for trb<0 corr in trcrad) !! REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmltot2 ! -+ REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlres2 ! | working arrays to diagnose the trends REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmltrdm2 ! | associated with the time meaned ML REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlatf2 ! | passive tracers REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlrad2 ! | (-> for trb<0 corr in trcrad) REAL(wp), DIMENSION(jpi,jpj,jpltrd_trc,jptra) :: ztmltrd2 ! -+ !! REAL(wp), DIMENSION(jpi,jpj) :: z2d ! temporary array, used for eiv arrays CHARACTER (LEN= 5) :: clvar #if defined key_dimgout INTEGER :: iyear,imon,iday CHARACTER(LEN=80) :: cltext, clmode #endif !!---------------------------------------------------------------------- IF( llwarn ) THEN ! warnings IF( ( nittrc000 /= nit000 ) & .OR.( ndttrc /= 1 ) ) THEN WRITE(numout,*) 'Be careful, trends diags never validated' STOP 'Uncomment this line to proceed' ENDIF ENDIF ! ====================================================================== ! I. Diagnose the purely vertical (K_z) diffusion trend ! ====================================================================== ! ... These terms can be estimated by flux computation at the lower boundary of the ML ! (we compute (-1/h) * K_z * d_z( tracer )) IF( ln_trcldf_iso ) THEN ! DO jj = 1,jpj DO ji = 1,jpi ik = nmld_trc(ji,jj) zavt = avt(ji,jj,ik) DO jn = 1, jptra IF( luttrd(jn) ) & tmltrd_trc(ji,jj,jpmld_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & & * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) ) & & / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1) END DO END DO END DO DO jn = 1, jptra ! ... Remove this K_z trend from the iso-neutral diffusion term (if any) IF( luttrd(jn) ) & tmltrd_trc(:,:,jpmld_trc_ldf,jn) = tmltrd_trc(:,:,jpmld_trc_ldf,jn) - tmltrd_trc(:,:,jpmld_trc_zdf,jn) END DO ! ENDIF #if ! defined key_gyre ! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm. ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.) DO jn = 1, jptra IF( luttrd(jn) ) THEN DO jl = 1, jpltrd_trc CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions END DO ENDIF END DO #endif ! ====================================================================== ! II. Cumulate the trends over the analysis window ! ====================================================================== ztmltrd2(:,:,:,:) = 0.e0 ; ztmltot2(:,:,:) = 0.e0 ! <<< reset arrays to zero ztmlres2(:,:,:) = 0.e0 ; ztmlatf2(:,:,:) = 0.e0 ztmlrad2(:,:,:) = 0.e0 ! II.1 Set before values of vertically averages passive tracers ! ------------------------------------------------------------- IF( kt > nittrc000 ) THEN DO jn = 1, jptra IF( luttrd(jn) ) THEN tmlb_trc (:,:,jn) = tml_trc (:,:,jn) tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_atf,jn) tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_radb,jn) ENDIF END DO ENDIF ! II.2 Vertically averaged passive tracers ! ---------------------------------------- tml_trc(:,:,:) = 0.e0 DO jk = 1, jpktrd_trc ! - 1 ??? DO jn = 1, jptra IF( luttrd(jn) ) & tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn) END DO END DO ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window ! ------------------------------------------------------------------------ IF( kt == 2 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) ??? ! DO jn = 1, jptra IF( luttrd(jn) ) THEN tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn) tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn) tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0 ; tmltrd_atf_sumb_trc (:,:,jn) = 0.e0 tmltrd_rad_sumb_trc (:,:,jn) = 0.e0 ENDIF END DO rmldbn_trc(:,:) = rmld_trc(:,:) ! ENDIF ! II.4 Cumulated trends over the analysis period ! ---------------------------------------------- ! ! [ 1rst analysis window ] [ 2nd analysis window ] ! ! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps ! ntrd 2*ntrd etc. ! 1 2 3 4 =5 e.g. =10 ! IF( ( kt >= 2 ).OR.( lrsttr ) ) THEN ! ??? ! nmoymltrd = nmoymltrd + 1 ! ... Cumulate over BOTH physical contributions AND over time steps DO jn = 1, jptra IF( luttrd(jn) ) THEN DO jl = 1, jpltrd_trc tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn) END DO ENDIF END DO DO jn = 1, jptra IF( luttrd(jn) ) THEN ! ... Special handling of the Asselin trend tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn) tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn) ! ... Trends associated with the time mean of the ML passive tracers tmltrd_sum_trc (:,:,:,jn) = tmltrd_sum_trc (:,:,:,jn) + tmltrd_trc (:,:,:,jn) tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn) tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) + tml_trc (:,:,jn) ENDIF ENDDO rmld_sum_trc (:,:) = rmld_sum_trc (:,:) + rmld_trc (:,:) ! ENDIF ! ====================================================================== ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD) ! ====================================================================== ! Convert to appropriate physical units tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * ucf_trc MODULO_NTRD : IF( MOD( kt, ntrd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd_trc ! ztmltot (:,:,:) = 0.e0 ! reset arrays to zero ztmlres (:,:,:) = 0.e0 ztmltot2(:,:,:) = 0.e0 ztmlres2(:,:,:) = 0.e0 zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn ! III.1 Prepare fields for output ("instantaneous" diagnostics) ! ------------------------------------------------------------- DO jn = 1, jptra IF( luttrd(jn) ) THEN !-- Compute total trends (use rdttrc instead of rdt ???) IF ( ln_trcadv_smolar .OR. ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt ELSE ! LEAP-FROG schemes ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt) ENDIF !-- Compute residuals ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) & & - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) ) !-- Diagnose Asselin trend over the analysis window ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) !-- Lateral boundary conditions #if ! defined key_gyre CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. ) CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. ) #endif #if defined key_diainstant STOP 'tmltrd_trc : key_diainstant was never checked within trdmld. Comment this to proceed.' #endif ENDIF END DO ! III.2 Prepare fields for output ("mean" diagnostics) ! ---------------------------------------------------- !-- Update the ML depth time sum (to build the Leap-Frog time mean) rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:) !-- Compute passive tracer total trends DO jn = 1, jptra IF( luttrd(jn) ) THEN tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn) ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec ENDIF END DO !-- Compute passive tracer residuals DO jn = 1, jptra IF( luttrd(jn) ) THEN ! DO jl = 1, jpltrd_trc ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn) END DO ztmltrdm2(:,:,jn) = 0.e0 DO jl = 1, jpltrd_trc ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn) END DO ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - & & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) & & - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) ) ! !-- Diagnose Asselin trend over the analysis window ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) & & + tmltrd_atf_sumb_trc(:,:,jn) ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) & & + tmltrd_rad_sumb_trc(:,:,jn) !-- Lateral boundary conditions #if ! defined key_gyre CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. ) CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. ) DO jl = 1, jpltrd_trc CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file END DO #endif ENDIF END DO ! * Debugging information * IF( lldebug ) THEN ! WRITE(numout,*) 'trd_mld_trc : write trends in the Mixed Layer for debugging process:' WRITE(numout,*) '~~~~~~~~~~~ ' WRITE(numout,*) WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd DO jn = 1, jptra IF( luttrd(jn) ) THEN WRITE(numout, *) WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<' WRITE(numout, *) WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn)) !CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ') WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn))) WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- ' WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn)) WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- ' WRITE(numout, *) WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn))) WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ ' WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_atf,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_radb,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) ) WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- ' WRITE(numout, *) WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn)) WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- ' WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn)) WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- ' WRITE(numout, *) WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn)) WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn)) WRITE(numout, *) DO jl = 1, jpltrd_trc WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmld_trc_xxx = ', jl, & & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn)) END DO WRITE(numout,*) WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** ' WRITE(numout,*) WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn) WRITE(numout,*) WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn) WRITE(numout,*) WRITE(numout,*) ' *********************** ZTMLRES *********************** ' WRITE(numout,*) WRITE(numout,*) '...................................................' WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : ' DO jj = 5, 1, -1 WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 ) END DO WRITE(numout,*) WRITE(numout,*) ' *********************** ZTMLRES2 *********************** ' WRITE(numout,*) WRITE(numout,*) '...................................................' WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : ' DO jj = 5, 1, -1 WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 ) END DO ! ENDIF ! END DO 97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10) 98 FORMAT(a10, i3, 2x, a30, 2x, g20.10) 99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x)) WRITE(numout,*) ! ENDIF ! III.3 Time evolution array swap ! ------------------------------- ! ML depth rmldbn_trc(:,:) = rmld_trc(:,:) rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum DO jn = 1, jptra IF( luttrd(jn) ) THEN ! For passive tracer instantaneous diagnostics tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn) tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn) ! For passive tracer mean diagnostics tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn) tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn) tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) ! III.4 Convert to appropriate physical units ! ------------------------------------------- ztmltot (:,:,jn) = ztmltot (:,:,jn) * ucf_trc/zfn ! instant diags ztmlres (:,:,jn) = ztmlres (:,:,jn) * ucf_trc/zfn ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * ucf_trc/zfn ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * ucf_trc/zfn tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * ucf_trc/zfn2 ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * ucf_trc/zfn2 ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * ucf_trc/zfn2 ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * ucf_trc/zfn2 ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * ucf_trc/zfn2 ENDIF END DO ! ENDIF MODULO_NTRD ! ====================================================================== ! IV. Write trends in the NetCDF file ! ====================================================================== ! IV.1 Code for dimg mpp output ! ----------------------------- # if defined key_dimgout STOP 'Not implemented' # else ! IV.2 Code for IOIPSL/NetCDF output ! ---------------------------------- IF( lwp .AND. MOD( kt , ntrd_trc ) == 0 ) THEN WRITE(numout,*) ' ' WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :' WRITE(numout,*) '~~~~~~~~~~~ ' WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd WRITE(numout,*) ' ' ENDIF it = kt - nit000 + 1 NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags ! DO jn = 1, jptra ! IF( luttrd(jn) ) THEN !-- Specific treatment for EIV trends ! WARNING : When eiv is switched on but key_diaeiv is not, we do NOT diagnose ! u_eiv, v_eiv, and w_eiv : the exact eiv advective trends thus cannot be computed, ! only their sum makes sense => mask directional contrib. to avoid confusion z2d(:,:) = tmltrd_trc(:,:,jpmld_trc_xei,jn) + tmltrd_trc(:,:,jpmld_trc_yei,jn) & & + tmltrd_trc(:,:,jpmld_trc_zei,jn) #if ( defined key_trcldf_eiv && defined key_diaeiv ) tmltrd_trc(:,:,jpmld_trc_xei,jn) = -999. tmltrd_trc(:,:,jpmld_trc_yei,jn) = -999. tmltrd_trc(:,:,jpmld_trc_zei,jn) = -999. #endif CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 ) !-- Output the fields clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc. CALL histwrite( nidtrd(jn), clvar , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) DO jl = 1, jpltrd_trc - 2 CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), & & it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 ) END DO CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1 & it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), & ! now Asselin : jpltrd_trc & it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)), & ! now total EIV : jpltrd_trc + 1 & it, z2d(:,:), ndimtrd1, ndextrd1 ) ! ENDIF END DO IF( kt == nitend ) THEN DO jn = 1, jptra IF( luttrd(jn) ) CALL histclo( nidtrd(jn) ) END DO ENDIF ELSE ! <<< write the trends for passive tracer mean diagnostics DO jn = 1, jptra ! IF( luttrd(jn) ) THEN !-- Specific treatment for EIV trends ! WARNING : see above z2d(:,:) = ztmltrd2(:,:,jpmld_trc_xei,jn) + ztmltrd2(:,:,jpmld_trc_yei,jn) & & + ztmltrd2(:,:,jpmld_trc_zei,jn) #if ( defined key_trcldf_eiv && defined key_diaeiv ) ztmltrd2(:,:,jpmld_trc_xei,jn) = -999. ztmltrd2(:,:,jpmld_trc_yei,jn) = -999. ztmltrd2(:,:,jpmld_trc_zei,jn) = -999. #endif CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) !-- Output the fields clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc. CALL histwrite( nidtrd(jn), clvar , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) DO jl = 1, jpltrd_trc - 2 CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), & & it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 ) END DO CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1 & it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), & ! now Asselin : jpltrd_trc & it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 ) CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)), & ! now total EIV : jpltrd_trc + 1 & it, z2d(:,:), ndimtrd1, ndextrd1 ) ENDIF ! END DO IF( kt == nitend ) THEN DO jn = 1, jptra IF( luttrd(jn) ) CALL histclo( nidtrd(jn) ) END DO ENDIF ! ENDIF NETCDF_OUTPUT ! Compute the control surface (for next time step) : flag = on icount = 1 # endif /* key_dimgout */ IF( MOD( kt, ntrd_trc ) == 0 ) THEN ! ! Reset cumulative arrays to zero ! ------------------------------- nmoymltrd = 0 tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0 tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0 tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 rmld_sum_trc (:,:) = 0.e0 ! ENDIF ! ====================================================================== ! V. Write restart file ! ====================================================================== IF( lrst_trc ) CALL trd_mld_trc_rst_write( kt ) ! this must be after the array swap above (III.3) END SUBROUTINE trd_mld_trc REAL FUNCTION sum2d( ztab ) !!---------------------------------------------------------------------- !! CD ??? prevoir d'utiliser plutot prtctl !!---------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab !!---------------------------------------------------------------------- sum2d = SUM(ztab(2:jpi-1,2:jpj-1)) END FUNCTION sum2d SUBROUTINE trd_mld_trc_init !!---------------------------------------------------------------------- !! *** ROUTINE trd_mld_init *** !! !! ** Purpose : computation of vertically integrated T and S budgets !! from ocean surface down to control surface (NetCDF output) !! !! ** Method/usage : !! !!---------------------------------------------------------------------- INTEGER :: ilseq, jl, jn REAL(wp) :: zjulian, zsto, zout CHARACTER (LEN=21) :: clold ='OLD' ! open specifier (direct access files) CHARACTER (LEN=21) :: clunf ='UNFORMATTED' ! open specifier (direct access files) CHARACTER (LEN=21) :: clseq ='SEQUENTIAL' ! open specifier (direct access files) CHARACTER (LEN=40) :: clop, cleiv CHARACTER (LEN=15) :: csuff CHARACTER (LEN=12) :: clmxl CHARACTER (LEN=16) :: cltrcu CHARACTER (LEN= 5) :: clvar CHARACTER (LEN=80) :: clname NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, & ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd !!---------------------------------------------------------------------- ! ====================================================================== ! I. initialization ! ====================================================================== IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' trd_mld_trc_init : Mixed-layer trends for passive tracers ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~' WRITE(numout,*) ENDIF ! I.1 Check consistency of user defined preferences ! ------------------------------------------------- #if defined key_trcldf_eiv IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN WRITE(numout,cform_war) WRITE(numout,*) ' You asked for ML diagnostics with iso-neutral diffusion ' WRITE(numout,*) ' and eiv physics. ' WRITE(numout,*) ' Yet, key_diaeiv is NOT switched on, so the eddy induced ' WRITE(numout,*) ' velocity is not diagnosed. ' WRITE(numout,*) ' Therefore, we cannot deduce the eiv advective trends. ' WRITE(numout,*) ' Only THE SUM of the i,j,k directional contributions then ' WRITE(numout,*) ' makes sense => To avoid any confusion, we choosed to mask ' WRITE(numout,*) ' these i,j,k directional contributions (with -999.) ' nwarn = nwarn + 1 ENDIF # endif IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_trc ) /= 0 ) ) THEN WRITE(numout,cform_err) WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend WRITE(numout,*) ' is no multiple of the trends diagnostics frequency ' WRITE(numout,*) ' you defined, ntrd_trc = ', ntrd_trc WRITE(numout,*) ' This will not allow you to restart from this simulation. ' WRITE(numout,*) ' You should reconsider this choice. ' WRITE(numout,*) WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a ' WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) ' nstop = nstop + 1 ENDIF IF( ( lk_trdmld_trc ) .AND. ( n_cla == 1 ) ) THEN WRITE(numout,cform_war) WRITE(numout,*) ' You set n_cla = 1. Note that the Mixed-Layer diagnostics ' WRITE(numout,*) ' are not exact along the corresponding straits. ' nwarn = nwarn + 1 ENDIF ! * Debugging information * IF( lldebug ) THEN WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl WRITE(numout,*) ' ln_trcadv_smolar = ' , ln_trcadv_smolar WRITE(numout,*) ' ln_trdmld_trc_instant = ', ln_trdmld_trc_instant ENDIF IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN WRITE(numout,cform_err) WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer Smolark. ' WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. ' WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and ' WRITE(numout,*) ' Smolarkiewicz scheme is Euler forward. ' WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However ' WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.' WRITE(numout,*) nstop = nstop + 1 ENDIF IF( ln_trcadv_muscl .AND. .NOT. ln_trdmld_trc_instant ) THEN WRITE(numout,cform_err) WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL ' WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. ' WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and ' WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note ' WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). ' WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However ' WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.' WRITE(numout,*) nstop = nstop + 1 ENDIF IF( ln_trcadv_muscl2 .AND. .NOT. ln_trdmld_trc_instant ) THEN WRITE(numout,cform_err) WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL2 ' WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. ' WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and ' WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note ' WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). ' WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However ' WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.' WRITE(numout,*) nstop = nstop + 1 ENDIF ! I.2 Initialize arrays to zero or read a restart file ! ---------------------------------------------------- nmoymltrd = 0 rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst. tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0 tmlradm_trc (:,:,:) = 0.e0 tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0 IF( lrsttr .AND. ln_trdmld_trc_restart ) THEN CALL trd_mld_trc_rst_read ELSE tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst. tmlbn_trc (:,:,:) = 0.e0 tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0 ENDIF ilseq = 1 ; icount = 1 ; ionce = 1 ! open specifier ! I.3 Read control surface from file ctlsurf_idx ! ---------------------------------------------- IF( nctls_trc == 1 ) THEN clname = 'ctlsurf_idx' CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 ) REWIND( numbol ) READ ( numbol ) nbol_trc ENDIF ! ====================================================================== ! II. netCDF output initialization ! ====================================================================== #if defined key_dimgout ??? #else ! clmxl = legend root for netCDF output IF( nctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion clmxl = 'Mixed Layer ' ELSE IF( nctls_trc == 1 ) THEN ! control surface = read index from file clmxl = ' Bowl ' ELSE IF( nctls_trc >= 2 ) THEN ! control surface = model level WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls_trc ENDIF ! II.1 Define frequency of output and means ! ----------------------------------------- # if defined key_diainstant IF( .NOT. ln_trdmld_trc_instant ) THEN STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...' ENDIF zsto = ntrd_trc * rdt clop ="inst(only(x))" # else IF( ln_trdmld_trc_instant ) THEN zsto = rdt ! inst. diags : we use IOIPSL time averaging ELSE zsto = ntrd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging ENDIF clop ="ave(only(x))" # endif zout = ntrd_trc * rdt IF(lwp) WRITE (numout,*) ' netCDF initialization' ! II.2 Compute julian date from starting date of the run ! ------------------------------------------------------ CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian ) IF(lwp) WRITE(numout,*)' ' IF(lwp) WRITE(numout,*)' Date 0 used :', nit000 & & ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday & & ,'Julian day : ', zjulian ! II.3 Define the T grid trend file (nidtrd) ! ------------------------------------------ !-- Define long and short names for the NetCDF output variables ! ==> choose them according to trdmld_trc_oce.F90 <== #if defined key_diaeiv cleiv = " (*** only total EIV is meaningful ***)" ! eiv advec. trends require u_eiv, v_eiv #else cleiv = " " #endif ctrd_trc(jpmld_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmld_trc_xad ,2) = "_xad" ctrd_trc(jpmld_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmld_trc_yad ,2) = "_yad" ctrd_trc(jpmld_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmld_trc_zad ,2) = "_zad" ctrd_trc(jpmld_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmld_trc_ldf ,2) = "_ldf" ctrd_trc(jpmld_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmld_trc_zdf ,2) = "_zdf" ctrd_trc(jpmld_trc_xei ,1) = " Zonal EIV advection"//cleiv ; ctrd_trc(jpmld_trc_xei ,2) = "_xei" ctrd_trc(jpmld_trc_yei ,1) = " Merid. EIV advection"//cleiv ; ctrd_trc(jpmld_trc_yei ,2) = "_yei" ctrd_trc(jpmld_trc_zei ,1) = " Vertical EIV advection"//cleiv ; ctrd_trc(jpmld_trc_zei ,2) = "_zei" ctrd_trc(jpmld_trc_bbc ,1) = " Geothermal flux" ; ctrd_trc(jpmld_trc_bbc ,2) = "_bbc" ctrd_trc(jpmld_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmld_trc_bbl ,2) = "_bbl" ctrd_trc(jpmld_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmld_trc_dmp ,2) = "_dmp" ctrd_trc(jpmld_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmld_trc_sbc ,2) = "_sbc" #if defined key_lobster ctrd_trc(jpmld_trc_sms_sed,1) = " Sources minus sinks : sed" ; ctrd_trc(jpmld_trc_sms_sed,2) = "_sms_sed" ctrd_trc(jpmld_trc_sms_bio,1) = " Sources minus sinks : bio" ; ctrd_trc(jpmld_trc_sms_bio,2) = "_sms_bio" ctrd_trc(jpmld_trc_sms_exp,1) = " Sources minus sinks : exp" ; ctrd_trc(jpmld_trc_sms_exp,2) = "_sms_exp" #else ctrd_trc(jpmld_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmld_trc_sms ,2) = "_sms" #endif ctrd_trc(jpmld_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmld_trc_radb ,2) = "_radb" ctrd_trc(jpmld_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmld_trc_radn ,2) = "_radn" ctrd_trc(jpmld_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmld_trc_atf ,2) = "_atf" ctrd_trc(jpltrd_trc+1 ,1) = " Total EIV"//cleiv ; ctrd_trc(jpltrd_trc+1 ,2) = "_tei" DO jn = 1, jptra !-- Create a NetCDF file and enter the define mode IF( luttrd(jn) ) THEN csuff="TD_"//ctrcnm(jn) CALL dia_nam( clhstnam, ntrd_trc, csuff ) CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & & 1, jpi, 1, jpj, 0, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom ) !-- Define the ML depth variable CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", & & jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ENDIF END DO !-- Define physical units IF( ucf_trc == 1. ) THEN cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit ELSEIF ( ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3) cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf ELSE cltrcu = "unknown?" ENDIF !-- Define miscellaneous passive tracer mixed-layer variables IF( jpltrd_trc /= jpmld_trc_atf .OR. jpltrd_trc - 1 /= jpmld_trc_radb ) THEN STOP 'Error : jpltrd_trc /= jpmld_trc_atf .OR. jpltrd_trc - 1 /= jpmld_trc_radb' ! see below ENDIF #if defined key_lobster IF( lldebug ) THEN DO jn = 1, jptra WRITE(numout, *) 'TRC jpdet=', jpdet, ' jpnh4=', jpnh4 WRITE(numout, *) 'TRC short title ctrcnm jn=", jn, " : ', ctrcnm(jn) WRITE(numout, *) 'TRC trim(ctrcnm(jn))//"_tot" = ', trim(ctrcnm(jn))//"ml_tot" ! tml_tot -> detml_tot END DO CALL flush(numout) ENDIF #else !! Error : this is not ready (PISCES) #endif DO jn = 1, jptra ! IF( luttrd(jn) ) THEN clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc. CALL histdef(nidtrd(jn), clvar, clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", & & "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) CALL histdef(nidtrd(jn), clvar//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", & & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) CALL histdef(nidtrd(jn), clvar//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", & & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmld_trc_atf CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), & & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean END DO ! if zsto=rdt above CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_radb,1), & & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1), & & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpltrd_trc+1,2)), clmxl//" "//clvar//ctrd_trc(jpltrd_trc+1 ,1), & & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! Total EIV ! ENDIF END DO !-- Leave IOIPSL/NetCDF define mode DO jn = 1, jptra IF( luttrd(jn) ) CALL histend( nidtrd(jn) ) END DO #endif /* key_dimgout */ END SUBROUTINE trd_mld_trc_init #else !!---------------------------------------------------------------------- !! Default option : Empty module !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trd_mld_trc( kt ) ! Empty routine INTEGER, INTENT( in) :: kt WRITE(*,*) 'trd_mld_trc: You should not have seen this print! error?', kt END SUBROUTINE trd_mld_trc SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn ) INTEGER, INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmld ! passive trc trend WRITE(*,*) 'trd_mld_trc_zint: You should not have seen this print! error?', ptrc_trdmld(1,1,1) WRITE(*,*) ' " " : You should not have seen this print! error?', ctype WRITE(*,*) ' " " : You should not have seen this print! error?', ktrd END SUBROUTINE trd_mld_trc_zint SUBROUTINE trd_mld_trc_init ! Empty routine WRITE(*,*) 'trd_mld_trc_init: You should not have seen this print! error?' END SUBROUTINE trd_mld_trc_init #endif !!====================================================================== END MODULE trdmld_trc