- Timestamp:
- 2017-10-16T16:19:11+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r4990 r8627 18 18 USE oce ! ocean dynamics and tracers 19 19 USE dom_oce ! ocean space and time domain 20 USE phycst ! physical constants 20 21 USE ldfdyn_oce ! ocean dynamics: lateral physics 21 22 ! 22 23 USE in_out_manager ! I/O manager 24 USE iom ! I/O library 23 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 26 USE wrk_nemo ! Memory Allocation … … 75 77 INTEGER, INTENT(in) :: kt ! ocean time-step index 76 78 ! 77 INTEGER :: ji, jj, jk ! dummy loop indices 78 REAL(wp) :: zua, zva, zbt, ze2u, ze2v ! temporary scalar 79 REAL(wp), POINTER, DIMENSION(:,: ) :: zcu, zcv 80 REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv 79 INTEGER :: ji, jj, jk ! dummy loop indices 80 REAL(wp) :: zua, zva, zbt, ze2u, ze2v, zzz ! local scalar 81 REAL(wp), POINTER, DIMENSION(:,: ) :: zcu, zcv 82 REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace 81 84 !!---------------------------------------------------------------------- 82 85 ! … … 112 115 DO jj = 2, jpjm1 113 116 DO ji = fs_2, fs_jpim1 ! vector opt. 114 zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) -zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) &115 & + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) /e1u(ji,jj)117 zlu(ji,jj,jk) = - ( zuf(ji ,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) & 118 & + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj ,jk) ) / e1u(ji,jj) 116 119 117 zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) -zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) &118 & + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) /e2v(ji,jj)120 zlv(ji,jj,jk) = + ( zuf(ji,jj ,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) & 121 & + ( hdivb(ji,jj+1,jk) - hdivb(ji ,jj,jk) ) / e2v(ji,jj) 119 122 END DO 120 123 END DO … … 123 126 DO ji = fs_2, fs_jpim1 ! vector opt. 124 127 zlu(ji,jj,jk) = - ( rotb (ji ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj) & 125 & + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj ,jk) ) / e1u(ji,jj)128 & + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj ,jk) ) / e1u(ji,jj) 126 129 127 130 zlv(ji,jj,jk) = + ( rotb (ji,jj ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj) & 128 & + ( hdivb(ji,jj+1,jk) - hdivb(ji ,jj,jk) ) / e2v(ji,jj)131 & + ( hdivb(ji,jj+1,jk) - hdivb(ji ,jj,jk) ) / e2v(ji,jj) 129 132 END DO 130 133 END DO … … 133 136 CALL lbc_lnk( zlu, 'U', -1. ) ; CALL lbc_lnk( zlv, 'V', -1. ) ! Boundary conditions 134 137 135 138 IF( iom_use('dispkexyfo') ) THEN ! ocean kinetic energy dissipation per unit area 139 ! ! due to xy friction (xy=lateral) 140 ! see NEMO_book appendix C, §C.7.2 (N.B. here averaged at t-points) 141 ! local dissipation of KE at t-point due to bilaplacian operator is given by : 142 ! + ahmu mi( zlu**2 ) + mj( ahmv zlv**2 ) 143 ! Note that a sign + is used as in v3.6 ahm is negative for bilaplacian viscosity 144 ! 145 ! NB: ahm is negative when bilaplacian is used 146 ALLOCATE( z2d(jpi,jpj) ) 147 z2d(:,:) = 0._wp 148 DO jk = 1, jpkm1 149 DO jj = 2, jpjm1 150 DO ji = 2, jpim1 151 z2d(ji,jj) = z2d(ji,jj) & 152 & + ( fsahmu(ji,jj,jk)*zlu(ji,jj,jk)**2 + fsahmu(ji-1,jj,jk)*zlu(ji-1,jj,jk)**2 & 153 & + fsahmv(ji,jj,jk)*zlv(ji,jj,jk)**2 + fsahmv(ji,jj-1,jk)*zlv(ji,jj-1,jk)**2 ) * tmask(ji,jj,jk) 154 END DO 155 END DO 156 END DO 157 zzz = 0.5_wp * rau0 158 z2d(:,:) = zzz * z2d(:,:) 159 CALL lbc_lnk( z2d,'T', 1. ) 160 CALL iom_put( 'dispkexyfo', z2d ) 161 DEALLOCATE( z2d ) 162 ENDIF 163 164 165 ! Third derivative 166 ! ---------------- 167 ! 136 168 DO jk = 1, jpkm1 137 138 ! Third derivative 139 ! ---------------- 140 169 ! 141 170 ! Multiply by the eddy viscosity coef. (at u- and v-points) 142 zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 143 144 zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 145 171 zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk) 172 zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk) 173 ! 146 174 ! Contravariant "laplacian" 147 175 zcu(:,:) = e1u(:,:) * zlu(:,:,jk) … … 152 180 DO ji = 1, fs_jpim1 ! vector opt. 153 181 zuf(ji,jj,jk) = fmask(ji,jj,jk) * ( zcv(ji+1,jj ) - zcv(ji,jj) & 154 & - zcu(ji ,jj+1) + zcu(ji,jj) ) &155 & * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj))182 & - zcu(ji ,jj+1) + zcu(ji,jj) ) & 183 & * fse3f(ji,jj,jk) * r1_e12f(ji,jj) 156 184 END DO 157 185 END DO … … 175 203 END DO 176 204 177 178 205 ! boundary conditions on the laplacian curl and div (zuf,zut) 179 206 !!bug gm no need to do this 2 following lbc... … … 181 208 CALL lbc_lnk( zut, 'T', 1. ) 182 209 183 DO jk = 1, jpkm1 184 185 ! Bilaplacian 186 ! ----------- 187 210 DO jk = 1, jpkm1 ! Bilaplacian 188 211 DO jj = 2, jpjm1 189 212 DO ji = fs_2, fs_jpim1 ! vector opt. … … 197 220 & + ( zut(ji,jj+1,jk) - zut(ji ,jj,jk) ) / e2v(ji,jj) 198 221 ! add it to the general momentum trends 199 ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 200 va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 201 END DO 202 END DO 203 222 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 223 va(ji,jj,jk) = va(ji,jj,jk) + zva 224 END DO 225 END DO 204 226 ! ! =============== 205 227 END DO ! End of slab
Note: See TracChangeset
for help on using the changeset viewer.
