- 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/dynadv_cen2.F90
r4990 r6808 10 10 11 11 !!---------------------------------------------------------------------- 12 !! dyn_adv_cen2 : flux form momentum advection (ln_dynadv_cen2=T) 13 !! trends using a 2nd order centred scheme 12 !! dyn_adv_cen2 : flux form momentum advection (ln_dynadv_cen2=T) using a 2nd order centred scheme 14 13 !!---------------------------------------------------------------------- 15 14 USE oce ! ocean dynamics and tracers … … 30 29 31 30 !! * Substitutions 32 # include "domzgr_substitute.h90"33 31 # include "vectopt_loop_substitute.h90" 34 32 !!---------------------------------------------------------------------- … … 53 51 ! 54 52 INTEGER :: ji, jj, jk ! dummy loop indices 55 REAL(wp) :: zbu, zbv ! local scalars56 53 REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw 57 54 REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfv … … 60 57 IF( nn_timing == 1 ) CALL timing_start('dyn_adv_cen2') 61 58 ! 62 CALL wrk_alloc( jpi, jpj, jpk,zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )59 CALL wrk_alloc( jpi,jpj,jpk, zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 63 60 ! 64 61 IF( kt == nit000 .AND. lwp ) THEN … … 68 65 ENDIF 69 66 ! 70 IF( l_trddyn ) THEN ! Save ua and vatrends67 IF( l_trddyn ) THEN ! trends: store the input trends 71 68 zfu_uw(:,:,:) = ua(:,:,:) 72 69 zfv_vw(:,:,:) = va(:,:,:) 73 70 ENDIF 74 75 ! ! ====================== ! 76 ! ! Horizontal advection ! 77 DO jk = 1, jpkm1 ! ====================== ! 78 ! ! horizontal volume fluxes 79 zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 80 zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 81 ! 82 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point 71 ! 72 ! !== Horizontal advection ==! 73 ! 74 DO jk = 1, jpkm1 ! horizontal transport 75 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 76 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 77 DO jj = 1, jpjm1 ! horizontal momentum fluxes (at T- and F-point) 83 78 DO ji = 1, fs_jpim1 ! vector opt. 84 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj 85 zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj 86 zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji 87 zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji 79 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) 80 zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji ,jj+1,jk) ) 81 zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) ) 82 zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 88 83 END DO 89 84 END DO 90 DO jj = 2, jpjm1 85 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 91 86 DO ji = fs_2, fs_jpim1 ! vector opt. 92 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 93 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 94 ! 95 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 96 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 97 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 98 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 87 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 88 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 89 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 90 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 99 91 END DO 100 92 END DO 101 93 END DO 102 94 ! 103 IF( l_trddyn ) THEN ! save the horizontal advection trendfor diagnostic95 IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic 104 96 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 105 97 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) … … 109 101 ENDIF 110 102 ! 111 112 ! ! ==================== ! 113 ! ! Vertical advection ! 114 DO jk = 1, jpkm1 ! ==================== ! 115 ! ! Vertical volume fluxesÊ 116 zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk) 117 ! 118 IF( jk == 1 ) THEN ! surface/bottom advective fluxes 119 zfu_uw(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 120 zfv_vw(:,:,jpk) = 0.e0 121 ! ! Surface value : 122 IF( lk_vvl ) THEN ! variable volume : flux set to zero 123 zfu_uw(:,:, 1 ) = 0.e0 124 zfv_vw(:,:, 1 ) = 0.e0 125 ELSE ! constant volume : advection through the surface 126 DO jj = 2, jpjm1 127 DO ji = fs_2, fs_jpim1 128 zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj ,1) ) * un(ji,jj,1) 129 zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji ,jj+1,1) ) * vn(ji,jj,1) 130 END DO 131 END DO 132 ENDIF 133 ELSE ! interior fluxes 134 DO jj = 2, jpjm1 135 DO ji = fs_2, fs_jpim1 ! vector opt. 136 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 137 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 138 END DO 103 ! !== Vertical advection ==! 104 ! 105 DO jj = 2, jpjm1 ! surface/bottom advective fluxes set to zero 106 DO ji = fs_2, fs_jpim1 107 zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(jj,jj,jpk) = 0._wp 108 zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(jj,jj, 1 ) = 0._wp 109 END DO 110 END DO 111 IF( ln_linssh ) THEN ! linear free surface: advection through the surface 112 DO jj = 2, jpjm1 113 DO ji = fs_2, fs_jpim1 114 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 115 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 139 116 END DO 140 ENDIF 117 END DO 118 ENDIF 119 DO jk = 2, jpkm1 ! interior advective fluxes 120 DO jj = 2, jpjm1 ! 1/4 * Vertical transport 121 DO ji = fs_2, fs_jpim1 122 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 123 END DO 124 END DO 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 ! vector opt. 127 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 128 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 129 END DO 130 END DO 141 131 END DO 142 DO jk = 1, jpkm1 132 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence 143 133 DO jj = 2, jpjm1 144 134 DO ji = fs_2, fs_jpim1 ! vector opt. 145 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 146 & / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 147 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 148 & / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 135 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 136 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 149 137 END DO 150 138 END DO 151 139 END DO 152 140 ! 153 IF( l_trddyn ) THEN ! save the vertical advection trendfor diagnostic141 IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic 154 142 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 155 143 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 156 144 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 157 145 ENDIF 158 ! 146 ! ! Control print 159 147 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask, & 160 148 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
Note: See TracChangeset
for help on using the changeset viewer.
