Changeset 503 for trunk/NEMO/OPA_SRC/TRD/trdmld.F90
- Timestamp:
- 2006-09-27T10:52:29+02:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/TRD/trdmld.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/TRD/trdmld.F90
r352 r503 4 4 !! Ocean diagnostics: mixed layer T-S trends 5 5 !!===================================================================== 6 !! History : ! 95-04 (J. Vialard) Original code 7 !! ! 97-02 (E. Guilyardi) Adaptation global + base cmo 8 !! ! 99-09 (E. Guilyardi) Re-writing + netCDF output 9 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 10 !! 9.0 ! 04-08 (C. Talandier) New trends organization 11 !! ! 05-05 (C. Deltel) Diagnose trends of time averaged ML T & S 12 !!---------------------------------------------------------------------- 6 13 #if defined key_trdmld || defined key_esopa 7 14 !!---------------------------------------------------------------------- 8 15 !! 'key_trdmld' mixed layer trend diagnostics 16 !!---------------------------------------------------------------------- 9 17 !!---------------------------------------------------------------------- 10 18 !! trd_mld : T and S cumulated trends averaged over the mixed layer … … 12 20 !! trd_mld_init : initialization step 13 21 !!---------------------------------------------------------------------- 14 !! * Modules used15 22 USE oce ! ocean dynamics and tracers variables 16 23 USE dom_oce ! ocean space and time domain variables … … 28 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 36 USE diadimg ! dimg direct access file format output 37 USE trdmld_rst , ONLY : trd_mld_rst_read ! restart for diagnosing the ML trends 38 USE prtctl ! Print control 30 39 31 40 IMPLICIT NONE 32 41 PRIVATE 33 42 34 !! * Accessibility 35 PUBLIC trd_mld ! routine called by step.F90 36 PUBLIC trd_mld_init ! routine called by opa.F90 37 PUBLIC trd_mld_zint ! routine called by tracers routines 38 39 !! * Shared module variables 40 LOGICAL, PUBLIC, PARAMETER :: lk_trdmld = .TRUE. !: momentum trend flag 41 42 !! * Module variables 43 INTEGER :: & 44 nh_t, nmoymltrd, & ! ??? 45 nidtrd, & 46 ndextrd1(jpi*jpj), & 47 ndimtrd1 48 INTEGER, SAVE :: & 49 ionce, icount, & 50 idebug ! (0/1) set it to 1 in case of problem to have more print 51 52 INTEGER, DIMENSION(jpi,jpj) :: & 53 nmld, & ! mixed layer depth 54 nbol 55 56 REAL(wp), DIMENSION(jpi,jpj) :: & 57 rmld , & ! mld depth (m) corresponding to nmld 58 tml , sml , & ! average T and S over mixed layer 59 tmlb , smlb , & ! before tml and sml (kt-1) 60 tmlbb , smlbb, & ! tml and sml at begining of the nwrite-1 61 ! ! timestep averaging period 62 tmlbn , smlbn, & ! after tml and sml at time step after the 63 ! ! begining of the NWRITE-1 timesteps 64 tmltrdm, smltrdm ! 65 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 67 tmltrd , & ! total cumulative trends of temperature and 68 smltrd , & ! salinity over nwrite-1 time steps 69 wkx 70 71 CHARACTER(LEN=80) :: clname 43 PUBLIC trd_mld ! routine called by step.F90 44 PUBLIC trd_mld_init ! routine called by opa.F90 45 PUBLIC trd_mld_zint ! routine called by tracers routines 46 47 CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file 48 INTEGER :: nh_t, nmoymltrd 49 INTEGER :: nidtrd, ndextrd1(jpi*jpj) 50 INTEGER :: ndimtrd1 51 INTEGER, SAVE :: ionce, icount 72 52 73 53 !! * Substitutions … … 83 63 CONTAINS 84 64 85 SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )65 SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype ) 86 66 !!---------------------------------------------------------------------- 87 67 !! *** ROUTINE trd_mld_zint *** 88 68 !! 89 !! ** Purpose : computation of vertically integrated T and S budgets 90 !! from ocean surface down to control surface 69 !! ** Purpose : Compute the vertical average of the 3D fields given as arguments 70 !! to the subroutine. This vertical average is performed from ocean 71 !! surface down to a chosen control surface. 91 72 !! 92 73 !! ** Method/usage : 93 !! integration done over nwrite-1 time steps 94 !! Control surface can be either a mixed layer depth (time varying) 74 !! The control surface can be either a mixed layer depth (time varying) 95 75 !! or a fixed surface (jk level or bowl). 96 !! Choose control surface with nctls in namelist NAM DIA.97 !! nctls = 0 : use mixed layer with density criterion98 !! nctls = 1 : read index from file 'ctlsurf_idx'99 !! nctls > 1 : use fixed level surface jk = nctls76 !! Choose control surface with nctls in namelist NAMTRD : 77 !! nctls = 0 : use mixed layer with density criterion 78 !! nctls = 1 : read index from file 'ctlsurf_idx' 79 !! nctls > 1 : use fixed level surface jk = nctls 100 80 !! Note: in the remainder of the routine, the volume between the 101 81 !! surface and the control surface is called "mixed-layer" 102 !! Method check : if the control surface is fixed, the residual dh/dt103 !! entrainment should be zero104 !!105 !! ** Action :106 !! /commld/ : rmld mld depth corresponding to nmld107 !! tml average T over mixed layer108 !! tmlb tml at kt-1109 !! tmlbb tml at begining of the NWRITE-1110 !! time steps averaging period111 !! tmlbn tml at time step after the112 !! begining of the NWRITE-1 time113 !! steps averaging period114 !!115 !! mixed layer trends :116 !!117 !! tmltrd (,,1) = zonal advection118 !! tmltrd (,,2) = meridional advection119 !! tmltrd (,,3) = vertical advection120 !! tmltrd (,,4) = lateral diffusion (horiz. component+Beckman)121 !! tmltrd (,,5) = forcing122 !! tmltrd (,,6) = entrainment due to vertical diffusion (TKE)123 !! if iso tmltrd (,,7) = lateral diffusion (vertical component)124 !! tmltrd (,,8) = eddy induced zonal advection125 !! tmltrd (,,9) = eddy induced meridional advection126 !! tmltrd (,,10) = eddy induced vertical advection127 !!128 !! tmltrdm(,) : total cumulative trends over nwrite-1 time steps129 !! ztmltot(,) : dT/dt over the NWRITE-1 time steps130 !! averaging period (including Asselin131 !! terms)132 !! ztmlres(,) : residual = dh/dt entrainment133 !!134 !! trends output in netCDF format using ioipsl135 !!136 !! History :137 !! ! 95-04 (J. Vialard) Original code138 !! ! 97-02 (E. Guilyardi) Adaptation global + base cmo139 !! ! 99-09 (E. Guilyardi) Re-writing + netCDF output140 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module141 !! 9.0 ! 04-08 (C. Talandier) New trends organization142 82 !!---------------------------------------------------------------------- 143 !! * Arguments 144 INTEGER, INTENT( in ) :: ktrd ! ocean trend index 145 146 CHARACTER(len=2), INTENT( in ) :: & 147 ctype ! surface/bottom (2D arrays) or 148 ! interior (3D arrays) physics 149 150 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & 151 pttrdmld, & ! Temperature trend 152 pstrdmld ! Salinity trend 153 154 !! * Local declarations 83 INTEGER, INTENT( in ) :: ktrd ! ocean trend index 84 CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D arrays) or 85 ! ! interior (3D arrays) physics 86 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pttrdmld ! temperature trend 87 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pstrdmld ! salinity trend 155 88 INTEGER :: ji, jj, jk, isum 156 # if defined key_trabbl_dif 157 INTEGER :: ikb 158 # endif 159 160 REAL(wp), DIMENSION(jpi,jpj) :: & 161 zvlmsk 89 REAL(wp), DIMENSION(jpi,jpj) :: zvlmsk 162 90 !!---------------------------------------------------------------------- 163 91 92 ! I. Definition of control surface and associated fields 93 ! ------------------------------------------------------ 94 ! ==> only once per time step <== 95 164 96 IF( icount == 1 ) THEN 165 166 zvlmsk(:,:) = 0.e0 167 tmltrd(:,:,:) = 0.e0 168 smltrd(:,:,:) = 0.e0 169 170 ! This computation should be done only once per time step 171 172 ! ======================================================== 173 ! I. definition of control surface and associated fields 174 ! ======================================================== 175 176 ! I.1 set nmld(ji,jj) = index of first T point below control surface 177 ! ------------------- or outside mixed-layer 178 179 IF( nctls == 0 ) THEN 180 ! control surface = mixed-layer with density criterion 181 ! (array nmln computed in zdfmxl.F90) 182 nmld(:,:) = nmln(:,:) 183 ELSE IF( nctls == 1 ) THEN 184 ! control surface = read index from file 97 ! 98 tmltrd(:,:,:) = 0.e0 ; smltrd(:,:,:) = 0.e0 ! <<< reset trend arrays to zero 99 100 ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer 101 IF( nctls == 0 ) THEN ! * control surface = mixed-layer with density criterion 102 nmld(:,:) = nmln(:,:) ! array nmln computed in zdfmxl.F90 103 ELSE IF( nctls == 1 ) THEN ! * control surface = read index from file 185 104 nmld(:,:) = nbol(:,:) 186 ELSE IF( nctls >= 2 ) THEN 187 ! control surface = model level 105 ELSE IF( nctls >= 2 ) THEN ! * control surface = model level 188 106 nctls = MIN( nctls, jpktrd - 1 ) 189 107 nmld(:,:) = nctls + 1 190 108 ENDIF 191 109 192 IF( ionce == 1 ) THEN ! compute ndextrd1 and ndimtrd1 only once 193 ! Check of validity : nmld(ji,jj) =< jpktrd 194 isum = 0 110 ! ... Compute ndextrd1 and ndimtrd1 only once 111 IF( ionce == 1 ) THEN 112 ! 113 ! Check of validity : nmld(ji,jj) <= jpktrd 114 isum = 0 115 zvlmsk(:,:) = 0.e0 195 116 196 117 IF( jpktrd < jpk ) THEN … … 215 136 ENDIF 216 137 217 ! no more pass here 218 ionce = 0 219 220 ENDIF 221 222 IF( idebug /= 0 ) THEN 223 ! CALL prihre (zvlmsk,jpi,jpj,1,jpi,2,1,jpj,2,3,numout) 224 WRITE(numout,*) ' debuging trd_mld_zint: I.1 done ' 225 CALL FLUSH(numout) 226 ENDIF 227 228 229 ! I.2 probability density function of presence in mixed-layer 230 ! -------------------------------- 231 ! (i.e. weight of each grid point in vertical integration : wkx(ji,jj,jk) 232 233 234 ! initialize wkx with vertical scale factor in mixed-layer 235 138 ionce = 0 ! no more pass here 139 ! 140 END IF 141 142 ! ... Weights for vertical averaging 236 143 wkx(:,:,:) = 0.e0 237 DO jk = 1, jpktrd 144 DO jk = 1, jpktrd ! initialize wkx with vertical scale factor in mixed-layer 238 145 DO jj = 1,jpj 239 146 DO ji = 1,jpi 240 IF( jk - nmld(ji,jj) < 0. ) wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)147 IF( jk - nmld(ji,jj) < 0.e0 ) wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk) 241 148 END DO 242 149 END DO 243 150 END DO 244 151 245 ! compute mixed-layer depth : rmld 246 247 rmld(:,:) = 0. 152 rmld(:,:) = 0.e0 ! compute mixed-layer depth : rmld 248 153 DO jk = 1, jpktrd 249 154 rmld(:,:) = rmld(:,:) + wkx(:,:,jk) 250 155 END DO 251 156 252 ! compute PDF 253 157 DO jk = 1, jpktrd ! compute integration weights 158 wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1., rmld(:,:) ) 159 END DO 160 161 icount = 0 ! <<< flag = off : control surface & integr. weights 162 ! ! computed only once per time step 163 END IF 164 165 ! II. Vertical integration of trends in the mixed-layer 166 ! ----------------------------------------------------- 167 168 SELECT CASE (ctype) 169 CASE ( '3D' ) ! mean T/S trends in the mixed-layer 254 170 DO jk = 1, jpktrd 255 wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1., rmld(:,:) ) 256 END DO 257 258 IF( idebug /= 0 ) THEN 259 WRITE(numout,*) ' debuging trd_mld_zint: I.2 done ' 260 CALL FLUSH(numout) 261 ENDIF 262 263 ! Set counter icount to 0 to avoid this part at each time step 264 icount = 0 265 266 ENDIF 267 268 269 ! ==================================================== 270 ! II. vertical integration of trends in mixed-layer 271 ! ==================================================== 272 273 ! II.1 vertical integration of 3D and 2D trends 274 ! --------------------------------------------- 275 276 SELECT CASE (ctype) 277 278 CASE ('3D') ! 3D treatment 279 280 ! trends terms in the mixed-layer 281 DO jk = 1, jpktrd 282 ! Temperature 283 tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,jk) * wkx(:,:,jk) 284 285 ! Salinity 286 smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,jk) * wkx(:,:,jk) 287 ENDDO 288 289 CASE ('2D') ! 2D treatment 290 291 SELECT CASE (ktrd) 292 293 CASE (jpmldldf) 294 295 # if defined key_trabbl_dif 296 ! trends terms from Beckman over-flow parameterization 297 DO jj = 1,jpj 298 DO ji = 1,jpi 299 ikb = MAX( mbathy(ji,jj)-1, 1 ) 300 ! beckmann component -> horiz. part of lateral diffusion 301 tmltrd(ji,jj,ktrd) = tmltrd(ji,jj,ktrd) + pttrdmld(ji,jj,1) * wkx(ji,jj,ikb) 302 smltrd(ji,jj,ktrd) = smltrd(ji,jj,ktrd) + pstrdmld(ji,jj,1) * wkx(ji,jj,ikb) 303 END DO 304 END DO 305 # endif 306 307 CASE DEFAULT 308 309 ! trends terms at upper boundary of mixed-layer 310 311 ! forcing term (non penetrative) 312 ! Temperature 313 tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,1) * wkx(:,:,1) 314 315 ! forcing term 316 ! Salinity 317 smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1) 318 319 END SELECT 320 171 tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,jk) * wkx(:,:,jk) ! temperature 172 smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,jk) * wkx(:,:,jk) ! salinity 173 END DO 174 CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer 175 tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,1) * wkx(:,:,1) ! non penetrative 176 smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1) 321 177 END SELECT 322 323 IF( idebug /= 0 ) THEN 324 IF(lwp) WRITE(numout,*) ' debuging trd_mld_zint: II.1 done' 325 CALL FLUSH(numout) 326 ENDIF 327 178 ! 328 179 END SUBROUTINE trd_mld_zint 329 330 180 331 181 332 182 SUBROUTINE trd_mld( kt ) … … 334 184 !! *** ROUTINE trd_mld *** 335 185 !! 336 !! ** Purpose : computation of cumulated trends over analysis period337 !! and make outputs (NetCDF or DIMG format)186 !! ** Purpose : Compute and cumulate the mixed layer trends over an analysis 187 !! period, and write NetCDF (or dimg) outputs. 338 188 !! 339 189 !! ** Method/usage : 190 !! The stored trends can be chosen twofold (according to the ln_trdmld_instant 191 !! logical namelist variable) : 192 !! 1) to explain the difference between initial and final 193 !! mixed-layer T & S (where initial and final relate to the 194 !! current analysis window, defined by ntrd in the namelist) 195 !! 2) to explain the difference between the current and previous 196 !! TIME-AVERAGED mixed-layer T & S (where time-averaging is 197 !! performed over each analysis window). 340 198 !! 341 !! History : 342 !! 9.0 ! 04-08 (C. Talandier) New trends organization 199 !! ** Consistency check : 200 !! If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt 201 !! entrainment) should be zero, at machine accuracy. Note that in the case 202 !! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO 203 !! over the first two analysis windows (except if restart). 204 !! N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8 205 !! for checking residuals. 206 !! On a NEC-SX5 computer, this typically leads to: 207 !! O(1.e-20) temp. residuals (tml_res) when ln_trdmld_instant=.false. 208 !! O(1.e-21) temp. residuals (tml_res) when ln_trdmld_instant=.true. 209 !! 210 !! ** Action : 211 !! At each time step, mixed-layer averaged trends are stored in the 212 !! tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx). 213 !! This array is known when trd_mld is called, at the end of the stp subroutine, 214 !! except for the purely vertical K_z diffusion term, which is embedded in the 215 !! lateral diffusion trend. 216 !! 217 !! In I), this K_z term is diagnosed and stored, thus its contribution is removed 218 !! from the lateral diffusion trend. 219 !! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative 220 !! arrays are updated. 221 !! In III), called only once per analysis window, we compute the total trends, 222 !! along with the residuals and the Asselin correction terms. 223 !! In IV), the appropriate trends are written in the trends NetCDF file. 224 !! 225 !! References : 226 !! - Vialard & al. 227 !! - See NEMO documentation (in preparation) 343 228 !!---------------------------------------------------------------------- 344 !! * Arguments345 229 INTEGER, INTENT( in ) :: kt ! ocean time-step index 346 347 !! * Local declarations 230 !! 348 231 INTEGER :: ji, jj, jk, jl, ik, it 349 350 REAL(wp) :: zmean, zavt 351 352 REAL(wp) ,DIMENSION(jpi,jpj) :: & 353 ztmltot, ztmlres, & 354 zsmltot, zsmlres, & 355 z2d 356 232 LOGICAL :: lldebug = .TRUE. 233 REAL(wp) :: zavt, zfn, zfn2 234 REAL(wp) ,DIMENSION(jpi,jpj) :: & 235 ztmltot, zsmltot, & ! dT/dt over the anlysis window (including Asselin) 236 ztmlres, zsmlres, & ! residual = dh/dt entrainment term 237 ztmlatf, zsmlatf, & ! needed for storage only 238 ztmltot2, ztmlres2, ztmltrdm2, & ! \ working arrays to diagnose the trends 239 zsmltot2, zsmlres2, zsmltrdm2, & ! > associated with the time meaned ML T & S 240 ztmlatf2, zsmlatf2 ! / 241 REAL(wp), DIMENSION(jpi,jpj,jpltrd) :: & 242 ztmltrd2, zsmltrd2 ! only needed for mean diagnostics 357 243 #if defined key_dimgout 358 244 INTEGER :: iyear,imon,iday … … 361 247 !!---------------------------------------------------------------------- 362 248 363 ! I. trends terms at lower boundary of mixed-layer 364 ! ------------------------------------------------ 365 249 ! ====================================================================== 250 ! I. Diagnose the purely vertical (K_z) diffusion trend 251 ! ====================================================================== 252 253 ! ... These terms can be estimated by flux computation at the lower boundary of the ML 254 ! (we compute (-1/h) * K_z * d_z( T ) and (-1/h) * K_z * d_z( S )) 366 255 DO jj = 1,jpj 367 256 DO ji = 1,jpi 368 369 257 ik = nmld(ji,jj) 370 371 ! Temperature372 ! entrainment due to vertical diffusion373 ! - due to vertical mixing scheme (TKE)374 258 zavt = avt(ji,jj,ik) 375 tmltrd(ji,jj,jpmldevd) = - 1. * zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & 376 & * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) ) & 377 & / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1) 378 ! Salinity 379 ! entrainment due to vertical diffusion 380 ! - due to vertical mixing scheme (TKE) 259 tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & 260 & * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) ) & 261 & / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1) 381 262 zavt = fsavs(ji,jj,ik) 382 smltrd(ji,jj,jpmld evd) = -1. * zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)&383 & * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) )&384 & / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)263 smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & 264 & * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) ) & 265 & / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1) 385 266 END DO 386 267 END DO 387 268 269 ! ... Remove this K_z trend from the iso-neutral diffusion term (if any) 388 270 IF( ln_traldf_iso ) THEN 389 ! We substract to the TOTAL vertical diffusion tmltrd(:,:,jpmldzdf) 390 ! computed in subroutines trazdf_iso.F90 or trazdf_imp.F90 391 ! the vertical part du to the Kz in order to keep only the vertical 392 ! isopycnal diffusion (i.e the isopycnal diffusion componant on the vertical): 393 tmltrd(:,:,jpmldzdf) = tmltrd(:,:,jpmldzdf) - tmltrd(:,:,jpmldevd) ! - due to isopycnal mixing scheme (implicit part) 394 smltrd(:,:,jpmldzdf) = smltrd(:,:,jpmldzdf) - smltrd(:,:,jpmldevd) ! - due to isopycnal mixing scheme (implicit part) 395 ENDIF 396 397 ! Boundary conditions 398 CALL lbc_lnk( tmltrd, 'T', 1. ) 399 CALL lbc_lnk( smltrd, 'T', 1. ) 400 401 IF( idebug /= 0 ) THEN 402 WRITE(numout,*) ' debuging trd_mld: I. done' 403 CALL FLUSH(numout) 404 ENDIF 405 406 ! ================================= 407 ! II. Cumulated trends 408 ! ================================= 409 410 ! II.1 set before values of vertically average T and S 411 ! --------------------------------------------------- 412 271 tmltrd(:,:,jpmld_ldf) = tmltrd(:,:,jpmld_ldf) - tmltrd(:,:,jpmld_zdf) 272 smltrd(:,:,jpmld_ldf) = smltrd(:,:,jpmld_ldf) - smltrd(:,:,jpmld_zdf) 273 END IF 274 275 ! ... Lateral boundary conditions 276 DO jl = 1, jpltrd 277 CALL lbc_lnk( tmltrd(:,:,jl), 'T', 1. ) 278 CALL lbc_lnk( smltrd(:,:,jl), 'T', 1. ) 279 END DO 280 281 ! ====================================================================== 282 ! II. Cumulate the trends over the analysis window 283 ! ====================================================================== 284 285 ztmltrd2(:,:,:) = 0.e0 ; zsmltrd2(:,:,:) = 0.e0 ! <<< reset arrays to zero 286 ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0 287 ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0 288 ztmlatf2(:,:) = 0.e0 ; zsmlatf2(:,:) = 0.e0 289 290 ! II.1 Set before values of vertically average T and S 291 ! ---------------------------------------------------- 413 292 IF( kt > nit000 ) THEN 414 tmlb(:,:) = tml(:,:) 415 smlb(:,:) = sml(:,:) 416 ENDIF 417 418 ! II.2 vertically integrated T and S 419 ! --------------------------------- 420 421 tml(:,:) = 0. 422 sml(:,:) = 0. 423 293 ! ... temperature ... ... salinity ... 294 tmlb (:,:) = tml (:,:) ; smlb (:,:) = sml (:,:) 295 tmlatfn(:,:) = tmltrd(:,:,jpmld_atf) ; smlatfn(:,:) = smltrd(:,:,jpmld_atf) 296 END IF 297 298 ! II.2 Vertically averaged T and S 299 ! -------------------------------- 300 tml(:,:) = 0.e0 ; sml(:,:) = 0.e0 424 301 DO jk = 1, jpktrd - 1 425 302 tml(:,:) = tml(:,:) + wkx(:,:,jk) * tn(:,:,jk) … … 427 304 END DO 428 305 429 IF(idebug /= 0) THEN 430 WRITE(numout,*) ' debuging trd_mld: II.2 done' 431 CALL FLUSH(numout) 432 ENDIF 433 434 ! II.3 set `before' mixed layer values for kt = nit000+1 435 ! -------------------------------------------------------- 436 437 IF( kt == nit000+1 ) THEN 438 tmlbb(:,:) = tmlb(:,:) 439 tmlbn(:,:) = tml (:,:) 440 smlbb(:,:) = smlb(:,:) 441 smlbn(:,:) = sml (:,:) 442 ENDIF 443 444 IF( idebug /= 0 ) THEN 445 WRITE(numout,*) ' debuging trd_mld: II.3 done' 446 CALL FLUSH(numout) 447 ENDIF 448 449 ! II.4 cumulated trends over analysis period (kt=2 to nwrite) 450 ! ----------------------------------------------------------- 451 452 ! trends cumulated over nwrite-2 time steps 453 454 IF( kt >= nit000+2 ) THEN 306 ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window 307 ! ------------------------------------------------------------------------ 308 IF( kt == 2 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) 309 ! 310 ! ... temperature ... ... salinity ... 311 tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:) 312 tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:) 313 tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:) 314 315 tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0 316 tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0 317 318 rmldbn(:,:) = rmld(:,:) 319 320 IF( ln_ctl ) THEN 321 WRITE(numout,*) ' we reach kt == nit000 + 1 = ', nit000+1 322 CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1) 323 CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1) 324 CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1) 325 END IF 326 ! 327 END IF 328 329 IF( ( ln_rstart ) .AND. ( kt == nit000 ) .AND. ( ln_ctl ) ) THEN 330 IF( ln_trdmld_instant ) THEN 331 WRITE(numout,*) ' restart from kt == nit000 = ', nit000 332 CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1) 333 CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1) 334 CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1) 335 ELSE 336 WRITE(numout,*) ' restart from kt == nit000 = ', nit000 337 CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1) 338 CALL prt_ctl(tab2d_1=rmldbn , clinfo1=' rmldbn - : ', mask1=tmask, ovlap=1) 339 CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1) 340 CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1) 341 CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1) 342 END IF 343 END IF 344 345 ! II.4 Cumulated trends over the analysis period 346 ! ---------------------------------------------- 347 ! 348 ! [ 1rst analysis window ] [ 2nd analysis window ] 349 ! 350 ! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps 351 ! ntrd 2*ntrd etc. 352 ! 1 2 3 4 =5 e.g. =10 353 ! 354 IF( ( kt >= 2 ).OR.( ln_rstart ) ) THEN 355 ! 455 356 nmoymltrd = nmoymltrd + 1 357 358 ! ... Cumulate over BOTH physical contributions AND over time steps 456 359 DO jl = 1, jpltrd 457 360 tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl) 458 361 smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl) 459 362 END DO 460 ENDIF 461 462 IF( idebug /= 0 ) THEN 463 WRITE(numout,*) ' debuging trd_mld: II.4 done' 464 CALL FLUSH(numout) 465 ENDIF 466 467 ! ============================================= 468 ! III. Output in netCDF + residual computation 469 ! ============================================= 470 471 ztmltot(:,:) = 0. 472 zsmltot(:,:) = 0. 473 ztmlres(:,:) = 0. 474 zsmlres(:,:) = 0. 475 476 IF( MOD( kt - nit000+1, nwrite ) == 0 ) THEN 477 478 ! III.1 compute total trend 479 ! ------------------------ 480 481 zmean = float(nmoymltrd) 482 483 ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / (zmean * 2. * rdt) 484 zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / (zmean * 2. * rdt) 485 486 IF(idebug /= 0) THEN 487 WRITE(numout,*) ' zmean = ',zmean 488 WRITE(numout,*) ' debuging trd_mld: III.1 done' 489 CALL FLUSH(numout) 490 ENDIF 491 492 493 ! III.2 compute residual 494 ! --------------------- 495 496 ztmlres(:,:) = ztmltot(:,:) - tmltrdm(:,:) / zmean 497 zsmlres(:,:) = zsmltot(:,:) - smltrdm(:,:) / zmean 498 499 500 ! Boundary conditions 501 502 CALL lbc_lnk( ztmltot, 'T', 1. ) 503 CALL lbc_lnk( ztmlres, 'T', 1. ) 504 CALL lbc_lnk( zsmltot, 'T', 1. ) 505 CALL lbc_lnk( zsmlres, 'T', 1. ) 506 507 IF( idebug /= 0 ) THEN 508 WRITE(numout,*) ' debuging trd_mld: III.2 done' 509 CALL FLUSH(numout) 510 ENDIF 511 512 513 ! III.3 time evolution array swap 514 ! ------------------------------ 515 516 tmlbb(:,:) = tmlb(:,:) 517 tmlbn(:,:) = tml (:,:) 518 smlbb(:,:) = smlb(:,:) 519 smlbn(:,:) = sml (:,:) 520 521 IF( idebug /= 0 ) THEN 522 WRITE(numout,*) ' debuging trd_mld: III.3 done' 523 CALL FLUSH(numout) 524 ENDIF 525 526 527 ! III.4 zero cumulative array 528 ! --------------------------- 529 530 nmoymltrd = 0 531 532 tmltrdm(:,:) = 0. 533 smltrdm(:,:) = 0. 534 535 IF(idebug /= 0) THEN 536 WRITE(numout,*) ' debuging trd_mld: III.4 done' 537 CALL FLUSH(numout) 538 ENDIF 539 540 ENDIF 541 542 ! III.5 write trends to output 543 ! --------------------------- 363 364 ! ... Special handling of the Asselin trend 365 tmlatfm(:,:) = tmlatfm(:,:) + tmlatfn(:,:) 366 smlatfm(:,:) = smlatfm(:,:) + smlatfn(:,:) 367 368 ! ... Trends associated with the time mean of the ML T/S 369 tmltrd_sum (:,:,:) = tmltrd_sum (:,:,:) + tmltrd (:,:,:) ! tem 370 tmltrd_csum_ln(:,:,:) = tmltrd_csum_ln(:,:,:) + tmltrd_sum(:,:,:) 371 tml_sum (:,:) = tml_sum (:,:) + tml (:,:) 372 smltrd_sum (:,:,:) = smltrd_sum (:,:,:) + smltrd (:,:,:) ! sal 373 smltrd_csum_ln(:,:,:) = smltrd_csum_ln(:,:,:) + smltrd_sum(:,:,:) 374 sml_sum (:,:) = sml_sum (:,:) + sml (:,:) 375 rmld_sum (:,:) = rmld_sum (:,:) + rmld (:,:) ! rmld 376 ! 377 END IF 378 379 ! ====================================================================== 380 ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD) 381 ! ====================================================================== 382 383 ! Convert to appropriate physical units 384 ! N.B. It may be useful to check IOIPSL time averaging with : 385 ! tmltrd (:,:,:) = 1. ; smltrd (:,:,:) = 1. 386 tmltrd(:,:,:) = tmltrd(:,:,:) * ucf ! (actually needed for 1:jpltrd-1, but trdmld(:,:,jpltrd) 387 smltrd(:,:,:) = smltrd(:,:,:) * ucf ! is no longer used, and is reset to 0. at next time step) 388 389 MODULO_NTRD : IF( MOD( kt, ntrd ) == 0 ) THEN ! nitend MUST be multiple of ntrd 390 ! 391 ztmltot (:,:) = 0.e0 ; zsmltot (:,:) = 0.e0 ! reset arrays to zero 392 ztmlres (:,:) = 0.e0 ; zsmlres (:,:) = 0.e0 393 ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0 394 ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0 395 396 zfn = float(nmoymltrd) ; zfn2 = zfn * zfn 397 398 ! III.1 Prepare fields for output ("instantaneous" diagnostics) 399 ! ------------------------------------------------------------- 400 401 !-- Compute total trends 402 ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / ( 2.*rdt ) 403 zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / ( 2.*rdt ) 404 405 !-- Compute residuals 406 ztmlres(:,:) = ztmltot(:,:) - ( tmltrdm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) ) 407 zsmlres(:,:) = zsmltot(:,:) - ( smltrdm(:,:) - smlatfn(:,:) + smlatfb(:,:) ) 408 409 !-- Diagnose Asselin trend over the analysis window 410 ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) 411 zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:) 412 413 !-- Lateral boundary conditions 414 ! ... temperature ... ... salinity ... 415 CALL lbc_lnk( ztmltot , 'T', 1. ) ; CALL lbc_lnk( zsmltot , 'T', 1. ) 416 CALL lbc_lnk( ztmlres , 'T', 1. ) ; CALL lbc_lnk( zsmlres , 'T', 1. ) 417 CALL lbc_lnk( ztmlatf , 'T', 1. ) ; CALL lbc_lnk( zsmlatf , 'T', 1. ) 418 419 #if defined key_diainstant 420 CALL ctl_stop( 'tml_trd : key_diainstant was never checked within trdmld. Comment this to proceed.') 421 #endif 422 ! III.2 Prepare fields for output ("mean" diagnostics) 423 ! ---------------------------------------------------- 424 425 !-- Update the ML depth time sum (to build the Leap-Frog time mean) 426 rmld_sum(:,:) = rmldbn(:,:) + 2 * ( rmld_sum(:,:) - rmld(:,:) ) + rmld(:,:) 427 428 !-- Compute temperature total trends 429 tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:) 430 ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / ( 2.*rdt ) ! now in degC/s 431 432 !-- Compute salinity total trends 433 sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:) 434 zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / ( 2.*rdt ) ! now in psu/s 435 436 !-- Compute temperature residuals 437 DO jl = 1, jpltrd 438 ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl) 439 END DO 440 441 ztmltrdm2(:,:) = 0.e0 442 DO jl = 1, jpltrd 443 ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl) 444 END DO 445 446 ztmlres2(:,:) = ztmltot2(:,:) - & 447 ( ztmltrdm2(:,:) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:) ) 448 449 !-- Compute salinity residuals 450 DO jl = 1, jpltrd 451 zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl) 452 END DO 453 454 zsmltrdm2(:,:) = 0. 455 DO jl = 1, jpltrd 456 zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl) 457 END DO 458 459 zsmlres2(:,:) = zsmltot2(:,:) - & 460 ( zsmltrdm2(:,:) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:) ) 461 462 !-- Diagnose Asselin trend over the analysis window 463 ztmlatf2(:,:) = ztmltrd2(:,:,jpmld_atf) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:) 464 zsmlatf2(:,:) = zsmltrd2(:,:,jpmld_atf) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:) 465 466 !-- Lateral boundary conditions 467 ! ... temperature ... ... salinity ... 468 CALL lbc_lnk( ztmltot2, 'T', 1. ) ; CALL lbc_lnk( zsmltot2, 'T', 1. ) 469 CALL lbc_lnk( ztmlres2, 'T', 1. ) ; CALL lbc_lnk( zsmlres2, 'T', 1. ) 470 DO jl = 1, jpltrd 471 CALL lbc_lnk( ztmltrd2(:,:,jl), 'T', 1. ) ! \ these will be output 472 CALL lbc_lnk( zsmltrd2(:,:,jl), 'T', 1. ) ! / in the NetCDF trends file 473 END DO 474 475 ! III.3 Time evolution array swap 476 ! ------------------------------- 477 478 ! For T/S instantaneous diagnostics 479 ! ... temperature ... ... salinity ... 480 tmlbb (:,:) = tmlb (:,:) ; smlbb (:,:) = smlb (:,:) 481 tmlbn (:,:) = tml (:,:) ; smlbn (:,:) = sml (:,:) 482 tmlatfb(:,:) = tmlatfn(:,:) ; smlatfb(:,:) = smlatfn(:,:) 483 484 ! For T mean diagnostics 485 tmltrd_csum_ub (:,:,:) = zfn * tmltrd_sum(:,:,:) - tmltrd_csum_ln(:,:,:) 486 tml_sumb (:,:) = tml_sum(:,:) 487 tmltrd_atf_sumb(:,:) = tmltrd_sum(:,:,jpmld_atf) 488 489 ! For S mean diagnostics 490 smltrd_csum_ub (:,:,:) = zfn * smltrd_sum(:,:,:) - smltrd_csum_ln(:,:,:) 491 sml_sumb (:,:) = sml_sum(:,:) 492 smltrd_atf_sumb(:,:) = smltrd_sum(:,:,jpmld_atf) 493 494 ! ML depth 495 rmldbn (:,:) = rmld (:,:) 496 497 IF( ln_ctl ) THEN 498 IF( ln_trdmld_instant ) THEN 499 CALL prt_ctl(tab2d_1=tmlbb , clinfo1=' tmlbb - : ', mask1=tmask, ovlap=1) 500 CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1) 501 CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb - : ', mask1=tmask, ovlap=1) 502 ELSE 503 CALL prt_ctl(tab2d_1=tmlbn , clinfo1=' tmlbn - : ', mask1=tmask, ovlap=1) 504 CALL prt_ctl(tab2d_1=rmldbn , clinfo1=' rmldbn - : ', mask1=tmask, ovlap=1) 505 CALL prt_ctl(tab2d_1=tml_sumb , clinfo1=' tml_sumb - : ', mask1=tmask, ovlap=1) 506 CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb - : ', mask1=tmask, ovlap=1) 507 CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub - : ', mask1=tmask, ovlap=1, kdim=1) 508 END IF 509 END IF 510 511 ! III.4 Convert to appropriate physical units 512 ! ------------------------------------------- 513 514 ! ... temperature ... ... salinity ... 515 ztmltot (:,:) = ztmltot(:,:) * ucf/zfn ; zsmltot (:,:) = zsmltot(:,:) * ucf/zfn 516 ztmlres (:,:) = ztmlres(:,:) * ucf/zfn ; zsmlres (:,:) = zsmlres(:,:) * ucf/zfn 517 ztmlatf (:,:) = ztmlatf(:,:) * ucf/zfn ; zsmlatf (:,:) = zsmlatf(:,:) * ucf/zfn 518 519 tml_sum (:,:) = tml_sum (:,:) / (2*zfn) ; sml_sum (:,:) = sml_sum (:,:) / (2*zfn) 520 ztmltot2(:,:) = ztmltot2(:,:) * ucf/zfn2 ; zsmltot2(:,:) = zsmltot2(:,:) * ucf/zfn2 521 ztmltrd2(:,:,:) = ztmltrd2(:,:,:)* ucf/zfn2 ; zsmltrd2(:,:,:) = zsmltrd2(:,:,:)* ucf/zfn2 522 ztmlatf2(:,:) = ztmlatf2(:,:) * ucf/zfn2 ; zsmlatf2(:,:) = zsmlatf2(:,:) * ucf/zfn2 523 ztmlres2(:,:) = ztmlres2(:,:) * ucf/zfn2 ; zsmlres2(:,:) = zsmlres2(:,:) * ucf/zfn2 524 525 rmld_sum(:,:) = rmld_sum(:,:) / (2*zfn) ! similar to tml_sum and sml_sum 526 527 ! * Debugging information * 528 IF( lldebug ) THEN 529 ! 530 WRITE(numout,*) 531 WRITE(numout,*) 'trd_mld : write trends in the Mixed Layer for debugging process:' 532 WRITE(numout,*) '~~~~~~~ ' 533 WRITE(numout,*) ' TRA kt = ', kt, 'nmoymltrd = ', nmoymltrd 534 WRITE(numout,*) 535 WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA TEMPERATURE <<<<<<<<<<<<<<<<<<' 536 WRITE(numout,*) ' TRA ztmlres : ', SUM(ztmlres(:,:)) 537 WRITE(numout,*) ' TRA ztmltot : ', SUM(ztmltot(:,:)) 538 WRITE(numout,*) ' TRA tmltrdm : ', SUM(tmltrdm(:,:)) 539 WRITE(numout,*) ' TRA tmlatfb : ', SUM(tmlatfb(:,:)) 540 WRITE(numout,*) ' TRA tmlatfn : ', SUM(tmlatfn(:,:)) 541 DO jl = 1, jpltrd 542 WRITE(numout,*) ' * TRA TREND INDEX jpmld_xxx = jl = ', jl, & 543 & ' tmltrd : ', SUM(tmltrd(:,:,jl)) 544 END DO 545 WRITE(numout,*) ' TRA ztmlres (jpi/2,jpj/2) : ', ztmlres (jpi/2,jpj/2) 546 WRITE(numout,*) ' TRA ztmlres2(jpi/2,jpj/2) : ', ztmlres2(jpi/2,jpj/2) 547 WRITE(numout,*) 548 WRITE(numout,*) ' >>>>>>>>>>>>>>>>>> TRA SALINITY <<<<<<<<<<<<<<<<<<' 549 WRITE(numout,*) ' TRA zsmlres : ', SUM(zsmlres(:,:)) 550 WRITE(numout,*) ' TRA zsmltot : ', SUM(zsmltot(:,:)) 551 WRITE(numout,*) ' TRA smltrdm : ', SUM(smltrdm(:,:)) 552 WRITE(numout,*) ' TRA smlatfb : ', SUM(smlatfb(:,:)) 553 WRITE(numout,*) ' TRA smlatfn : ', SUM(smlatfn(:,:)) 554 DO jl = 1, jpltrd 555 WRITE(numout,*) ' * TRA TREND INDEX jpmld_xxx = jl = ', jl, & 556 & ' smltrd : ', SUM(smltrd(:,:,jl)) 557 END DO 558 WRITE(numout,*) ' TRA zsmlres (jpi/2,jpj/2) : ', zsmlres (jpi/2,jpj/2) 559 WRITE(numout,*) ' TRA zsmlres2(jpi/2,jpj/2) : ', zsmlres2(jpi/2,jpj/2) 560 ! 561 END IF 562 ! 563 END IF MODULO_NTRD 564 565 ! ====================================================================== 566 ! IV. Write trends in the NetCDF file 567 ! ====================================================================== 568 569 ! IV.1 Code for dimg mpp output 570 ! ----------------------------- 544 571 545 572 #if defined key_dimgout 546 ! code for dimg mpp output 547 IF ( MOD(kt,nwrite) == 0 ) THEN 548 WRITE(clmode,'(f5.1,a)' ) nwrite*rdt/86400.,' days average' 549 iyear = ndastp/10000 550 imon = (ndastp-iyear*10000)/100 551 iday = ndastp - imon*100 - iyear*10000 573 574 IF( MOD( kt, ntrd ) == 0 ) THEN 575 iyear = ndastp/10000 576 imon = (ndastp-iyear*10000)/100 577 iday = ndastp - imon*100 - iyear*10000 552 578 WRITE(clname,9000) TRIM(cexper),'MLDiags',iyear,imon,iday 553 cltext=TRIM(cexper)//' mld diags'//TRIM(clmode) 579 WRITE(clmode,'(f5.1,a)') ntrd*rdt/86400.,' days average' 580 cltext = TRIM(cexper)//' mld diags'//TRIM(clmode) 554 581 CALL dia_wri_dimg (clname, cltext, smltrd, jpltrd, '2') 555 9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc") 556 END IF 582 END IF 583 584 9000 FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc") 557 585 558 586 #else 559 IF( kt >= nit000+1 ) THEN 560 561 ! define time axis 562 it= kt-nit000+1 563 IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN 564 WRITE(numout,*) ' trd_mld : write NetCDF fields' 565 ENDIF 566 567 CALL histwrite( nidtrd,"somlttml",it,rmld ,ndimtrd1,ndextrd1) ! Mixed-layer depth 568 569 ! Temperature trends 570 ! ------------------ 571 CALL histwrite( nidtrd,"somltemp",it,tml ,ndimtrd1,ndextrd1) ! Mixed-layer temperature 572 CALL histwrite( nidtrd,"somlttto",it,ztmltot ,ndimtrd1,ndextrd1) ! total 573 CALL histwrite( nidtrd,"somlttax",it,tmltrd(:,:, 1),ndimtrd1,ndextrd1) ! i- adv. 574 CALL histwrite( nidtrd,"somlttay",it,tmltrd(:,:, 2),ndimtrd1,ndextrd1) ! j- adv. 575 CALL histwrite( nidtrd,"somlttaz",it,tmltrd(:,:, 3),ndimtrd1,ndextrd1) ! vertical adv. 576 CALL histwrite( nidtrd,"somlttdh",it,tmltrd(:,:, 4),ndimtrd1,ndextrd1) ! hor. lateral diff. 577 CALL histwrite( nidtrd,"somlttfo",it,tmltrd(:,:, 5),ndimtrd1,ndextrd1) ! forcing 578 579 CALL histwrite( nidtrd,"somlbtdz",it,tmltrd(:,:, 6),ndimtrd1,ndextrd1) ! vert. diffusion 580 CALL histwrite( nidtrd,"somlbtdt",it,ztmlres ,ndimtrd1,ndextrd1) ! dh/dt entrainment (residual) 581 IF( ln_traldf_iso ) THEN 582 CALL histwrite( nidtrd,"somlbtdv",it,tmltrd(:,:, 7),ndimtrd1,ndextrd1) ! vert. lateral diff. 583 ENDIF 584 #if defined key_traldf_eiv 585 CALL histwrite( nidtrd,"somlgtax",it,tmltrd(:,:, 8),ndimtrd1,ndextrd1) ! i- adv. (eiv) 586 CALL histwrite( nidtrd,"somlgtay",it,tmltrd(:,:, 9),ndimtrd1,ndextrd1) ! j- adv. (eiv) 587 CALL histwrite( nidtrd,"somlgtaz",it,tmltrd(:,:,10),ndimtrd1,ndextrd1) ! vert. adv. (eiv) 588 z2d(:,:) = tmltrd(:,:,8) + tmltrd(:,:,9) + tmltrd(:,:,10) 589 CALL histwrite( nidtrd,"somlgtat",it,z2d ,ndimtrd1,ndextrd1) ! total adv. (eiv) 590 #endif 591 592 ! Salinity trends 593 ! --------------- 594 CALL histwrite( nidtrd,"somlsalt",it,sml ,ndimtrd1,ndextrd1) ! Mixed-layer salinity 595 CALL histwrite( nidtrd,"somltsto",it,zsmltot ,ndimtrd1,ndextrd1) ! total 596 CALL histwrite( nidtrd,"somltsax",it,smltrd(:,:, 1),ndimtrd1,ndextrd1) ! i- adv. 597 CALL histwrite( nidtrd,"somltsay",it,smltrd(:,:, 2),ndimtrd1,ndextrd1) ! j- adv. 598 CALL histwrite( nidtrd,"somltsaz",it,smltrd(:,:, 3),ndimtrd1,ndextrd1) ! vert. adv. 599 CALL histwrite( nidtrd,"somltsdh",it,smltrd(:,:, 4),ndimtrd1,ndextrd1) ! hor. lateral diff. 600 CALL histwrite( nidtrd,"somltsfo",it,smltrd(:,:, 5),ndimtrd1,ndextrd1) ! forcing 601 CALL histwrite( nidtrd,"somlbsdz",it,smltrd(:,:, 6),ndimtrd1,ndextrd1) ! vert. diff. 602 CALL histwrite( nidtrd,"somlbsdt",it,zsmlres ,ndimtrd1,ndextrd1) ! dh/dt entrainment (residual) 603 IF( ln_traldf_iso ) THEN 604 CALL histwrite( nidtrd,"somlbsdv",it,smltrd(:,:, 7),ndimtrd1,ndextrd1) ! vert. lateral diff; 605 ENDIF 606 #if defined key_traldf_eiv 607 CALL histwrite( nidtrd,"somlgsax",it,smltrd(:,:, 8),ndimtrd1,ndextrd1) ! i-adv. (eiv) 608 CALL histwrite( nidtrd,"somlgsay",it,smltrd(:,:, 9),ndimtrd1,ndextrd1) ! j-adv. (eiv) 609 CALL histwrite( nidtrd,"somlgsaz",it,smltrd(:,:,10),ndimtrd1,ndextrd1) ! vert. adv. (eiv) 610 z2d(:,:) = smltrd(:,:,8) + smltrd(:,:,9) + smltrd(:,:,10) 611 CALL histwrite( nidtrd,"somlgsat",it,z2d ,ndimtrd1,ndextrd1) ! total adv. (eiv) 587 588 ! IV.2 Code for IOIPSL/NetCDF output 589 ! ---------------------------------- 590 591 IF( lwp .AND. MOD( kt , ntrd ) == 0 ) THEN 592 WRITE(numout,*) ' ' 593 WRITE(numout,*) 'trd_mld : write trends in the NetCDF file :' 594 WRITE(numout,*) '~~~~~~~ ' 595 WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt 596 WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd 597 WRITE(numout,*) ' ' 598 END IF 599 600 it = kt - nit000 + 1 601 602 !-- Write the trends for T/S instantaneous diagnostics 603 IF( ln_trdmld_instant ) THEN 604 605 CALL histwrite( nidtrd, "mxl_depth", it, rmld(:,:), ndimtrd1, ndextrd1 ) 606 607 !................................. ( ML temperature ) ................................... 608 609 !-- Output the fields 610 CALL histwrite( nidtrd, "tml" , it, tml (:,:), ndimtrd1, ndextrd1 ) 611 CALL histwrite( nidtrd, "tml_tot" , it, ztmltot(:,:), ndimtrd1, ndextrd1 ) 612 CALL histwrite( nidtrd, "tml_res" , it, ztmlres(:,:), ndimtrd1, ndextrd1 ) 613 614 DO jl = 1, jpltrd - 1 615 CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)), & 616 & it, tmltrd (:,:,jl), ndimtrd1, ndextrd1 ) 617 END DO 618 619 CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)), & 620 & it, ztmlatf(:,:), ndimtrd1, ndextrd1 ) 621 622 !.................................. ( ML salinity ) ..................................... 623 624 !-- Output the fields 625 CALL histwrite( nidtrd, "sml" , it, sml (:,:), ndimtrd1, ndextrd1 ) 626 CALL histwrite( nidtrd, "sml_tot" , it, zsmltot(:,:), ndimtrd1, ndextrd1 ) 627 CALL histwrite( nidtrd, "sml_res" , it, zsmlres(:,:), ndimtrd1, ndextrd1 ) 628 629 DO jl = 1, jpltrd - 1 630 CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)), & 631 & it, smltrd(:,:,jl), ndimtrd1, ndextrd1 ) 632 END DO 633 634 CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)), & 635 & it, zsmlatf(:,:), ndimtrd1, ndextrd1 ) 636 637 IF( kt == nitend ) CALL histclo( nidtrd ) 638 639 !-- Write the trends for T/S mean diagnostics 640 ELSE 641 642 CALL histwrite( nidtrd, "mxl_depth", it, rmld_sum(:,:), ndimtrd1, ndextrd1 ) 643 644 !................................. ( ML temperature ) ................................... 645 646 !-- Output the fields 647 CALL histwrite( nidtrd, "tml" , it, tml_sum (:,:), ndimtrd1, ndextrd1 ) 648 CALL histwrite( nidtrd, "tml_tot" , it, ztmltot2(:,:), ndimtrd1, ndextrd1 ) 649 CALL histwrite( nidtrd, "tml_res" , it, ztmlres2(:,:), ndimtrd1, ndextrd1 ) 650 651 DO jl = 1, jpltrd - 1 652 CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)), & 653 & it, ztmltrd2(:,:,jl), ndimtrd1, ndextrd1 ) 654 END DO 655 656 CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)), & 657 & it, ztmlatf2(:,:), ndimtrd1, ndextrd1 ) 658 659 !.................................. ( ML salinity ) ..................................... 660 661 !-- Output the fields 662 CALL histwrite( nidtrd, "sml" , it, sml_sum (:,:), ndimtrd1, ndextrd1 ) 663 CALL histwrite( nidtrd, "sml_tot" , it, zsmltot2(:,:), ndimtrd1, ndextrd1 ) 664 CALL histwrite( nidtrd, "sml_res" , it, zsmlres2(:,:), ndimtrd1, ndextrd1 ) 665 666 DO jl = 1, jpltrd - 1 667 CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)), & 668 & it, zsmltrd2(:,:,jl), ndimtrd1, ndextrd1 ) 669 END DO 670 671 CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)), & 672 & it, zsmlatf2(:,:), ndimtrd1, ndextrd1 ) 673 674 IF( kt == nitend ) CALL histclo( nidtrd ) 675 676 END IF 677 678 ! Compute the control surface (for next time step) : flag = on 679 icount = 1 680 ! 612 681 #endif 613 682 614 IF( idebug /= 0 ) THEN 615 WRITE(numout,*) ' debuging trd_mld: III.5 done' 616 CALL FLUSH(numout) 617 ENDIF 618 619 ! set counter icount to one to allow the calculation 620 ! of the surface control in the next time step in the trd_mld_zint subroutine 621 icount = 1 622 623 ENDIF 624 625 ! At the end of the 1st time step, set icount to 1 to be 626 ! able to compute the surface control at the beginning of 627 ! the second time step 628 IF( kt == nit000 ) icount = 1 629 630 IF( kt == nitend ) CALL histclo( nidtrd ) 631 #endif 683 IF( MOD( kt , ntrd ) == 0 ) THEN 684 ! 685 ! III.5 Reset cumulative arrays to zero 686 ! ------------------------------------- 687 nmoymltrd = 0 688 689 ! ... temperature ... ... salinity ... 690 tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0 691 tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0 692 tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0 693 tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0 694 tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0 695 696 rmld_sum (:,:) = 0.e0 697 ! 698 END IF 632 699 633 700 END SUBROUTINE trd_mld … … 642 709 !! from ocean surface down to control surface (NetCDF output) 643 710 !! 644 !! ** Method/usage :645 !!646 !! History :647 !! ! 95-04 (J. Vialard) Original code648 !! ! 97-02 (E. Guilyardi) Adaptation global + base cmo649 !! ! 99-09 (E. Guilyardi) Re-writing + netCDF output650 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module651 !! 9.0 ! 04-08 (C. Talandier) New trends organization652 711 !!---------------------------------------------------------------------- 653 712 !! * Local declarations 654 INTEGER :: ilseq 713 INTEGER :: ilseq, jl 655 714 656 715 REAL(wp) :: zjulian, zsto, zout 657 716 658 CHARACTER (LEN=21) :: &717 CHARACTER (LEN=21) :: & 659 718 clold ='OLD' , & ! open specifier (direct access files) 660 719 clunf ='UNFORMATTED', & ! open specifier (direct access files) 661 720 clseq ='SEQUENTIAL' ! open specifier (direct access files) 662 CHARACTER (LEN=40) :: clhstnam663 721 CHARACTER (LEN=40) :: clop 664 CHARACTER (LEN=12) :: clmxl 665 666 NAMELIST/namtrd/ ntrd, nctls 722 CHARACTER (LEN=12) :: clmxl, cltu, clsu 723 667 724 !!---------------------------------------------------------------------- 668 725 669 ! =================== 670 ! I. initialization 671 ! =================== 672 673 ! Open specifier 674 ilseq = 1 675 idebug = 0 ! set it to 1 in case of problem to have more print 676 icount = 1 677 ionce = 1 678 679 ! namelist namtrd : trend diagnostic 680 REWIND( numnam ) 681 READ ( numnam, namtrd ) 726 ! ====================================================================== 727 ! I. initialization 728 ! ====================================================================== 682 729 683 730 IF(lwp) THEN 684 WRITE(numout,*) ' ' 685 WRITE(numout,*) 'trd_mld_init: mixed layer heat & freshwater budget trends' 686 WRITE(numout,*) '~~~~~~~~~~~~~' 687 WRITE(numout,*) ' ' 688 WRITE(numout,*) ' Namelist namtrd : ' 689 WRITE(numout,*) ' control surface for trends nctls = ',nctls 690 WRITE(numout,*) ' ' 691 ENDIF 692 693 ! cumulated trends array init 731 WRITE(numout,*) 732 WRITE(numout,*) ' trd_mld_init : Mixed-layer trends' 733 WRITE(numout,*) ' ~~~~~~~~~~~~~' 734 WRITE(numout,*) ' namelist namtrd read in trd_mod_init ' 735 WRITE(numout,*) 736 END IF 737 738 ! I.1 Check consistency of user defined preferences 739 ! ------------------------------------------------- 740 741 IF( ( lk_trdmld ) .AND. ( MOD( nitend, ntrd ) /= 0 ) ) THEN 742 WRITE(numout,cform_err) 743 WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend 744 WRITE(numout,*) ' is no multiple of the trends diagnostics frequency ' 745 WRITE(numout,*) ' you defined, ntrd = ', ntrd 746 WRITE(numout,*) ' This will not allow you to restart from this simulation. ' 747 WRITE(numout,*) ' You should reconsider this choice. ' 748 WRITE(numout,*) 749 WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a ' 750 WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) ' 751 nstop = nstop + 1 752 END IF 753 754 IF( ( lk_trdmld ) .AND. ( n_cla == 1 ) ) THEN 755 WRITE(numout,cform_war) 756 WRITE(numout,*) ' You set n_cla = 1. Note that the Mixed-Layer diagnostics ' 757 WRITE(numout,*) ' are not exact along the corresponding straits. ' 758 nwarn = nwarn + 1 759 END IF 760 761 ! I.2 Initialize arrays to zero or read a restart file 762 ! ---------------------------------------------------- 763 694 764 nmoymltrd = 0 695 tmltrdm(:,:) = 0.e0 696 smltrdm(:,:) = 0.e0 697 698 ! read control surface from file ctlsurf_idx 699 765 766 ! ... temperature ... ... salinity ... 767 tml (:,:) = 0.e0 ; sml (:,:) = 0.e0 ! inst. 768 tmltrdm (:,:) = 0.e0 ; smltrdm (:,:) = 0.e0 769 tmlatfm (:,:) = 0.e0 ; smlatfm (:,:) = 0.e0 770 tml_sum (:,:) = 0.e0 ; sml_sum (:,:) = 0.e0 ! mean 771 tmltrd_sum (:,:,:) = 0.e0 ; smltrd_sum (:,:,:) = 0.e0 772 tmltrd_csum_ln (:,:,:) = 0.e0 ; smltrd_csum_ln (:,:,:) = 0.e0 773 774 rmld (:,:) = 0.e0 775 rmld_sum (:,:) = 0.e0 776 777 IF( ln_rstart .AND. ln_trdmld_restart ) THEN 778 CALL trd_mld_rst_read 779 ELSE 780 ! ... temperature ... ... salinity ... 781 tmlb (:,:) = 0.e0 ; smlb (:,:) = 0.e0 ! inst. 782 tmlbb (:,:) = 0.e0 ; smlbb (:,:) = 0.e0 783 tmlbn (:,:) = 0.e0 ; smlbn (:,:) = 0.e0 784 tml_sumb (:,:) = 0.e0 ; sml_sumb (:,:) = 0.e0 ! mean 785 tmltrd_csum_ub (:,:,:) = 0.e0 ; smltrd_csum_ub (:,:,:) = 0.e0 786 tmltrd_atf_sumb(:,:) = 0.e0 ; smltrd_atf_sumb(:,:) = 0.e0 787 END IF 788 789 ilseq = 1 ; icount = 1 ; ionce = 1 ! open specifier 790 791 ! I.3 Read control surface from file ctlsurf_idx 792 ! ---------------------------------------------- 793 700 794 IF( nctls == 1 ) THEN 701 clname ='ctlsurf_idx' 702 CALL ctlopn(numbol,clname,clold,clunf,clseq, & 703 ilseq,numout,lwp,1) 704 REWIND (numbol) 705 READ(numbol) nbol 706 ENDIF 707 708 709 IF( idebug /= 0 ) THEN 710 WRITE(numout,*) ' debuging trd_mld_init: 0. done ' 711 CALL FLUSH(numout) 712 ENDIF 713 714 ! =================================== 715 ! II. netCDF output initialization 716 ! =================================== 795 clname = 'ctlsurf_idx' 796 CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 ) 797 REWIND( numbol ) 798 READ ( numbol ) nbol 799 END IF 800 801 ! ====================================================================== 802 ! II. netCDF output initialization 803 ! ====================================================================== 717 804 718 805 #if defined key_dimgout 719 806 ??? 720 807 #else 721 ! clmxl = legend root for netCDF output 722 IF( nctls == 0 ) THEN 723 ! control surface = mixed-layer with density criterion 724 ! (array nmln computed in zdfmxl.F90) 725 clmxl = 'Mixed Layer ' 726 ELSE IF( nctls == 1 ) THEN 727 ! control surface = read index from file 808 ! clmxl = legend root for netCDF output 809 IF( nctls == 0 ) THEN ! control surface = mixed-layer with density criterion 810 clmxl = 'Mixed Layer ' ! (array nmln computed in zdfmxl.F90) 811 ELSE IF( nctls == 1 ) THEN ! control surface = read index from file 728 812 clmxl = ' Bowl ' 729 ELSE IF( nctls >= 2 ) THEN 730 ! control surface = model level 731 WRITE(clmxl,'(A9,I2,1X)') 'Levels 1-', nctls 732 ENDIF 733 734 !----------------------------------------- 813 ELSE IF( nctls >= 2 ) THEN ! control surface = model level 814 WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls 815 END IF 816 735 817 ! II.1 Define frequency of output and means 736 818 ! ----------------------------------------- 737 738 #if defined key_diainstant 739 zsto = nwrite*rdt 740 clop ="inst(x)" 741 #else 742 zsto = rdt 743 clop ="ave(x)" 744 #endif 745 zout = nwrite*rdt 746 747 IF(lwp) WRITE (numout,*) ' trdmld_ncinit: netCDF initialization' 819 # if defined key_diainstant 820 IF( .NOT. ln_trdmld_instant ) THEN 821 CALL ctl_stop( 'trd_mld : this was never checked. Comment this line to proceed...' ) 822 END IF 823 zsto = ntrd * rdt 824 clop ="inst(only(x))" 825 # else 826 IF( ln_trdmld_instant ) THEN 827 zsto = rdt ! inst. diags : we use IOIPSL time averaging 828 ELSE 829 zsto = ntrd * rdt ! mean diags : we DO NOT use any IOIPSL time averaging 830 END IF 831 clop ="ave(only(x))" 832 # endif 833 zout = ntrd * rdt 834 835 IF(lwp) WRITE (numout,*) ' netCDF initialization' 748 836 749 837 ! II.2 Compute julian date from starting date of the run 750 ! ------------------------ 751 838 ! ------------------------------------------------------ 752 839 CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian ) 753 IF 754 IF (lwp) WRITE(numout,*)' Date 0 used :',nit000&755 ,' YEAR ', nyear,' MONTH ', nmonth,' DAY ', nday&756 ,'Julian day : ', zjulian840 IF(lwp) WRITE(numout,*)' ' 841 IF(lwp) WRITE(numout,*)' Date 0 used :',nit000, & 842 & ' YEAR ', nyear,' MONTH ' , nmonth, & 843 & ' DAY ' , nday, 'Julian day : ', zjulian 757 844 758 845 759 846 ! II.3 Define the T grid trend file (nidtrd) 760 ! --------------------------------- 761 762 CALL dia_nam( clhstnam, nwrite, 'trends' ) ! filename 847 ! ------------------------------------------ 848 !-- Define long and short names for the NetCDF output variables 849 ! ==> choose them according to trdmld_oce.F90 <== 850 851 ctrd(jpmld_xad,1) = " Zonal advection" ; ctrd(jpmld_xad,2) = "_xad" 852 ctrd(jpmld_yad,1) = " Meridional advection" ; ctrd(jpmld_yad,2) = "_yad" 853 ctrd(jpmld_zad,1) = " Vertical advection" ; ctrd(jpmld_zad,2) = "_zad" 854 ctrd(jpmld_ldf,1) = " Lateral diffusion" ; ctrd(jpmld_ldf,2) = "_ldf" 855 ctrd(jpmld_for,1) = " Forcing" ; ctrd(jpmld_for,2) = "_for" 856 ctrd(jpmld_zdf,1) = " Vertical diff. (Kz)" ; ctrd(jpmld_zdf,2) = "_zdf" 857 ctrd(jpmld_bbc,1) = " Geothermal flux" ; ctrd(jpmld_bbc,2) = "_bbc" 858 ctrd(jpmld_bbl,1) = " Adv/diff. Bottom boundary layer" ; ctrd(jpmld_bbl,2) = "_bbl" 859 ctrd(jpmld_dmp,1) = " Tracer damping" ; ctrd(jpmld_dmp,2) = "_dmp" 860 ctrd(jpmld_npc,1) = " Non penetrative convec. adjust." ; ctrd(jpmld_npc,2) = "_npc" 861 ctrd(jpmld_atf,1) = " Asselin time filter" ; ctrd(jpmld_atf,2) = "_atf" 862 863 !-- Create a NetCDF file and enter the define mode 864 CALL dia_nam( clhstnam, ntrd, 'trends' ) 763 865 IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam 764 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,1, jpi, & ! Horizontal grid : glamt and gphit 765 & 1, jpj, 0, zjulian, rdt, nh_t, nidtrd, domain_id=nidom ) 766 767 ! Declare output fields as netCDF variables 768 769 ! Mixed layer Depth 770 CALL histdef( nidtrd, "somlttml", clmxl//"Depth" , "m" , & ! hmlp 771 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 772 773 ! Temperature 774 CALL histdef( nidtrd, "somltemp", clmxl//"Temperature" , "C" , & ! ??? 775 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 776 ! Temperature trends 777 CALL histdef( nidtrd, "somlttto", clmxl//"T Total" , "C/s", & ! total 778 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) 779 CALL histdef( nidtrd, "somlttax", clmxl//"T Zonal Advection", "C/s", & ! i-adv. 780 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 781 CALL histdef( nidtrd, "somlttay", clmxl//"T Meridional Advection", "C/s", & ! j-adv. 782 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 783 CALL histdef( nidtrd, "somlttaz", clmxl//"T Vertical Advection", "C/s", & ! vert. adv. 784 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 785 CALL histdef( nidtrd, "somlttdh", clmxl//"T Horizontal Diffusion ", "C/s", & ! hor. lateral diff. 786 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 787 CALL histdef( nidtrd, "somlttfo", clmxl//"T Forcing", "C/s", & ! forcing 788 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 789 CALL histdef( nidtrd, "somlbtdz", clmxl//"T Vertical Diffusion", "C/s", & ! vert. diff. 790 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 791 CALL histdef( nidtrd, "somlbtdt", clmxl//"T dh/dt Entrainment (Residual)", "C/s", & ! T * dh/dt 792 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) 793 IF( ln_traldf_iso ) THEN 794 CALL histdef( nidtrd, "somlbtdv", clmxl//"T Vert. lateral Diffusion","C/s", & ! vertical diffusion entrainment (ISO) 795 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 796 ENDIF 797 #if defined key_traldf_eiv 798 CALL histdef( nidtrd, "somlgtax", clmxl//"T Zonal EIV Advection", "C/s", & ! i-adv. (eiv) 799 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 800 CALL histdef( nidtrd, "somlgtay", clmxl//"T Meridional EIV Advection", "C/s", & ! j-adv. (eiv) 801 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 802 CALL histdef( nidtrd, "somlgtaz", clmxl//"T Vertical EIV Advection", "C/s", & ! vert. adv. (eiv) 803 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 804 CALL histdef( nidtrd, "somlgtat", clmxl//"T Total EIV Advection", "C/s", & ! total advection (eiv) 805 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 806 #endif 807 ! Salinity 808 CALL histdef( nidtrd, "somlsalt", clmxl//"Salinity", "PSU", & ! Mixed-layer salinity 809 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 810 ! Salinity trends 811 CALL histdef( nidtrd, "somltsto", clmxl//"S Total", "PSU/s", & ! total 812 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 813 CALL histdef( nidtrd, "somltsax", clmxl//"S Zonal Advection", "PSU/s", & ! i-advection 814 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 815 CALL histdef( nidtrd, "somltsay", clmxl//"S Meridional Advection", "PSU/s", & ! j-advection 816 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 817 CALL histdef( nidtrd, "somltsaz", clmxl//"S Vertical Advection", "PSU/s", & ! vertical advection 818 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 819 CALL histdef( nidtrd, "somltsdh", clmxl//"S Horizontal Diffusion ", "PSU/s", & ! hor. lat. diff. 820 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 821 CALL histdef( nidtrd, "somltsfo", clmxl//"S Forcing", "PSU/s", & ! forcing 822 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 823 824 CALL histdef( nidtrd, "somlbsdz", clmxl//"S Vertical Diffusion", "PSU/s", & ! vert. diff. 825 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 826 CALL histdef( nidtrd, "somlbsdt", clmxl//"S dh/dt Entrainment (Residual)", "PSU/s", & ! S * dh/dt 827 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 828 IF( ln_traldf_iso ) THEN 829 ! vertical diffusion entrainment (ISO) 830 CALL histdef( nidtrd, "somlbsdv", clmxl//"S Vertical lateral Diffusion", "PSU/s", & ! vert. lat. diff. 831 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 832 ENDIF 833 #if defined key_traldf_eiv 834 CALL histdef( nidtrd, "somlgsax", clmxl//"S Zonal EIV Advection", "PSU/s", & ! i-advection (eiv) 835 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 836 CALL histdef( nidtrd, "somlgsay", clmxl//"S Meridional EIV Advection", "PSU/s", & ! j-advection (eiv) 837 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 838 CALL histdef( nidtrd, "somlgsaz", clmxl//"S Vertical EIV Advection", "PSU/s", & ! vert. adv. (eiv) 839 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 840 CALL histdef( nidtrd, "somlgsat", clmxl//"S Total EIV Advection", "PSU/s", & ! total adv. (eiv) 841 & jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 842 #endif 866 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & 867 & 1, jpi, 1, jpj, 0, zjulian, rdt, nh_t, nidtrd, domain_id=nidom ) 868 869 !-- Define the ML depth variable 870 CALL histdef(nidtrd, "mxl_depth", clmxl//" Mixed Layer Depth" , "m", & 871 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 872 873 !-- Define physical units 874 IF( ucf == 1. ) THEN 875 cltu = "degC/s" ; clsu = "p.s.u./s" 876 ELSEIF ( ucf == 3600.*24.) THEN 877 cltu = "degC/day" ; clsu = "p.s.u./day" 878 ELSE 879 cltu = "unknown?" ; clsu = "unknown?" 880 END IF 881 882 !-- Define miscellaneous T and S mixed-layer variables 883 884 IF( jpltrd /= jpmld_atf ) CALL ctl_stop( 'Error : jpltrd /= jpmld_atf' ) ! see below 885 886 !................................. ( ML temperature ) ................................... 887 888 CALL histdef(nidtrd, "tml" , clmxl//" T Mixed Layer Temperature" , "C", & 889 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 890 CALL histdef(nidtrd, "tml_tot", clmxl//" T Total trend" , cltu, & 891 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) 892 CALL histdef(nidtrd, "tml_res", clmxl//" T dh/dt Entrainment (Resid.)" , cltu, & 893 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) 894 895 DO jl = 1, jpltrd - 1 ! <== only true if jpltrd == jpmld_atf 896 CALL histdef(nidtrd, trim("tml"//ctrd(jl,2)), clmxl//" T"//ctrd(jl,1), cltu, & 897 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean 898 END DO ! if zsto=rdt above 899 900 CALL histdef(nidtrd, trim("tml"//ctrd(jpmld_atf,2)), clmxl//" T"//ctrd(jpmld_atf,1), cltu, & 901 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean 902 903 !.................................. ( ML salinity ) ..................................... 904 905 CALL histdef(nidtrd, "sml" , clmxl//" S Mixed Layer Salinity" , "p.s.u.", & 906 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 907 CALL histdef(nidtrd, "sml_tot", clmxl//" S Total trend" , clsu, & 908 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) 909 CALL histdef(nidtrd, "sml_res", clmxl//" S dh/dt Entrainment (Resid.)" , clsu, & 910 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) 911 912 DO jl = 1, jpltrd - 1 ! <== only true if jpltrd == jpmld_atf 913 CALL histdef(nidtrd, trim("sml"//ctrd(jl,2)), clmxl//" S"//ctrd(jl,1), clsu, & 914 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean 915 END DO ! if zsto=rdt above 916 917 CALL histdef(nidtrd, trim("sml"//ctrd(jpmld_atf,2)), clmxl//" S"//ctrd(jpmld_atf,1), clsu, & 918 jpi, jpj, nh_t, 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean 919 920 !-- Leave IOIPSL/NetCDF define mode 843 921 CALL histend( nidtrd ) 844 #endif 845 846 IF( idebug /= 0 ) THEN 847 WRITE(numout,*) ' debuging trd_mld_init: II. done' 848 CALL FLUSH(numout) 849 ENDIF 850 851 852 END SUBROUTINE trd_mld_init 922 923 #endif /* key_dimgout */ 924 END SUBROUTINE trd_mld_init 853 925 854 926 #else … … 856 928 !! Default option : Empty module 857 929 !!---------------------------------------------------------------------- 858 LOGICAL, PUBLIC, PARAMETER :: lk_trdmld = .FALSE. !: momentum trend flag859 930 CONTAINS 860 931 SUBROUTINE trd_mld( kt ) ! Empty routine
Note: See TracChangeset
for help on using the changeset viewer.
