- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r3625 r6808 8 8 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 !! 3.7 ! 2015-11 (J. Chanut) output velocities instead of trends 10 11 !!---------------------------------------------------------------------- 11 12 12 13 !!---------------------------------------------------------------------- 13 !! dyn_zdf_exp : update the momentum trend with the vertical diffu-14 !! sion using an explicit time-stepping scheme.14 !! dyn_zdf_exp : update the momentum trend with the vertical diffusion using a split-explicit scheme 15 !! and perform the Leap-Frog time integration. 15 16 !!---------------------------------------------------------------------- 16 USE oce 17 USE dom_oce 18 USE phycst 19 USE zdf_oce 20 USE sbc_oce ! surface boundary condition: ocean21 USE lib_mpp ! MPP library22 USE in_out_manager ! I/O manager23 USE lib_mpp ! MPP library24 USE wrk_nemo ! Memory Allocation25 USE timing ! Timing26 17 USE oce ! ocean dynamics and tracers 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 20 USE zdf_oce ! ocean vertical physics 21 USE dynadv , ONLY: ln_dynadv_vec ! Momentum advection form 22 USE sbc_oce ! surface boundary condition: ocean 23 ! 24 USE in_out_manager ! I/O manager 25 USE lib_mpp ! MPP library 26 USE wrk_nemo ! Memory Allocation 27 USE timing ! Timing 27 28 28 29 IMPLICIT NONE … … 32 33 33 34 !! * Substitutions 34 # include "domzgr_substitute.h90"35 35 # include "vectopt_loop_substitute.h90" 36 36 !!---------------------------------------------------------------------- 37 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)37 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 38 38 !! $Id$ 39 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 46 46 !! 47 47 !! ** Purpose : Compute the trend due to the vert. momentum diffusion 48 !! and perform the Leap-Frog time stepping. 48 49 !! 49 !! ** Method : Explicit forward time stepping with a time splitting50 !! technique. The vertical diffusionof momentum is given by:50 !! ** Method : - Split-explicit forward time stepping. 51 !! The vertical mixing of momentum is given by: 51 52 !! diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) ) 52 53 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) … … 54 55 !! Add this trend to the general trend ua : 55 56 !! ua = ua + dz( avmu dz(u) ) 57 !! - Leap-Frog time stepping (Asselin filter will be applied in dyn_nxt) 58 !! ua = ub + 2*dt * ua vector form or linear free surf. 59 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 56 60 !! 57 !! ** Action : - Update (ua,va) with the vertical diffusive trend61 !! ** Action : - (ua,va) after velocity 58 62 !!--------------------------------------------------------------------- 59 63 INTEGER , INTENT(in) :: kt ! ocean time-step index 60 64 REAL(wp), INTENT(in) :: p2dt ! time-step 61 65 ! 62 INTEGER :: ji, jj, jk, jl ! dummy loop indices66 INTEGER :: ji, jj, jk, jl ! dummy loop indices 63 67 REAL(wp) :: zlavmr, zua, zva ! local scalars 64 68 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwz, zww 65 69 !!---------------------------------------------------------------------- 66 70 ! 67 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_exp')71 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_exp') 68 72 ! 69 CALL wrk_alloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )73 CALL wrk_alloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) 70 74 ! 71 75 IF( kt == nit000 .AND. lwp ) THEN … … 74 78 WRITE(numout,*) '~~~~~~~~~~~ ' 75 79 ENDIF 76 80 ! 81 ! !== vertical mixing trend ==! 82 ! 77 83 zlavmr = 1. / REAL( nn_zdfexp ) 78 79 80 DO jj = 2, jpjm1 ! Surface boundary condition 84 ! 85 DO jj = 2, jpjm1 ! Surface boundary condition 81 86 DO ji = 2, jpim1 82 87 zwy(ji,jj,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * r1_rau0 … … 84 89 END DO 85 90 END DO 86 DO jk = 1, jpk 91 DO jk = 1, jpk ! Initialization of x, z and contingently trends array 87 92 DO jj = 2, jpjm1 88 93 DO ji = 2, jpim1 … … 93 98 END DO 94 99 ! 95 DO jl = 1, nn_zdfexp 100 DO jl = 1, nn_zdfexp ! Time splitting loop 96 101 ! 97 DO jk = 2, jpk 102 DO jk = 2, jpk ! First vertical derivative 98 103 DO jj = 2, jpjm1 99 104 DO ji = 2, jpim1 100 zwy(ji,jj,jk) = avmu(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) / fse3uw(ji,jj,jk)101 zww(ji,jj,jk) = avmv(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) / fse3vw(ji,jj,jk)105 zwy(ji,jj,jk) = avmu(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) / e3uw_n(ji,jj,jk) 106 zww(ji,jj,jk) = avmv(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) / e3vw_n(ji,jj,jk) 102 107 END DO 103 108 END DO 104 109 END DO 105 DO jk = 1, jpkm1 110 DO jk = 1, jpkm1 ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 106 111 DO jj = 2, jpjm1 107 112 DO ji = 2, jpim1 108 zua = zlavmr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) / fse3u(ji,jj,jk)109 zva = zlavmr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) / fse3v(ji,jj,jk)113 zua = zlavmr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) / e3u_n(ji,jj,jk) 114 zva = zlavmr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) / e3v_n(ji,jj,jk) 110 115 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 111 116 va(ji,jj,jk) = va(ji,jj,jk) + zva … … 116 121 END DO 117 122 END DO 118 ! 119 END DO ! End of time splitting 123 END DO ! End of time splitting 120 124 ! 121 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )122 125 ! 123 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_exp') 126 ! !== Leap-Frog time integration ==! 127 ! 128 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 129 DO jk = 1, jpkm1 130 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 131 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 132 END DO 133 ELSE ! applied on thickness weighted velocity 134 DO jk = 1, jpkm1 135 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 136 & + p2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 137 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 138 & + p2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 139 END DO 140 ENDIF 141 ! 142 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) 143 ! 144 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_exp') 124 145 ! 125 146 END SUBROUTINE dyn_zdf_exp
Note: See TracChangeset
for help on using the changeset viewer.
