- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r5930 r7351 12 12 13 13 !!---------------------------------------------------------------------- 14 !! dyn_zdf_exp : update the momentum trend with the vertical diffu-15 !! 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. 16 16 !!---------------------------------------------------------------------- 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 USE lib_mpp ! MPP library 24 USE in_out_manager ! I/O manager 25 USE lib_mpp ! MPP library 26 USE wrk_nemo ! Memory Allocation 27 USE timing ! Timing 28 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 29 28 30 29 IMPLICIT NONE … … 34 33 35 34 !! * Substitutions 36 # include "domzgr_substitute.h90"37 35 # include "vectopt_loop_substitute.h90" 38 36 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)37 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 40 38 !! $Id$ 41 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 48 46 !! 49 47 !! ** Purpose : Compute the trend due to the vert. momentum diffusion 48 !! and perform the Leap-Frog time stepping. 50 49 !! 51 !! ** Method : Explicit forward time stepping with a time splitting52 !! technique. The vertical diffusionof momentum is given by:50 !! ** Method : - Split-explicit forward time stepping. 51 !! The vertical mixing of momentum is given by: 53 52 !! diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) ) 54 53 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) … … 56 55 !! Add this trend to the general trend ua : 57 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 58 60 !! 59 !! ** Action : - Update (ua,va) with the vertical diffusive trend61 !! ** Action : - (ua,va) after velocity 60 62 !!--------------------------------------------------------------------- 61 63 INTEGER , INTENT(in) :: kt ! ocean time-step index 62 64 REAL(wp), INTENT(in) :: p2dt ! time-step 63 65 ! 64 INTEGER :: ji, jj, jk, jl ! dummy loop indices66 INTEGER :: ji, jj, jk, jl ! dummy loop indices 65 67 REAL(wp) :: zlavmr, zua, zva ! local scalars 66 68 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwz, zww 67 69 !!---------------------------------------------------------------------- 68 70 ! 69 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_exp')71 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_exp') 70 72 ! 71 CALL wrk_alloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )73 CALL wrk_alloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) 72 74 ! 73 75 IF( kt == nit000 .AND. lwp ) THEN … … 76 78 WRITE(numout,*) '~~~~~~~~~~~ ' 77 79 ENDIF 78 80 ! 81 ! !== vertical mixing trend ==! 82 ! 79 83 zlavmr = 1. / REAL( nn_zdfexp ) 80 81 82 DO jj = 2, jpjm1 ! Surface boundary condition 84 ! 85 DO jj = 2, jpjm1 ! Surface boundary condition 83 86 DO ji = 2, jpim1 84 87 zwy(ji,jj,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * r1_rau0 … … 86 89 END DO 87 90 END DO 88 DO jk = 1, jpk 91 DO jk = 1, jpk ! Initialization of x, z and contingently trends array 89 92 DO jj = 2, jpjm1 90 93 DO ji = 2, jpim1 … … 95 98 END DO 96 99 ! 97 DO jl = 1, nn_zdfexp 100 DO jl = 1, nn_zdfexp ! Time splitting loop 98 101 ! 99 DO jk = 2, jpk 102 DO jk = 2, jpk ! First vertical derivative 100 103 DO jj = 2, jpjm1 101 104 DO ji = 2, jpim1 102 zwy(ji,jj,jk) = avmu(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) / fse3uw(ji,jj,jk)103 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) 104 107 END DO 105 108 END DO 106 109 END DO 107 DO jk = 1, jpkm1 110 DO jk = 1, jpkm1 ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 108 111 DO jj = 2, jpjm1 109 112 DO ji = 2, jpim1 110 zua = zlavmr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) / fse3u(ji,jj,jk)111 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) 112 115 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 113 116 va(ji,jj,jk) = va(ji,jj,jk) + zva … … 118 121 END DO 119 122 END DO 120 ! 121 END DO ! End of time splitting 122 123 ! Time step momentum here to be compliant with what is done in the implicit case 123 END DO ! End of time splitting 124 124 ! 125 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 125 ! 126 ! !== Leap-Frog time integration ==! 127 ! 128 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 126 129 DO jk = 1, jpkm1 127 130 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 128 131 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 129 132 END DO 130 ELSE 133 ELSE ! applied on thickness weighted velocity 131 134 DO jk = 1, jpkm1 132 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 133 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 134 & / fse3u_a(:,:,jk) * umask(:,:,jk) 135 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 136 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 137 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 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) 138 139 END DO 139 140 ENDIF 140 141 ! 141 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )142 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) 142 143 ! 143 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_exp')144 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_exp') 144 145 ! 145 146 END SUBROUTINE dyn_zdf_exp
Note: See TracChangeset
for help on using the changeset viewer.
