Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
r5836 r6140 20 20 USE oce ! ocean dynamics and tracers 21 21 USE dom_oce ! ocean space and time domain 22 USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean22 USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean 23 23 USE sbcrnf ! river runoff 24 24 USE sbcisf ! ice shelf 25 USE iscplhsb ! ice sheet / ocean coupling 26 USE iscplini ! ice sheet / ocean coupling 25 27 ! 26 28 USE in_out_manager ! I/O manager … … 36 38 37 39 !! * Substitutions 38 # include "domzgr_substitute.h90"39 40 # include "vectopt_loop_substitute.h90" 40 41 !!---------------------------------------------------------------------- … … 74 75 DO jj = 2, jpjm1 75 76 DO ji = fs_2, fs_jpim1 ! vector opt. 76 hdivn(ji,jj,jk) = ( e2u(ji ,jj) * fse3u_n(ji ,jj,jk) * un(ji ,jj,jk) &77 & - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * un(ji-1,jj,jk) &78 & + e1v(ji,jj ) * fse3v_n(ji,jj ,jk) * vn(ji,jj ,jk) &79 & - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk) ) &80 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )77 hdivn(ji,jj,jk) = ( e2u(ji ,jj) * e3u_n(ji ,jj,jk) * un(ji ,jj,jk) & 78 & - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * un(ji-1,jj,jk) & 79 & + e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vn(ji,jj ,jk) & 80 & - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk) ) & 81 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 81 82 END DO 82 83 END DO … … 89 90 END DO 90 91 ! 91 IF( ln_rnf 92 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) !== runoffs ==! (update hdivn field) 92 93 ! 93 IF( ln_ divisf .AND. nn_isf > 0) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field)94 IF( ln_isf ) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field) 94 95 ! 95 CALL lbc_lnk( hdivn, 'T', 1. ) !== lateral boundary cond. ==! (no sign change) 96 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn ) !== ice sheet ==! (update hdivn field) 97 ! 98 CALL lbc_lnk( hdivn, 'T', 1. ) !== lateral boundary cond. ==! (no sign change) 96 99 ! 97 100 IF( nn_timing == 1 ) CALL timing_stop('div_hor') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r5836 r6140 39 39 40 40 !! * Substitutions 41 # include "domzgr_substitute.h90"42 41 # include "vectopt_loop_substitute.h90" 43 42 !!---------------------------------------------------------------------- … … 98 97 READ ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) 99 98 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_adv in reference namelist', lwp ) 100 99 ! 101 100 REWIND( numnam_cfg ) ! Namelist namdyn_adv in configuration namelist : Momentum advection scheme 102 101 READ ( numnam_cfg, namdyn_adv, IOSTAT = ios, ERR = 902 ) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r5836 r6140 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 = e1e2u(ji,jj) * fse3u(ji,jj,jk) 93 zbv = e1e2v(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 * e1e2t(:,:) * 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 & / ( e1e2u(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 & / ( e1e2v(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' ) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r5836 r6140 35 35 36 36 !! * Substitutions 37 # include "domzgr_substitute.h90"38 37 # include "vectopt_loop_substitute.h90" 39 38 !!---------------------------------------------------------------------- … … 71 70 !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 72 71 !!---------------------------------------------------------------------- 73 INTEGER, INTENT(in) :: kt ! ocean time-step index 74 ! 75 INTEGER :: ji, jj, jk ! dummy loop indices 76 REAL(wp) :: zbu, zbv ! temporary scalars 77 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! temporary scalars 72 INTEGER, INTENT(in) :: kt ! ocean time-step index 73 ! 74 INTEGER :: ji, jj, jk ! dummy loop indices 75 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars 78 76 REAL(wp), POINTER, DIMENSION(:,:,: ) :: zfu, zfv 79 77 REAL(wp), POINTER, DIMENSION(:,:,: ) :: zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw … … 83 81 IF( nn_timing == 1 ) CALL timing_start('dyn_adv_ubs') 84 82 ! 85 CALL wrk_alloc( jpi, jpj, jpk,zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )86 CALL wrk_alloc( jpi, jpj, jpk, jpts,zlu_uu, zlv_vv, zlu_uv, zlv_vu )83 CALL wrk_alloc( jpi,jpj,jpk, zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 84 CALL wrk_alloc( jpi,jpj,jpk,jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu ) 87 85 ! 88 86 IF( kt == nit000 ) THEN … … 101 99 zlu_uv(:,:,:,:) = 0._wp 102 100 zlv_vu(:,:,:,:) = 0._wp 103 104 IF( l_trddyn ) THEN ! Save ua and vatrends101 ! 102 IF( l_trddyn ) THEN ! trends: store the input trends 105 103 zfu_uw(:,:,:) = ua(:,:,:) 106 104 zfv_vw(:,:,:) = va(:,:,:) 107 105 ENDIF 108 109 106 ! ! =========================== ! 110 107 DO jk = 1, jpkm1 ! Laplacian of the velocity ! 111 108 ! ! =========================== ! 112 109 ! ! horizontal volume fluxes 113 zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)114 zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)110 zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 111 zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 115 112 ! 116 113 DO jj = 2, jpjm1 ! laplacian 117 114 DO ji = fs_2, fs_jpim1 ! vector opt. 118 !119 115 zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj ,jk) ) * umask(ji,jj,jk) 120 116 zlv_vv(ji,jj,jk,1) = ( vb (ji ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji ,jj-1,jk) ) * vmask(ji,jj,jk) … … 137 133 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 138 134 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 139 135 ! 140 136 ! ! ====================== ! 141 137 ! ! Horizontal advection ! 142 138 DO jk = 1, jpkm1 ! ====================== ! 143 139 ! ! horizontal volume fluxes 144 zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)145 zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)140 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 141 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 146 142 ! 147 143 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point … … 150 146 zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 151 147 ! 152 IF (zui > 0) THEN ; zl_u = zlu_uu(ji ,jj,jk,1)153 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1)154 ENDIF 155 IF (zvj > 0) THEN ; zl_v = zlv_vv(ji,jj ,jk,1)156 ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1)148 IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) 149 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) 150 ENDIF 151 IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,jk,1) 152 ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1) 157 153 ENDIF 158 154 ! … … 166 162 zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) ) 167 163 zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) ) 168 IF (zfuj > 0) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1)169 ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1)170 ENDIF 171 IF (zfvi > 0) THEN ; zl_u = zlu_uv( ji,jj ,jk,1)172 ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1)164 IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1) 165 ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1) 166 ENDIF 167 IF( zfvi > 0 ) THEN ; zl_u = zlu_uv( ji,jj ,jk,1) 168 ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1) 173 169 ENDIF 174 170 ! … … 181 177 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 182 178 DO ji = fs_2, fs_jpim1 ! vector opt. 183 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 184 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 185 ! 186 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 187 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 188 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 189 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 190 END DO 191 END DO 192 END DO 193 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 179 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 180 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 181 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 182 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 183 END DO 184 END DO 185 END DO 186 IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic 194 187 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 195 188 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) … … 198 191 zfv_t(:,:,:) = va(:,:,:) 199 192 ENDIF 200 201 193 ! ! ==================== ! 202 194 ! ! Vertical advection ! 203 DO jk = 1, jpkm1! ==================== !204 ! ! Vertical volume fluxesÊ205 zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk)206 !207 IF( jk == 1 ) THEN ! surface/bottom advective fluxes208 zfu_uw( :,:,jpk) = 0.e0 ! Bottom value : flux set to zero209 zfv_vw( :,:,jpk) = 0.e0210 ! ! Surface value :211 IF( lk_vvl ) THEN ! variable volume : flux set to zero212 zfu_uw(:,:, 1 ) = 0.e0213 zfv_vw(:,:, 1 ) = 0.e0214 ELSE ! constant volume : advection through the surface215 DO jj = 2, jpjm1216 DO ji = fs_2, fs_jpim1217 zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj ,1) ) * un(ji,jj,1)218 zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji ,jj+1,1) ) * vn(ji,jj,1)219 END DO220 END DO221 ENDIF222 ELSE ! interior fluxes223 DO jj = 2, jpjm1224 DO ji = fs_2, fs_jpim1 ! vector opt.225 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )226 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )227 END DO228 END DO229 ENDIF230 END DO231 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence232 DO jj = 2, jpjm1233 DO ji = fs_2, fs_jpim1 ! vector opt.234 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) &235 & / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )236 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) &237 & / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))238 END DO 239 END DO 240 END DO 241 ! 242 IF( l_trddyn ) THEN 195 ! ! ==================== ! 196 DO jj = 2, jpjm1 ! surface/bottom advective fluxes set to zero 197 DO ji = fs_2, fs_jpim1 198 zfu_uw(ji,jj,jpk) = 0._wp 199 zfv_vw(ji,jj,jpk) = 0._wp 200 zfu_uw(ji,jj, 1 ) = 0._wp 201 zfv_vw(ji,jj, 1 ) = 0._wp 202 END DO 203 END DO 204 IF( ln_linssh ) THEN ! constant volume : advection through the surface 205 DO jj = 2, jpjm1 206 DO ji = fs_2, fs_jpim1 207 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) 208 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) 209 END DO 210 END DO 211 ENDIF 212 DO jk = 2, jpkm1 ! interior fluxes 213 DO jj = 2, jpjm1 214 DO ji = fs_2, fs_jpim1 215 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 216 END DO 217 END DO 218 DO jj = 2, jpjm1 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 221 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 222 END DO 223 END DO 224 END DO 225 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence 226 DO jj = 2, jpjm1 227 DO ji = fs_2, fs_jpim1 ! vector opt. 228 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) 229 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) 230 END DO 231 END DO 232 END DO 233 ! 234 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 243 235 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 244 236 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 245 237 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 246 238 ENDIF 247 ! 239 ! ! Control print 248 240 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask, & 249 241 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 250 242 ! 251 CALL wrk_dealloc( jpi, jpj, jpk,zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )252 CALL wrk_dealloc( jpi, jpj, jpk, jpts,zlu_uu, zlv_vv, zlu_uv, zlv_vu )243 CALL wrk_dealloc( jpi,jpj,jpk, zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 244 CALL wrk_dealloc( jpi,jpj,jpk,jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu ) 253 245 ! 254 246 IF( nn_timing == 1 ) CALL timing_stop('dyn_adv_ubs') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r5120 r6140 29 29 30 30 !! * Substitutions 31 # include "domzgr_substitute.h90"32 # include "zdfddm_substitute.h90"33 31 # include "vectopt_loop_substitute.h90" 34 32 !!---------------------------------------------------------------------- … … 65 63 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 66 64 67 IF( l_trddyn ) THEN ! temporary save of ua and vatrends68 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )65 IF( l_trddyn ) THEN ! trends: store the input trends 66 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 69 67 ztrdu(:,:,:) = ua(:,:,:) 70 68 ztrdv(:,:,:) = va(:,:,:) … … 78 76 ! 79 77 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 80 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu)81 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv)78 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 79 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 82 80 END DO 83 81 END DO 84 85 IF ( ln_isfcav ) THEN82 ! 83 IF( ln_isfcav ) THEN ! ocean cavities 86 84 DO jj = 2, jpjm1 87 85 DO ji = 2, jpim1 … … 91 89 ! 92 90 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 93 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) &91 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) & 94 92 & * (1.-umask(ji,jj,1)) 95 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) &93 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) & 96 94 & * (1.-vmask(ji,jj,1)) 97 95 ! (ISF) … … 99 97 END DO 100 98 END IF 101 102 99 ! 103 IF( l_trddyn ) THEN ! save the vertical diffusive trendsfor further diagnostics100 IF( l_trddyn ) THEN ! trends: send trends to trddyn for further diagnostics 104 101 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 105 102 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 106 103 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 107 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )104 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 108 105 ENDIF 109 106 ! ! print mean trends (used for debugging) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5930 r6140 45 45 USE wrk_nemo ! Memory Allocation 46 46 USE timing ! Timing 47 USE iom 47 48 48 49 IMPLICIT NONE … … 52 53 PUBLIC dyn_hpg_init ! routine called by opa module 53 54 54 ! 55 LOGICAL , PUBLIC :: ln_hpg_zco 56 LOGICAL , PUBLIC :: ln_hpg_zps 57 LOGICAL , PUBLIC :: ln_hpg_sco 58 LOGICAL , PUBLIC :: ln_hpg_djc 59 LOGICAL , PUBLIC :: ln_hpg_prj 60 LOGICAL , PUBLIC :: ln_hpg_isf 55 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 56 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 57 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 58 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 59 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 60 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 61 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 61 62 62 63 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 63 64 64 65 !! * Substitutions 65 # include "domzgr_substitute.h90"66 66 # include "vectopt_loop_substitute.h90" 67 67 !!---------------------------------------------------------------------- … … 131 131 INTEGER :: ios ! Local integer output status for namelist read 132 132 !! 133 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 134 REAL(wp) :: znad 135 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd ! hypothesys on isf density 136 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 137 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 138 !! 133 139 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 134 140 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf … … 137 143 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient 138 144 READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) 139 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )140 145 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp ) 146 ! 141 147 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient 142 148 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 143 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )149 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 144 150 IF(lwm) WRITE ( numond, namdyn_hpg ) 145 151 ! … … 162 168 & either ln_hpg_sco or ln_hpg_prj instead') 163 169 ! 164 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )&165 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&166 & the standard jacobian formulation hpg_sco or&167 & the pressure jacobian formulation hpg_prj')170 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 171 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 172 & ' the standard jacobian formulation hpg_sco or ' , & 173 & ' the pressure jacobian formulation hpg_prj' ) 168 174 169 175 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & … … 190 196 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 191 197 ! 192 ! initialisation of ice load 193 riceload(:,:)=0.0 198 ! initialisation of ice shelf load 199 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 200 IF ( ln_isfcav ) THEN 201 CALL wrk_alloc( jpi,jpj, 2, ztstop) 202 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 203 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 204 ! 205 IF(lwp) WRITE(numout,*) 206 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 207 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 208 209 ! To use density and not density anomaly 210 znad=1._wp 211 212 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 213 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 214 215 ! compute density of the water displaced by the ice shelf 216 DO jk = 1, jpk 217 CALL eos(ztstop(:,:,:),gdept_n(:,:,jk),zrhd(:,:,jk)) 218 END DO 219 220 ! compute rhd at the ice/oce interface (ice shelf side) 221 CALL eos(ztstop,risfdep,zrhdtop_isf) 222 223 ! Surface value + ice shelf gradient 224 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 225 ! divided by 2 later 226 ziceload = 0._wp 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ikt=mikt(ji,jj) 230 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 231 DO jk=2,ikt-1 232 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 233 & * (1._wp - tmask(ji,jj,jk)) 234 END DO 235 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 236 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 237 END DO 238 END DO 239 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 240 241 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 242 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 243 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 244 END IF 194 245 ! 195 246 END SUBROUTINE dyn_hpg_init … … 213 264 !!---------------------------------------------------------------------- 214 265 INTEGER, INTENT(in) :: kt ! ocean time-step index 215 ! !266 ! 216 267 INTEGER :: ji, jj, jk ! dummy loop indices 217 268 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars … … 219 270 !!---------------------------------------------------------------------- 220 271 ! 221 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )272 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 222 273 ! 223 274 IF( kt == nit000 ) THEN … … 232 283 DO jj = 2, jpjm1 233 284 DO ji = fs_2, fs_jpim1 ! vector opt. 234 zcoef1 = zcoef0 * fse3w(ji,jj,1)285 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 235 286 ! hydrostatic pressure gradient 236 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) /e1u(ji,jj)237 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)287 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 288 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 238 289 ! add to the general momentum trend 239 290 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 247 298 DO jj = 2, jpjm1 248 299 DO ji = fs_2, fs_jpim1 ! vector opt. 249 zcoef1 = zcoef0 * fse3w(ji,jj,jk)300 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 250 301 ! hydrostatic pressure gradient 251 302 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 252 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &253 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)303 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 304 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 254 305 255 306 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 256 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &257 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)307 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 308 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 258 309 ! add to the general momentum trend 259 310 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 263 314 END DO 264 315 ! 265 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )316 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 266 317 ! 267 318 END SUBROUTINE hpg_zco … … 284 335 !!---------------------------------------------------------------------- 285 336 ! 286 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )337 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 287 338 ! 288 339 IF( kt == nit000 ) THEN … … 301 352 DO jj = 2, jpjm1 302 353 DO ji = fs_2, fs_jpim1 ! vector opt. 303 zcoef1 = zcoef0 * fse3w(ji,jj,1)354 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 304 355 ! hydrostatic pressure gradient 305 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) /e1u(ji,jj)306 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)356 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 357 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 307 358 ! add to the general momentum trend 308 359 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 310 361 END DO 311 362 END DO 312 313 363 314 364 ! interior value (2=<jk=<jpkm1) … … 316 366 DO jj = 2, jpjm1 317 367 DO ji = fs_2, fs_jpim1 ! vector opt. 318 zcoef1 = zcoef0 * fse3w(ji,jj,jk)368 zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 319 369 ! hydrostatic pressure gradient 320 370 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 321 371 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 322 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)372 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 323 373 324 374 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 325 375 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 326 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)376 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 327 377 ! add to the general momentum trend 328 378 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 331 381 END DO 332 382 END DO 333 334 383 335 384 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) … … 338 387 iku = mbku(ji,jj) 339 388 ikv = mbkv(ji,jj) 340 zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) )341 zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) )389 zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj ,iku) ) 390 zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) ) 342 391 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 343 392 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 344 393 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 345 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) /e1u(ji,jj)394 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 346 395 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 347 396 ENDIF … … 349 398 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 350 399 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 351 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) /e2v(ji,jj)400 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 352 401 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 353 402 ENDIF … … 355 404 END DO 356 405 ! 357 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )406 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 358 407 ! 359 408 END SUBROUTINE hpg_zps 409 360 410 361 411 SUBROUTINE hpg_sco( kt ) … … 391 441 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 392 442 ENDIF 393 394 ! Local constant initialization 443 ! 395 444 zcoef0 = - grav * 0.5_wp 396 ! To use density and not density anomaly 397 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 398 ELSE ; znad = 0._wp ! Fixed volume 445 IF ( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 446 ELSE ; znad = 1._wp ! Variable volume: density 399 447 ENDIF 400 448 ! 401 449 ! Surface value 402 450 DO jj = 2, jpjm1 403 451 DO ji = fs_2, fs_jpim1 ! vector opt. 404 452 ! hydrostatic pressure gradient along s-surfaces 405 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) )&406 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))407 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) )&408 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ))453 zhpi(ji,jj,1) = zcoef0 * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 454 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 455 zhpj(ji,jj,1) = zcoef0 * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 456 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 409 457 ! s-coordinate pressure gradient correction 410 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd(ji,jj,1) + 2._wp * znad ) &411 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) /e1u(ji,jj)412 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd(ji,jj,1) + 2._wp * znad ) &413 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) /e2v(ji,jj)458 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 459 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 460 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 461 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 414 462 ! add to the general momentum trend 415 463 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 423 471 DO ji = fs_2, fs_jpim1 ! vector opt. 424 472 ! hydrostatic pressure gradient along s-surfaces 425 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 /e1u(ji,jj) &426 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &427 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )428 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 /e2v(ji,jj) &429 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &430 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )473 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 474 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 475 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 476 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 477 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 478 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 431 479 ! s-coordinate pressure gradient correction 432 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd(ji,jj,jk) + 2._wp * znad ) &433 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) /e1u(ji,jj)434 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd(ji,jj,jk) + 2._wp * znad ) &435 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) /e2v(ji,jj)480 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 481 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 482 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 483 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 436 484 ! add to the general momentum trend 437 485 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 445 493 END SUBROUTINE hpg_sco 446 494 495 447 496 SUBROUTINE hpg_isf( kt ) 448 497 !!--------------------------------------------------------------------- 449 !! *** ROUTINE hpg_ sco***498 !! *** ROUTINE hpg_isf *** 450 499 !! 451 500 !! ** Method : s-coordinate case. Jacobian scheme. … … 466 515 INTEGER, INTENT(in) :: kt ! ocean time-step index 467 516 !! 468 INTEGER :: ji, jj, jk, ik u, ikv, ikt, iktp1i, iktp1j! dummy loop indices469 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars470 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd517 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 518 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 519 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 471 520 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 472 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 473 !!---------------------------------------------------------------------- 474 ! 475 CALL wrk_alloc( jpi,jpj, 2, ztstop) 476 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 477 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 478 ! 479 IF( kt == nit000 ) THEN 480 IF(lwp) WRITE(numout,*) 481 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 482 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 483 ENDIF 484 521 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 522 !!---------------------------------------------------------------------- 523 ! 524 CALL wrk_alloc( jpi,jpj, 2, ztstop) 525 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 526 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 527 ! 485 528 ! Local constant initialization 486 529 zcoef0 = - grav * 0.5_wp 530 487 531 ! To use density and not density anomaly 488 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume489 ! ELSE ; znad = 0._wp ! Fixed volume490 ! ENDIF491 532 znad=1._wp 533 492 534 ! iniitialised to 0. zhpi zhpi 493 535 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 536 495 ! Partial steps: top & bottom before horizontal gradient496 !jc CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi, &497 !jc & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &498 !jc & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi )499 500 !==================================================================================501 !=====Compute iceload and contribution of the half first wet layer =================502 !===================================================================================503 504 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)505 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp506 507 ! compute density of the water displaced by the ice shelf508 zrhd = rhd ! save rhd509 DO jk = 1, jpk510 zdept(:,:)=gdept_1d(jk)511 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))512 END DO513 WHERE ( tmask(:,:,:) == 1._wp)514 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd515 END WHERE516 517 ! compute rhd at the ice/oce interface (ice shelf side)518 CALL eos(ztstop,risfdep,zrhdtop_isf)519 520 537 ! compute rhd at the ice/oce interface (ocean side) 538 ! usefull to reduce residual current in the test case ISOMIP with no melting 521 539 DO ji=1,jpi 522 540 DO jj=1,jpj … … 526 544 END DO 527 545 END DO 528 CALL eos(ztstop,risfdep,zrhdtop_oce) 529 ! 530 ! Surface value + ice shelf gradient 531 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 532 ziceload = 0._wp 533 DO jj = 1, jpj 534 DO ji = 1, jpi ! vector opt. 535 ikt=mikt(ji,jj) 536 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 537 DO jk=2,ikt-1 538 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 539 & * (1._wp - tmask(ji,jj,jk)) 540 END DO 541 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 542 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 543 END DO 544 END DO 545 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 546 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 546 CALL eos( ztstop, risfdep, zrhdtop_oce ) 547 548 !================================================================================== 549 !===== Compute surface value ===================================================== 550 !================================================================================== 547 551 DO jj = 2, jpjm1 548 552 DO ji = fs_2, fs_jpim1 ! vector opt. 549 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 553 ikt = mikt(ji,jj) 554 iktp1i = mikt(ji+1,jj) 555 iktp1j = mikt(ji,jj+1) 550 556 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 551 557 ! we assume ISF is in isostatic equilibrium 552 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj,iktp1i) &553 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj) ) &554 & - 0.5_wp * fse3w(ji ,jj ,ikt )&555 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&556 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)))557 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j) &558 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) &559 & - 0.5_wp * fse3w(ji ,jj ,ikt )&560 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&561 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ))558 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i) & 559 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 560 & - 0.5_wp * e3w_n(ji,jj,ikt) & 561 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 562 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 563 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j) & 564 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 565 & - 0.5_wp * e3w_n(ji,jj,ikt) & 566 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 567 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 562 568 ! s-coordinate pressure gradient correction (=0 if z coordinate) 563 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd(ji,jj,1) + 2._wp * znad ) &564 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) /e1u(ji,jj)565 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd(ji,jj,1) + 2._wp * znad ) &566 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) /e2v(ji,jj)569 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 570 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 571 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 572 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 567 573 ! add to the general momentum trend 568 574 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) … … 571 577 END DO 572 578 !================================================================================== 573 !===== Compute partial cell contribution for the top cell =========================574 !==================================================================================575 DO jj = 2, jpjm1576 DO ji = fs_2, fs_jpim1 ! vector opt.577 iku = miku(ji,jj) ;578 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp579 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))580 ! u direction581 IF ( iku .GT. 1 ) THEN582 ! case iku583 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &584 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &585 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )586 ! corrective term ( = 0 if z coordinate )587 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)588 ! zhpi will be added in interior loop589 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap590 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure591 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)592 593 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)594 zhpiint = zcoef0 / e1u(ji,jj) &595 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &596 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &597 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &598 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )599 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint600 END IF601 602 ! v direction603 ikv = mikv(ji,jj)604 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))605 IF ( ikv .GT. 1 ) THEN606 ! case ikv607 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &608 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &609 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )610 ! corrective term ( = 0 if z coordinate )611 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)612 ! zhpi will be added in interior loop613 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap614 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure615 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)616 617 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)618 zhpjint = zcoef0 / e2v(ji,jj) &619 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &620 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &621 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &622 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )623 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint624 END IF625 END DO626 END DO627 628 !==================================================================================629 579 !===== Compute interior value ===================================================== 630 580 !================================================================================== 631 632 DO jj = 2, jpjm1 633 DO ji = fs_2, fs_jpim1 ! vector opt. 634 iku=miku(ji,jj); ikv=mikv(ji,jj) 635 DO jk = 2, jpkm1 581 ! interior value (2=<jk=<jpkm1) 582 DO jk = 2, jpkm1 583 DO jj = 2, jpjm1 584 DO ji = fs_2, fs_jpim1 ! vector opt. 636 585 ! hydrostatic pressure gradient along s-surfaces 637 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 638 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 639 & + zcoef0 / e1u(ji,jj) & 640 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 641 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 642 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 643 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 586 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 587 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 588 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 589 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 590 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 591 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 644 592 ! s-coordinate pressure gradient correction 645 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 646 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 647 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 648 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 649 650 ! hydrostatic pressure gradient along s-surfaces 651 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 652 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 653 & + zcoef0 / e2v(ji,jj) & 654 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 655 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 656 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 657 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 658 ! s-coordinate pressure gradient correction 659 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 660 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 661 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 593 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 594 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 595 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 596 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 662 597 ! add to the general momentum trend 663 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 598 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 599 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 664 600 END DO 665 601 END DO 666 602 END DO 667 668 !================================================================================== 669 !===== Compute bottom cell contribution (partial cell) ============================ 670 !================================================================================== 671 672 DO jj = 2, jpjm1 673 DO ji = 2, jpim1 674 iku = mbku(ji,jj) 675 ikv = mbkv(ji,jj) 676 677 IF (iku .GT. 1) THEN 678 ! remove old value (interior case) 679 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 680 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 681 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 682 ! put new value 683 ! -zpshpi to avoid double contribution of the partial step in the top layer 684 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 685 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 686 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 687 END IF 688 ! v direction 689 IF (ikv .GT. 1) THEN 690 ! remove old value (interior case) 691 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 692 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 693 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 694 ! put new value 695 ! -zpshpj to avoid double contribution of the partial step in the top layer 696 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 697 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 698 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 699 END IF 700 END DO 701 END DO 702 703 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 704 rhd = zrhd 705 ! 706 CALL wrk_dealloc( jpi,jpj,2, ztstop) 707 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 708 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 603 ! 604 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 605 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 606 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 709 607 ! 710 608 END SUBROUTINE hpg_isf … … 756 654 DO jj = 2, jpjm1 757 655 DO ji = fs_2, fs_jpim1 ! vector opt. 758 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd(ji,jj,jk-1)759 dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1)760 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd(ji,jj,jk )761 dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk )762 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd(ji,jj,jk )763 dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk )656 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 657 dzz (ji,jj,jk) = gde3w_n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1) 658 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 659 dzx (ji,jj,jk) = gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk ) 660 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 661 dzy (ji,jj,jk) = gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk ) 764 662 END DO 765 663 END DO … … 843 741 !------------------------------------------------------------- 844 742 845 !!bug gm : e3w- de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified846 ! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be743 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified 744 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 847 745 848 746 DO jj = 2, jpjm1 849 747 DO ji = fs_2, fs_jpim1 ! vector opt. 850 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &851 & * ( rhd(ji,jj,1) &852 & + 0.5_wp * ( rhd (ji,jj,2) - rhd(ji,jj,1) )&853 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )&854 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )748 rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) ) & 749 & * ( rhd(ji,jj,1) & 750 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 751 & * ( e3w_n (ji,jj,1) - gde3w_n(ji,jj,1) ) & 752 & / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) ) ) 855 753 END DO 856 754 END DO … … 863 761 DO ji = fs_2, fs_jpim1 ! vector opt. 864 762 865 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd(ji,jj,jk-1) ) &866 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) &867 & - grav * z1_10 * ( 868 & ( drhow (ji,jj,jk) - drhow(ji,jj,jk-1) ) &869 & * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &870 & - ( dzw (ji,jj,jk) - dzw(ji,jj,jk-1) ) &871 & * ( rhd (ji,jj,jk) - rhd(ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &763 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 764 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) & 765 & - grav * z1_10 * ( & 766 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 767 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 768 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 769 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 872 770 & ) 873 771 874 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd(ji,jj,jk) ) &875 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) &876 & - grav* z1_10 * ( 877 & ( drhou (ji+1,jj,jk) - drhou(ji,jj,jk) ) &878 & * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &879 & - ( dzu (ji+1,jj,jk) - dzu(ji,jj,jk) ) &880 & * ( rhd (ji+1,jj,jk) - rhd(ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &772 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 773 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) & 774 & - grav* z1_10 * ( & 775 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 776 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 777 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 778 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 881 779 & ) 882 780 883 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) )&884 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) &885 & - grav* z1_10 * ( 886 & ( drhov (ji,jj+1,jk) - drhov(ji,jj,jk) ) &887 & * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &888 & - ( dzv (ji,jj+1,jk) - dzv(ji,jj,jk) ) &889 & * ( rhd (ji,jj+1,jk) - rhd(ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &781 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 782 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) & 783 & - grav* z1_10 * ( & 784 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 785 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 786 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 787 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 890 788 & ) 891 789 … … 903 801 DO jj = 2, jpjm1 904 802 DO ji = fs_2, fs_jpim1 ! vector opt. 905 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) /e1u(ji,jj)906 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) /e2v(ji,jj)803 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 804 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 907 805 ! add to the general momentum trend 908 806 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 920 818 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 921 819 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 922 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) /e1u(ji,jj)820 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 923 821 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 924 822 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 925 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) /e2v(ji,jj)823 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 926 824 ! add to the general momentum trend 927 825 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 953 851 !! 954 852 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 955 REAL(wp) :: zcoef0, znad ! temporaryscalars956 ! !853 REAL(wp) :: zcoef0, znad ! local scalars 854 ! 957 855 !! The local variables for the correction term 958 856 INTEGER :: jk1, jis, jid, jjs, jjd … … 965 863 !!---------------------------------------------------------------------- 966 864 ! 967 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )968 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )969 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )865 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 866 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 867 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 970 868 ! 971 869 IF( kt == nit000 ) THEN … … 975 873 ENDIF 976 874 977 !!----------------------------------------------------------------------978 875 ! Local constant initialization 979 876 zcoef0 = - grav 980 znad = 0.0_wp981 IF( l k_vvl ) znad = 1._wp877 znad = 1._wp 878 IF( ln_linssh ) znad = 0._wp 982 879 983 880 ! Clean 3-D work arrays … … 989 886 DO ji = 1, jpi 990 887 jk = mbathy(ji,jj) 991 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp992 ELSE IF(jk == 1) THEN; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)993 ELSE IF(jk < jpkm1) THEN888 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 889 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 890 ELSEIF( jk < jpkm1 ) THEN 994 891 DO jkk = jk+1, jpk 995 zrhh(ji,jj,jkk) = interp1( fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),&996 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))892 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), & 893 & gde3w_n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 997 894 END DO 998 895 ENDIF … … 1003 900 DO jj = 1, jpj 1004 901 DO ji = 1, jpi 1005 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad902 zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 1006 903 END DO 1007 904 END DO … … 1010 907 DO jj = 1, jpj 1011 908 DO ji = 1, jpi 1012 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)909 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 1013 910 END DO 1014 911 END DO … … 1021 918 ! constrained cubic spline interpolation 1022 919 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 1023 CALL cspline( fsp,xsp,asp,bsp,csp,dsp,polynomial_type)920 CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1024 921 1025 922 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1026 923 DO jj = 2, jpj 1027 924 DO ji = 2, jpi 1028 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 1029 bsp(ji,jj,1), csp(ji,jj,1), & 1030 dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1) 925 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 926 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1031 927 1032 928 ! assuming linear profile across the top half surface layer 1033 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1929 zhpi(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) * zrhdt1 1034 930 END DO 1035 931 END DO … … 1039 935 DO jj = 2, jpj 1040 936 DO ji = 2, jpi 1041 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &1042 integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&1043 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), &1044 csp(ji,jj,jk-1), dsp(ji,jj,jk-1))937 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 938 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 939 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 940 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1045 941 END DO 1046 942 END DO … … 1052 948 DO jj = 2, jpjm1 1053 949 DO ji = 2, jpim1 950 !!gm BUG ? if it is ssh at u- & v-point then it should be: 951 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 952 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 953 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 954 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 955 !!gm not this: 1054 956 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1055 957 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp … … 1061 963 DO jj = 2, jpjm1 1062 964 DO ji = 2, jpim1 1063 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad)1064 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad)965 zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 966 zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 1065 967 END DO 1066 968 END DO … … 1069 971 DO jj = 2, jpjm1 1070 972 DO ji = 2, jpim1 1071 zu(ji,jj,jk) = zu(ji,jj,jk-1) - fse3u(ji,jj,jk)1072 zv(ji,jj,jk) = zv(ji,jj,jk-1) - fse3v(ji,jj,jk)973 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 974 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 1073 975 END DO 1074 976 END DO … … 1078 980 DO jj = 2, jpjm1 1079 981 DO ji = 2, jpim1 1080 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)1081 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)982 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 983 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 1082 984 END DO 1083 985 END DO … … 1087 989 DO jj = 2, jpjm1 1088 990 DO ji = 2, jpim1 1089 zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1090 zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1091 zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))1092 zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))991 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 992 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 993 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 994 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1093 995 END DO 1094 996 END DO … … 1148 1050 ! update the momentum trends in u direction 1149 1051 1150 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))1151 IF( lk_vvl) THEN1152 zdpdx2 = zcoef0 /e1u(ji,jj) * &1153 1052 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1053 IF( .NOT.ln_linssh ) THEN 1054 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1055 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1154 1056 ELSE 1155 zdpdx2 = zcoef0 /e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)1057 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1156 1058 ENDIF 1157 1158 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1159 &umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)1059 !!gm Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:) by definition 1060 !!gm in the line below only * umask(ji,jj,jk) is needed !! 1061 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1160 1062 ENDIF 1161 1063 … … 1205 1107 ! update the momentum trends in v direction 1206 1108 1207 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))1208 IF( lk_vvl) THEN1209 zdpdy2 = zcoef0 /e2v(ji,jj) * &1210 1109 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1110 IF( .NOT.ln_linssh ) THEN 1111 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1112 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1211 1113 ELSE 1212 zdpdy2 = zcoef0 /e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )1114 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1213 1115 ENDIF 1214 1215 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1216 &vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)1116 !!gm Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:) by definition 1117 !!gm in the line below only * vmask(ji,jj,jk) is needed !! 1118 va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1217 1119 ENDIF 1218 1219 1220 END DO 1221 END DO 1222 END DO 1223 ! 1224 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1225 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1226 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1120 ! 1121 END DO 1122 END DO 1123 END DO 1124 ! 1125 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1126 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1127 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1227 1128 ! 1228 1129 END SUBROUTINE hpg_prj 1229 1130 1230 1131 1231 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type)1132 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1232 1133 !!---------------------------------------------------------------------- 1233 1134 !! *** ROUTINE cspline *** … … 1239 1140 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1240 1141 !!---------------------------------------------------------------------- 1241 IMPLICIT NONE 1242 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1243 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1244 ! the interpoated function 1245 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1246 ! 2: Linear 1142 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate 1143 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1144 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1247 1145 ! 1248 1146 INTEGER :: ji, jj, jk ! dummy loop indices … … 1252 1150 REAL(wp) :: zdf(size(fsp,3)) 1253 1151 !!---------------------------------------------------------------------- 1254 1152 ! 1153 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1255 1154 jpi = size(fsp,1) 1256 1155 jpj = size(fsp,2) 1257 1156 jpkm1 = size(fsp,3) - 1 1258 1259 1157 ! 1260 1158 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1261 1159 DO ji = 1, jpi … … 1290 1188 1291 1189 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1292 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2)1190 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1293 1191 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1294 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1295 & 0.5_wp * zdf(jpkm1 - 1) 1192 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1296 1193 1297 1194 DO jk = 1, jpkm1 - 1 … … 1316 1213 END DO 1317 1214 1318 ELSE IF (polynomial_type == 2) THEN ! Linear1215 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1319 1216 DO ji = 1, jpi 1320 1217 DO jj = 1, jpj … … 1347 1244 !! extrapolation is also permitted (no value limit) 1348 1245 !!---------------------------------------------------------------------- 1349 IMPLICIT NONE1350 1246 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1351 1247 REAL(wp) :: f ! result of the interpolation (extrapolation) 1352 1248 REAL(wp) :: zdeltx 1353 1249 !!---------------------------------------------------------------------- 1354 1250 ! 1355 1251 zdeltx = xr - xl 1356 IF( abs(zdeltx) <= 10._wp * EPSILON(x)) THEN1357 f = 0.5_wp * (fl + fr)1252 IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN 1253 f = 0.5_wp * (fl + fr) 1358 1254 ELSE 1359 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx1255 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1360 1256 ENDIF 1361 1257 ! 1362 1258 END FUNCTION interp1 1363 1259 1364 1260 1365 FUNCTION interp2( x, a, b, c, d) RESULT(f)1261 FUNCTION interp2( x, a, b, c, d ) RESULT(f) 1366 1262 !!---------------------------------------------------------------------- 1367 1263 !! *** ROUTINE interp1 *** … … 1372 1268 !! 1373 1269 !!---------------------------------------------------------------------- 1374 IMPLICIT NONE1375 1270 REAL(wp), INTENT(in) :: x, a, b, c, d 1376 1271 REAL(wp) :: f ! value from the interpolation 1377 1272 !!---------------------------------------------------------------------- 1378 1273 ! 1379 1274 f = a + x* ( b + x * ( c + d * x ) ) 1380 1275 ! 1381 1276 END FUNCTION interp2 1382 1277 1383 1278 1384 FUNCTION interp3( x, a, b, c, d) RESULT(f)1279 FUNCTION interp3( x, a, b, c, d ) RESULT(f) 1385 1280 !!---------------------------------------------------------------------- 1386 1281 !! *** ROUTINE interp1 *** … … 1392 1287 !! 1393 1288 !!---------------------------------------------------------------------- 1394 IMPLICIT NONE1395 1289 REAL(wp), INTENT(in) :: x, a, b, c, d 1396 1290 REAL(wp) :: f ! value from the interpolation 1397 1291 !!---------------------------------------------------------------------- 1398 1292 ! 1399 1293 f = b + x * ( 2._wp * c + 3._wp * d * x) 1400 1294 ! 1401 1295 END FUNCTION interp3 1402 1296 1403 1297 1404 FUNCTION integ_spline( xl, xr, a, b, c, d) RESULT(f)1298 FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f) 1405 1299 !!---------------------------------------------------------------------- 1406 1300 !! *** ROUTINE interp1 *** … … 1411 1305 !! 1412 1306 !!---------------------------------------------------------------------- 1413 IMPLICIT NONE1414 1307 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1415 1308 REAL(wp) :: za1, za2, za3 1416 1309 REAL(wp) :: f ! integration result 1417 1310 !!---------------------------------------------------------------------- 1418 1311 ! 1419 1312 za1 = 0.5_wp * b 1420 1313 za2 = c / 3.0_wp 1421 1314 za3 = 0.25_wp * d 1422 1315 ! 1423 1316 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1424 1317 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1425 1318 ! 1426 1319 END FUNCTION integ_spline 1427 1320 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r5836 r6140 36 36 PUBLIC dyn_ldf_init ! called by opa module 37 37 38 ! ! Flagto control the type of lateral viscous operator38 ! ! Parameter to control the type of lateral viscous operator 39 39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in setting the operator 40 40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral viscous trend) … … 46 46 47 47 !! * Substitutions 48 # include "domzgr_substitute.h90"49 48 # include "vectopt_loop_substitute.h90" 50 49 !!---------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r5836 r6140 41 41 42 42 !! * Substitutions 43 # include "domzgr_substitute.h90"44 43 # include "vectopt_loop_substitute.h90" 45 44 !!---------------------------------------------------------------------- … … 135 134 DO jj = 2, jpjm1 136 135 DO ji = 2, jpim1 137 uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)138 vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)139 wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5140 wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5136 uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 137 vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 138 wslpi(ji,jj,jk) = - ( gdepw_b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 139 wslpj(ji,jj,jk) = - ( gdepw_b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 141 140 END DO 142 141 END DO … … 183 182 DO jj = 2, jpjm1 184 183 DO ji = fs_2, jpi ! vector opt. 185 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj)184 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u_n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) 186 185 187 186 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & … … 198 197 DO jj = 2, jpjm1 199 198 DO ji = fs_2, jpi ! vector opt. 200 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj)199 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj) 201 200 202 201 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & … … 215 214 DO jj = 1, jpjm1 216 215 DO ji = 1, fs_jpim1 ! vector opt. 217 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj)216 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj) 218 217 219 218 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & … … 236 235 DO jj = 2, jpjm1 237 236 DO ji = 1, fs_jpim1 ! vector opt. 238 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj)237 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj) 239 238 240 239 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & … … 253 252 DO jj = 2, jpj 254 253 DO ji = 1, fs_jpim1 ! vector opt. 255 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj)254 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v_n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) 256 255 257 256 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 268 267 DO jj = 2, jpj 269 268 DO ji = 1, fs_jpim1 ! vector opt. 270 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj)269 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj) 271 270 272 271 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 287 286 DO jj = 2, jpjm1 288 287 DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug 289 ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) &290 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk))291 va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) &292 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))288 ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 289 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 290 va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 291 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 293 292 END DO 294 293 END DO … … 403 402 DO jk = 1, jpkm1 404 403 DO ji = 2, jpim1 405 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk))406 va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))404 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 405 va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 407 406 END DO 408 407 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90
r5861 r6140 4 4 !! Ocean dynamics: lateral viscosity trend (laplacian and bilaplacian) 5 5 !!====================================================================== 6 !! History : OPA ! 1990-09 (G. Madec) Original code 7 !! 4.0 ! 1991-11 (G. Madec) 8 !! 6.0 ! 1996-01 (G. Madec) statement function for e3 and ahm 9 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 10 !! - ! 2004-08 (C. Talandier) New trends organization 11 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 12 !! ! add velocity dependent coefficient and optional read in file 6 !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian 13 7 !!---------------------------------------------------------------------- 14 8 … … 35 29 36 30 !! * Substitutions 37 # include "domzgr_substitute.h90"38 31 # include "vectopt_loop_substitute.h90" 39 32 !!---------------------------------------------------------------------- … … 87 80 DO ji = fs_2, jpi ! vector opt. 88 81 ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) 89 !!gm open question here : fse3f at before or now ? probably now...82 !!gm open question here : e3f at before or now ? probably now... 90 83 !!gm note that ahmf has already been multiplied by fmask 91 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * fse3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) &84 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 92 85 & * ( e2v(ji ,jj-1) * pvb(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) & 93 & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) 86 & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) 94 87 ! ! ahm * div (computed from 2 to jpi/jpj) 95 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / fse3t_b(ji,jj,jk) * tmask(ji,jj,jk) & 96 & * ( e2u(ji,jj)*fse3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*fse3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & 97 & + e1v(ji,jj)*fse3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*fse3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) 88 !!gm note that ahmt has already been multiplied by tmask 89 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) & 90 & * ( e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & 91 & + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) 98 92 END DO 99 93 END DO … … 101 95 DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) 102 96 DO ji = fs_2, fs_jpim1 ! vector opt. 103 pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( 104 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk)) &97 pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( & 98 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) & 105 99 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 106 100 ! 107 pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( 108 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk)) &101 pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( & 102 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) & 109 103 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 110 104 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5930 r6140 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 !! 3. 7! 2014-04 (G. Madec) add the diagnostic of the time filter trends20 !! 3.6 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 22 22 !!------------------------------------------------------------------------- 23 23 24 24 !!------------------------------------------------------------------------- 25 !! dyn_nxt : obtain the next (after) horizontal velocity25 !! dyn_nxt : obtain the next (after) horizontal velocity 26 26 !!------------------------------------------------------------------------- 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 USE phycst ! physical constants 31 USE dynadv ! dynamics: vector invariant versus flux form 32 USE domvvl ! variable volume 33 USE bdy_oce ! ocean open boundary conditions 34 USE bdydta ! ocean open boundary conditions 35 USE bdydyn ! ocean open boundary conditions 36 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 37 USE trd_oce ! trends: ocean variables 38 USE trddyn ! trend manager: dynamics 39 USE trdken ! trend manager: kinetic energy 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 USE phycst ! physical constants 31 USE dynadv ! dynamics: vector invariant versus flux form 32 USE dynspg_ts ! surface pressure gradient: split-explicit scheme 33 USE domvvl ! variable volume 34 USE bdy_oce ! ocean open boundary conditions 35 USE bdydta ! ocean open boundary conditions 36 USE bdydyn ! ocean open boundary conditions 37 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 38 USE trd_oce ! trends: ocean variables 39 USE trddyn ! trend manager: dynamics 40 USE trdken ! trend manager: kinetic energy 40 41 ! 41 USE in_out_manager 42 USE iom 43 USE lbclnk 44 USE lib_mpp 45 USE wrk_nemo 46 USE prtctl 47 USE timing 42 USE in_out_manager ! I/O manager 43 USE iom ! I/O manager library 44 USE lbclnk ! lateral boundary condition (or mpp link) 45 USE lib_mpp ! MPP library 46 USE wrk_nemo ! Memory Allocation 47 USE prtctl ! Print control 48 USE timing ! Timing 48 49 #if defined key_agrif 49 50 USE agrif_opa_interp … … 55 56 PUBLIC dyn_nxt ! routine called by step.F90 56 57 57 !! * Substitutions58 # include "domzgr_substitute.h90"59 58 !!---------------------------------------------------------------------- 60 59 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 85 84 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 86 85 !! (un,vn) = (ua,va). 87 !! Note that with flux form advection and variable volume layer88 !! (lk_vvl=T), the time filter is applied on thickness weighted89 !! velocity.86 !! Note that with flux form advection and non linear free surface, 87 !! the time filter is applied on thickness weighted velocity. 88 !! As a result, dyn_nxt MUST be called after tra_nxt. 90 89 !! 91 90 !! ** Action : ub,vb filtered before horizontal velocity of next time-step … … 95 94 ! 96 95 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: ik u, ikv! local integers98 REAL(wp) :: zue3a, zue3n, zue3b, zuf 96 INTEGER :: ikt ! local integers 97 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 99 98 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 100 99 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 104 103 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 105 104 ! 106 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 107 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 108 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva ) 105 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve) 106 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva) 109 107 ! 110 108 IF( kt == nit000 ) THEN … … 117 115 ! Ensure below that barotropic velocities match time splitting estimate 118 116 ! Compute actual transport and replace it with ts estimate at "after" time step 119 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)120 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)117 zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 118 zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 121 119 DO jk = 2, jpkm1 122 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)123 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)120 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 121 zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 124 122 END DO 125 123 DO jk = 1, jpkm1 126 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)127 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)124 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 125 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 128 126 END DO 129 130 IF 127 ! 128 IF( .NOT.ln_bt_fw ) THEN 131 129 ! Remove advective velocity from "now velocities" 132 130 ! prior to asselin filtering … … 151 149 # if defined key_bdy 152 150 ! !* BDY open boundaries 153 IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )154 IF( lk_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )151 IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 152 IF( lk_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 155 153 156 154 !!$ Do we need a call to bdy_vol here?? … … 184 182 vn(:,:,jk) = va(:,:,jk) 185 183 END DO 186 IF (lk_vvl) THEN184 IF(.NOT.ln_linssh ) THEN 187 185 DO jk = 1, jpkm1 188 fse3t_b(:,:,jk) = fse3t_n(:,:,jk)189 fse3u_b(:,:,jk) = fse3u_n(:,:,jk)190 fse3v_b(:,:,jk) = fse3v_n(:,:,jk)191 END DO186 e3t_b(:,:,jk) = e3t_n(:,:,jk) 187 e3u_b(:,:,jk) = e3u_n(:,:,jk) 188 e3v_b(:,:,jk) = e3v_n(:,:,jk) 189 END DO 192 190 ENDIF 193 191 ELSE !* Leap-Frog : Asselin filter and swap 194 192 ! ! =============! 195 IF( .NOT. lk_vvl ) THEN! Fixed volume !193 IF( ln_linssh ) THEN ! Fixed volume ! 196 194 ! ! =============! 197 195 DO jk = 1, jpkm1 … … 214 212 ! (used as a now filtered scale factor until the swap) 215 213 ! ---------------------------------------------------- 216 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 217 ! No asselin filtering on thicknesses if forward time splitting 218 fse3t_b(:,:,:) = fse3t_n(:,:,:) 214 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN ! No asselin filtering on thicknesses if forward time splitting 215 e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1) 219 216 ELSE 220 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 217 DO jk = 1, jpkm1 218 e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 219 END DO 221 220 ! Add volume filter correction: compatibility with tracer advection scheme 222 221 ! => time filter + conservation correction (only at the first level) 223 IF ( nn_isf == 0) THEN ! if no ice shelf melting 224 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 225 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 222 zcoef = atfp * rdt * r1_rau0 223 IF ( .NOT. ln_isf ) THEN ! if no ice shelf melting 224 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 225 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 226 226 ELSE ! if ice shelf melting 227 DO jj = 1,jpj 228 DO ji = 1,jpi 229 jk = mikt(ji,jj) 230 fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0 & 231 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 232 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 233 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk) 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ikt = mikt(ji,jj) 230 e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * ( emp_b (ji,jj) - emp (ji,jj) & 231 & - rnf_b (ji,jj) + rnf (ji,jj) & 232 & + fwfisf_b(ji,jj) - fwfisf(ji,jj) ) * tmask(ji,jj,ikt) 234 233 END DO 235 234 END DO … … 237 236 ENDIF 238 237 ! 239 IF( ln_dynadv_vec ) THEN 240 ! Before scale factor at (u/v)-points 241 ! ----------------------------------- 242 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 243 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 244 ! Leap-Frog - Asselin filter and swap: applied on velocity 245 ! ----------------------------------- 238 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 239 ! Before filtered scale factor at (u/v)-points 240 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 241 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 246 242 DO jk = 1, jpkm1 247 243 DO jj = 1, jpj … … 258 254 END DO 259 255 ! 260 ELSE 261 ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 262 !------------------------------------------------ 263 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 264 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 265 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 266 ! ----------------------------------- =========================== 256 ELSE ! Asselin filter applied on thickness weighted velocity 257 ! 258 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 259 ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 260 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 261 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 267 262 DO jk = 1, jpkm1 268 263 DO jj = 1, jpj 269 264 DO ji = 1, jpi 270 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)271 zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)272 zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)273 zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)274 zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)275 zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)265 zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) 266 zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) 267 zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) 268 zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) 269 zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) 270 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 276 271 ! 277 272 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) … … 285 280 END DO 286 281 END DO 287 fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 288 fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 282 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 283 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 284 ! 285 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 289 286 ENDIF 290 287 ! 291 288 ENDIF 292 289 ! 293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN290 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 294 291 ! Revert "before" velocities to time split estimate 295 292 ! Doing it here also means that asselin filter contribution is removed 296 zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)297 zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)293 zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 294 zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 298 295 DO jk = 2, jpkm1 299 zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)300 zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)296 zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 297 zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 301 298 END DO 302 299 DO jk = 1, jpkm1 303 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)304 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)300 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 301 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 305 302 END DO 306 303 ENDIF … … 313 310 ! 314 311 ! 315 IF (lk_vvl) THEN316 hu_b(:,:) = 0.317 hv_b(:,:) = 0.318 DO jk = 1, jpkm1319 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)320 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)312 IF(.NOT.ln_linssh ) THEN 313 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 314 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 315 DO jk = 2, jpkm1 316 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 317 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 321 318 END DO 322 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 323 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 324 ENDIF 325 ! 326 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 327 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 328 ! 329 DO jk = 1, jpkm1 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 333 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 334 ! 335 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 336 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 337 END DO 338 END DO 319 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 320 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 321 ENDIF 322 ! 323 un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 324 ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 325 vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) 326 vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 327 DO jk = 2, jpkm1 328 un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 329 ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 330 vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 331 vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 339 332 END DO 340 !341 !342 u n_b(:,:) = un_b(:,:) * hur_a(:,:)343 v n_b(:,:) = vn_b(:,:) * hvr_a(:,:)344 ub_b(:,:) = ub_b(:,:) * hur_b(:,:)345 vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)346 !347 !348 333 un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 334 vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) 335 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 336 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 337 ! 338 IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents 339 CALL iom_put( "ubar", un_b(:,:) ) 340 CALL iom_put( "vbar", vn_b(:,:) ) 341 ENDIF 349 342 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 350 343 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt … … 355 348 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 356 349 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 357 ! 358 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 359 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 360 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 350 ! 351 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 352 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 361 353 ! 362 354 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5930 r6140 6 6 !! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code 7 7 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 !! 3.7 ! 2015-11 (J. Chanut) Suppression of filtered free surface9 8 !!---------------------------------------------------------------------- 10 9 … … 21 20 USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine) 22 21 USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine) 23 USE sbctide 24 USE updtide 22 USE sbctide ! 23 USE updtide ! 25 24 USE trd_oce ! trends: ocean variables 26 25 USE trddyn ! trend manager: dynamics … … 32 31 USE timing ! Timing 33 32 34 35 33 IMPLICIT NONE 36 34 PRIVATE … … 39 37 PUBLIC dyn_spg_init ! routine called by opa module 40 38 41 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from ln_dynspg_... 39 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from lk_dynspg_... 40 41 ! ! Parameter to control the surface pressure gradient scheme 42 INTEGER, PARAMETER :: np_TS = 1 ! split-explicit time stepping (Time-Splitting) 43 INTEGER, PARAMETER :: np_EXP = 0 ! explicit time stepping 44 INTEGER, PARAMETER :: np_NO =-1 ! no surface pressure gradient, no scheme 42 45 43 46 !! * Substitutions 44 # include "domzgr_substitute.h90"45 47 # include "vectopt_loop_substitute.h90" 46 48 !!---------------------------------------------------------------------- … … 55 57 !! *** ROUTINE dyn_spg *** 56 58 !! 57 !! ** Purpose : achieve the momentum time stepping by computing the 58 !! last trend, the surface pressure gradient including the 59 !! atmospheric pressure forcing (ln_apr_dyn=T), and performing 60 !! the Leap-Frog integration. 61 !!gm In the current version only the filtered solution provide 62 !!gm the after velocity, in the 2 other (ua,va) are still the trends 59 !! ** Purpose : compute surface pressure gradient including the 60 !! atmospheric pressure forcing (ln_apr_dyn=T). 63 61 !! 64 62 !! ** Method : Two schemes: 65 !! - explicit computation: the spg is evaluated at now66 !! - split-explicit computation: a time splitting technique is used63 !! - explicit : the spg is evaluated at now 64 !! - split-explicit : a time splitting technique is used 67 65 !! 68 66 !! ln_apr_dyn=T : the atmospheric pressure forcing is applied … … 73 71 !! period is used to prevent the divergence of odd and even time step. 74 72 !!---------------------------------------------------------------------- 75 !76 73 INTEGER, INTENT(in ) :: kt ! ocean time-step index 77 74 ! … … 84 81 IF( nn_timing == 1 ) CALL timing_start('dyn_spg') 85 82 ! 86 87 !!gm NOTA BENE : the dynspg_exp and dynspg_ts should be modified so that88 !!gm they return the after velocity, not the trends (as in trazdf_imp...)89 !!gm In this case, change/simplify dynnxt90 91 92 83 IF( l_trddyn ) THEN ! temporary save of ta and sa trends 93 CALL wrk_alloc( jpi, jpj, jpk,ztrdu, ztrdv )84 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 94 85 ztrdu(:,:,:) = ua(:,:,:) 95 86 ztrdv(:,:,:) = va(:,:,:) 96 87 ENDIF 97 88 ! 98 89 IF( ln_apr_dyn & ! atmos. pressure 99 90 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting) … … 107 98 END DO 108 99 ! 109 IF( ln_apr_dyn .AND. (.NOT. ln_dynspg_ts) ) THEN!== Atmospheric pressure gradient (added later in time-split case) ==!100 IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==! 110 101 zg_2 = grav * 0.5 111 102 DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh … … 133 124 ! 134 125 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: Pressure gradient due to snow-ice mass ==! 135 CALL wrk_alloc( jpi, jpj,zpice )126 CALL wrk_alloc( jpi,jpj, zpice ) 136 127 ! 137 128 zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) … … 145 136 END DO 146 137 ! 147 CALL wrk_dealloc( jpi, jpj,zpice )138 CALL wrk_dealloc( jpi,jpj, zpice ) 148 139 ENDIF 149 140 ! 150 DO jk = 1, jpkm1 141 DO jk = 1, jpkm1 !== Add all terms to the general trend 151 142 DO jj = 2, jpjm1 152 143 DO ji = fs_2, fs_jpim1 ! vector opt. … … 156 147 END DO 157 148 END DO 158 149 ! 159 150 !!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ???? 160 161 ENDIF 162 163 SELECT CASE ( nspg ) ! compute surf. pressure gradient trend and add it to the general trend 164 ! 165 CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit 166 CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 167 ! 151 ! 152 ENDIF 153 ! 154 SELECT CASE ( nspg ) !== surface pressure gradient computed and add to the general trend ==! 155 CASE ( np_EXP ) ; CALL dyn_spg_exp( kt ) ! explicit 156 CASE ( np_TS ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 168 157 END SELECT 169 158 ! 170 IF( l_trddyn ) THEN 159 IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics 171 160 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 161 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 162 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 174 ! 175 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 176 ENDIF 177 ! ! print mean trends (used for debugging) 163 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 164 ENDIF 165 ! ! print mean trends (used for debugging) 178 166 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' spg - Ua: ', mask1=umask, & 179 167 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) … … 191 179 !! surface pressure gradient schemes 192 180 !!---------------------------------------------------------------------- 193 INTEGER :: ioptio, ios 194 ! 195 NAMELIST/namdyn_spg/ ln_dynspg_exp , ln_dynspg_ts,&196 & ln_bt_fw, ln_bt_av , ln_bt_auto,&197 & nn_baro , rn_bt_cmax, nn_bt_flt181 INTEGER :: ioptio, ios ! local integers 182 ! 183 NAMELIST/namdyn_spg/ ln_dynspg_exp , ln_dynspg_ts, & 184 & ln_bt_fw, ln_bt_av , ln_bt_auto , & 185 & nn_baro , rn_bt_cmax, nn_bt_flt 198 186 !!---------------------------------------------------------------------- 199 187 ! … … 202 190 REWIND( numnam_ref ) ! Namelist namdyn_spg in reference namelist : Free surface 203 191 READ ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901) 204 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist', lwp )205 192 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist', lwp ) 193 ! 206 194 REWIND( numnam_cfg ) ! Namelist namdyn_spg in configuration namelist : Free surface 207 195 READ ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 208 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp )196 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp ) 209 197 IF(lwm) WRITE ( numond, namdyn_spg ) 210 198 ! … … 216 204 WRITE(numout,*) ' Free surface with time splitting ln_dynspg_ts = ', ln_dynspg_ts 217 205 ENDIF 218 219 IF( ln_dynspg_ts ) THEN 220 CALL dyn_spg_ts_init( nit000 ) ! do it first, to set nn_baro, used to allocate some arrays later on 221 ! ! allocate dyn_spg arrays 222 IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays') 223 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 224 ENDIF 225 226 ! ! Control of surface pressure gradient scheme options 227 ioptio = 0 228 IF(ln_dynspg_exp) ioptio = ioptio + 1 229 IF(ln_dynspg_ts ) ioptio = ioptio + 1 230 ! 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 232 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 233 IF( ln_dynspg_ts .AND. ln_isfcav ) & 234 & CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' ) 235 ! 236 IF( ln_dynspg_exp) nspg = 0 237 IF( ln_dynspg_ts ) nspg = 1 206 ! ! Control of surface pressure gradient scheme options 207 ; nspg = np_NO ; ioptio = 0 208 IF( ln_dynspg_exp ) THEN ; nspg = np_EXP ; ioptio = ioptio + 1 ; ENDIF 209 IF( ln_dynspg_ts ) THEN ; nspg = np_TS ; ioptio = ioptio + 1 ; ENDIF 210 ! 211 IF( ioptio > 1 ) CALL ctl_stop( 'Choose only one surface pressure gradient scheme' ) 212 IF( ioptio == 0 ) CALL ctl_warn( 'NO surface pressure gradient trend in momentum Eqs.' ) 213 IF( ln_dynspg_exp .AND. ln_isfcav ) & 214 & CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' ) 238 215 ! 239 216 IF(lwp) THEN 240 217 WRITE(numout,*) 241 IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface' 242 IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme' 218 IF( nspg == np_EXP ) WRITE(numout,*) ' explicit free surface' 219 IF( nspg == np_TS ) WRITE(numout,*) ' free surface with time splitting scheme' 220 IF( nspg == np_NO ) WRITE(numout,*) ' No surface surface pressure gradient trend in momentum Eqs.' 221 ENDIF 222 ! 223 IF( nspg == np_TS ) THEN ! split-explicit scheme initialisation 224 CALL dyn_spg_ts_init ! do it first: set nn_baro used to allocate some arrays later on 225 IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' ) 226 IF( neuler/=0 .AND. ln_bt_fw ) CALL ts_rst( nit000, 'READ' ) 243 227 ENDIF 244 228 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r5930 r6140 2 2 !!====================================================================== 3 3 !! *** MODULE dynspg_exp *** 4 !! Ocean dynamics: surface pressure gradient trend 4 !! Ocean dynamics: surface pressure gradient trend, explicit scheme 5 5 !!====================================================================== 6 6 !! History : 2.0 ! 2005-11 (V. Garnier, G. Madec, L. Bessieres) Original code 7 7 !! 3.2 ! 2009-06 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 8 8 !!---------------------------------------------------------------------- 9 !!---------------------------------------------------------------------- 10 !! explicit free surface 9 11 10 !!---------------------------------------------------------------------- 12 11 !! dyn_spg_exp : update the momentum trend with the surface … … 26 25 USE timing ! Timing 27 26 28 29 27 IMPLICIT NONE 30 28 PRIVATE 31 29 32 PUBLIC dyn_spg_exp ! routine called by dyn _spg30 PUBLIC dyn_spg_exp ! routine called by dynspg.F90 33 31 34 32 !! * Substitutions 35 # include "domzgr_substitute.h90"36 33 # include "vectopt_loop_substitute.h90" 37 34 !!---------------------------------------------------------------------- … … 73 70 spgu(:,:) = 0._wp ; spgv(:,:) = 0._wp 74 71 ! 75 IF( lk_vvl .AND. lwp ) WRITE(numout,*) ' lk_vvl=T: spg is included in dynhpg'72 IF( .NOT.ln_linssh .AND. lwp ) WRITE(numout,*) ' non linear free surface: spg is included in dynhpg' 76 73 ENDIF 77 74 78 IF( .NOT. lk_vvl ) THEN !* fixed volume : add the surface pressure gradient trend75 IF( ln_linssh ) THEN !* linear free surface : add the surface pressure gradient trend 79 76 ! 80 77 DO jj = 2, jpjm1 ! now surface pressure gradient 81 78 DO ji = fs_2, fs_jpim1 ! vector opt. 82 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) /e1u(ji,jj)83 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) /e2v(ji,jj)79 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj) 80 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj) 84 81 END DO 85 82 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5930 r6140 1 1 MODULE dynspg_ts 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend, split-explicit scheme 2 5 !!====================================================================== 3 6 !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code … … 13 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 14 17 !!--------------------------------------------------------------------- 18 15 19 !!---------------------------------------------------------------------- 16 !! split explicit free surface17 !! ----------------------------------------------------------------------18 !! dyn_spg_ts : compute surface pressure gradient trend using a time-19 !! splitting scheme and add to the general trend20 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 21 !! dyn_spg_ts_init: initialisation of the time-splitting scheme 22 !! ts_wgt : set time-splitting weights for temporal averaging (or not) 23 !! ts_rst : read/write time-splitting fields in restart file 20 24 !!---------------------------------------------------------------------- 21 25 USE oce ! ocean dynamics and tracers 22 26 USE dom_oce ! ocean space and time domain 23 27 USE sbc_oce ! surface boundary condition: ocean 28 USE zdf_oce ! Bottom friction coefts 24 29 USE sbcisf ! ice shelf variable (fwfisf) 30 USE sbcapr ! surface boundary condition: atmospheric pressure 31 USE dynadv , ONLY: ln_dynadv_vec 25 32 USE phycst ! physical constants 26 33 USE dynvor ! vorticity term … … 30 37 USE sbctide ! tides 31 38 USE updtide ! tide potential 39 ! 40 USE in_out_manager ! I/O manager 32 41 USE lib_mpp ! distributed memory computing library 33 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 43 USE prtctl ! Print control 35 USE in_out_manager ! I/O manager36 44 USE iom ! IOM library 37 45 USE restart ! only for lrst_oce 38 USE zdf_oce ! Bottom friction coefts39 46 USE wrk_nemo ! Memory Allocation 40 47 USE timing ! Timing 41 USE sbcapr ! surface boundary condition: atmospheric pressure 42 USE dynadv, ONLY: ln_dynadv_vec 48 USE diatmb ! Top,middle,bottom output 43 49 #if defined key_agrif 44 50 USE agrif_opa_interp ! agrif … … 47 53 USE asminc ! Assimilation increment 48 54 #endif 55 49 56 50 57 IMPLICIT NONE … … 59 66 REAL(wp),SAVE :: rdtbt ! Barotropic time step 60 67 61 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 62 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 63 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 64 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 69 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff/h at F points 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) 73 74 !! Time filtered arrays at baroclinic time step: 75 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 68 76 69 77 !! * Substitutions 70 # include "domzgr_substitute.h90"71 78 # include "vectopt_loop_substitute.h90" 72 79 !!---------------------------------------------------------------------- … … 81 88 !! *** routine dyn_spg_ts_alloc *** 82 89 !!---------------------------------------------------------------------- 83 INTEGER :: ierr( 4)90 INTEGER :: ierr(3) 84 91 !!---------------------------------------------------------------------- 85 92 ierr(:) = 0 86 87 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 88 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 89 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 90 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT= ierr(1) ) 91 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 93 94 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 96 97 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 98 #if defined key_agrif 99 & ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , & 100 #endif 101 & STAT= ierr(4)) 102 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 93 ! 94 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 95 ! 96 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 97 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 98 ! 99 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) 100 ! 101 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 102 ! 105 103 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 104 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') … … 112 110 !!---------------------------------------------------------------------- 113 111 !! 114 !! ** Purpose : 115 !! -Compute the now trend due to the explicit time stepping116 !! of the quasi-linear barotropic system.112 !! ** Purpose : - Compute the now trend due to the explicit time stepping 113 !! of the quasi-linear barotropic system, and add it to the 114 !! general momentum trend. 117 115 !! 118 !! ** Method : 116 !! ** Method : - split-explicit schem (time splitting) : 119 117 !! Barotropic variables are advanced from internal time steps 120 118 !! "n" to "n+1" if ln_bt_fw=T … … 130 128 !! continuity equation taken at the baroclinic time steps. This 131 129 !! ensures tracers conservation. 132 !! - Update 3d trend (ua, va)with barotropic component.130 !! - (ua, va) momentum trend updated with barotropic component. 133 131 !! 134 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 135 !! The regional oceanic modeling system (ROMS): 136 !! a split-explicit, free-surface, 137 !! topography-following-coordinate oceanic model. 138 !! Ocean Modelling, 9, 347-404. 132 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 139 133 !!--------------------------------------------------------------------- 140 !141 134 INTEGER, INTENT(in) :: kt ! ocean time-step index 142 135 ! … … 145 138 INTEGER :: ji, jj, jk, jn ! dummy loop indices 146 139 INTEGER :: ikbu, ikbv, noffset ! local integers 140 INTEGER :: iktu, iktv ! local integers 141 REAL(wp) :: zmdi 147 142 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 143 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - … … 163 158 ! 164 159 ! !* Allocate temporary arrays 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 170 CALL wrk_alloc( jpi, jpj, zhf ) 171 ! 160 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 161 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 162 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 163 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 164 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 165 CALL wrk_alloc( jpi,jpj, zhf ) 166 ! 167 zmdi=1.e+20 ! missing data indicator for masking 172 168 ! !* Local constant initialization 173 169 z1_12 = 1._wp / 12._wp … … 176 172 z1_2 = 0.5_wp 177 173 zraur = 1._wp / rau0 178 ! 179 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 180 z2dt_bf = rdt 181 ELSE 182 z2dt_bf = 2.0_wp * rdt 174 ! ! reciprocal of baroclinic time step 175 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 176 ELSE ; z2dt_bf = 2.0_wp * rdt 183 177 ENDIF 184 178 z1_2dt_b = 1.0_wp / z2dt_bf 185 179 ! 186 ll_init = ln_bt_av! if no time averaging, then no specific restart180 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 181 ll_fw_start = .FALSE. 188 ! 189 ! time offset in steps for bdy data update 190 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 182 ! ! time offset in steps for bdy data update 183 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 184 ELSE ; noffset = 0 185 ENDIF 191 186 ! 192 187 IF( kt == nit000 ) THEN !* initialisation … … 197 192 IF(lwp) WRITE(numout,*) 198 193 ! 199 IF (neuler==0)ll_init=.TRUE.200 ! 201 IF (ln_bt_fw.OR.(neuler==0)) THEN202 ll_fw_start=.TRUE.203 noffset= 0194 IF( neuler == 0 ) ll_init=.TRUE. 195 ! 196 IF( ln_bt_fw .OR. neuler == 0 ) THEN 197 ll_fw_start =.TRUE. 198 noffset = 0 204 199 ELSE 205 ll_fw_start=.FALSE.200 ll_fw_start =.FALSE. 206 201 ENDIF 207 202 ! 208 203 ! Set averaging weights and cycle length: 209 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 210 ! 204 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 211 205 ! 212 206 ENDIF … … 219 213 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 220 214 ! 221 IF ( kt == nit000 .OR. lk_vvl) THEN222 IF ( ln_dynvor_een ) THEN!== EEN scheme ==!215 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 216 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 223 217 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 224 218 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 225 219 DO jj = 1, jpjm1 226 220 DO ji = 1, jpim1 227 zwz(ji,jj) = ( ht (ji ,jj+1) + ht(ji+1,jj+1) + &228 & ht (ji ,jj ) + ht(ji+1,jj ) ) / 4._wp221 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 222 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 229 223 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 224 END DO … … 233 227 DO jj = 1, jpjm1 234 228 DO ji = 1, jpim1 235 zwz(ji,jj) = ( ht (ji ,jj+1) + ht(ji+1,jj+1) + &236 & ht (ji ,jj ) + ht(ji+1,jj ) ) &229 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 230 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 237 231 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 232 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) … … 276 270 DO jk = 1, jpkm1 277 271 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)272 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 279 273 END DO 280 274 END DO … … 308 302 ! 309 303 DO jk = 1, jpkm1 310 zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)311 zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)304 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 305 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 312 306 END DO 313 307 ! 314 zu_frc(:,:) = zu_frc(:,:) * hur(:,:)315 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)308 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 309 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 316 310 ! 317 311 ! … … 327 321 ! !* barotropic Coriolis trends (vorticity scheme dependent) 328 322 ! ! -------------------------------------------------------- 329 zwx(:,:) = un_b(:,:) * hu (:,:) * e2u(:,:) ! now fluxes330 zwy(:,:) = vn_b(:,:) * hv (:,:) * e1v(:,:)323 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 324 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 331 325 ! 332 326 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 373 367 ! !* Right-Hand-Side of the barotropic momentum equation 374 368 ! ! ---------------------------------------------------- 375 IF( lk_vvl ) THEN! Variable volume : remove surface pressure gradient369 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 376 370 DO jj = 2, jpjm1 377 371 DO ji = fs_2, fs_jpim1 ! vector opt. … … 384 378 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 385 379 DO ji = fs_2, fs_jpim1 386 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)387 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)380 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 381 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 388 382 END DO 389 383 END DO … … 411 405 ! 412 406 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 413 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 414 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 407 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 408 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 409 ! 410 ! ! Add top stress contribution from baroclinic velocities: 411 IF (ln_bt_fw) THEN 412 DO jj = 2, jpjm1 413 DO ji = fs_2, fs_jpim1 ! vector opt. 414 iktu = miku(ji,jj) 415 iktv = mikv(ji,jj) 416 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 417 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 418 END DO 419 END DO 420 ELSE 421 DO jj = 2, jpjm1 422 DO ji = fs_2, fs_jpim1 ! vector opt. 423 iktu = miku(ji,jj) 424 iktv = mikv(ji,jj) 425 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 426 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 427 END DO 428 END DO 429 ENDIF 430 ! 431 ! Note that the "unclipped" top friction parameter is used even with explicit drag 432 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 433 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 415 434 ! 416 435 IF (ln_bt_fw) THEN ! Add wind forcing 417 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)418 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)436 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 437 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 419 438 ELSE 420 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)421 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)439 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 440 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 422 441 ENDIF 423 442 ! … … 484 503 ! 485 504 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 sshn_e(:,:) = sshn(:,:)487 un_e (:,:) = un_b(:,:)488 vn_e (:,:) = vn_b(:,:)489 ! 490 hu_e (:,:) = hu(:,:)491 hv_e (:,:) = hv(:,:)492 hur_e (:,:) = hur(:,:)493 hvr_e (:,:) = hvr(:,:)505 sshn_e(:,:) = sshn(:,:) 506 un_e (:,:) = un_b(:,:) 507 vn_e (:,:) = vn_b(:,:) 508 ! 509 hu_e (:,:) = hu_n(:,:) 510 hv_e (:,:) = hv_n(:,:) 511 hur_e (:,:) = r1_hu_n(:,:) 512 hvr_e (:,:) = r1_hv_n(:,:) 494 513 ELSE ! CENTRED integration: start from BEFORE fields 495 sshn_e(:,:) = sshb(:,:)496 un_e (:,:) = ub_b(:,:)497 vn_e (:,:) = vb_b(:,:)498 ! 499 hu_e (:,:) = hu_b(:,:)500 hv_e (:,:) = hv_b(:,:)501 hur_e (:,:) = hur_b(:,:)502 hvr_e (:,:) = hvr_b(:,:)514 sshn_e(:,:) = sshb(:,:) 515 un_e (:,:) = ub_b(:,:) 516 vn_e (:,:) = vb_b(:,:) 517 ! 518 hu_e (:,:) = hu_b(:,:) 519 hv_e (:,:) = hv_b(:,:) 520 hur_e (:,:) = r1_hu_b(:,:) 521 hvr_e (:,:) = r1_hv_b(:,:) 503 522 ENDIF 504 523 ! … … 518 537 ! Update only tidal forcing at open boundaries 519 538 #if defined key_tide 520 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1))521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset)539 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 540 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 522 541 #endif 523 542 ! … … 537 556 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 557 539 IF( lk_vvl ) THEN!* Update ocean depth (variable volume case only)558 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 540 559 ! ! ------------------ 541 560 ! Extrapolate Sea Level at step jit+0.5: … … 544 563 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 545 564 DO ji = 2, fs_jpim1 ! Vector opt. 546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &565 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 547 566 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 567 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &568 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 550 569 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 570 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) … … 557 576 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 558 577 ELSE 559 zhup2_e (:,:) = hu (:,:)560 zhvp2_e (:,:) = hv (:,:)578 zhup2_e (:,:) = hu_n(:,:) 579 zhvp2_e (:,:) = hv_n(:,:) 561 580 ENDIF 562 581 ! !* after ssh … … 569 588 ! 570 589 #if defined key_agrif 571 ! Set fluxes during predictor step to ensure 572 ! volume conservation 573 IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN 590 ! Set fluxes during predictor step to ensure volume conservation 591 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 574 592 IF((nbondi == -1).OR.(nbondi == 2)) THEN 575 593 DO jj=1,jpj … … 607 625 END DO 608 626 END DO 609 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)627 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 610 628 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 611 629 612 630 #if defined key_bdy 613 ! Duplicate sea level across open boundaries (this is only cosmetic if l k_vvl=.false.)614 IF (lk_bdy)CALL bdy_ssh( ssha_e )631 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 632 IF( lk_bdy ) CALL bdy_ssh( ssha_e ) 615 633 #endif 616 634 #if defined key_agrif 617 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )635 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 618 636 #endif 619 637 ! 620 638 ! Sea Surface Height at u-,v-points (vvl case only) 621 IF ( lk_vvl) THEN639 IF( .NOT.ln_linssh ) THEN 622 640 DO jj = 2, jpjm1 623 641 DO ji = 2, jpim1 ! NO Vector Opt. 624 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj)&625 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &626 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )627 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj)&628 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &629 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )642 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 643 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 644 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 645 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 646 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 647 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 630 648 END DO 631 649 END DO … … 651 669 za3=0.013_wp ! za3 = eps 652 670 ENDIF 653 671 ! 654 672 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 655 673 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 656 657 674 ! 658 675 ! Compute associated depths at U and V points: 659 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN676 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 660 677 ! 661 678 DO jj = 2, jpjm1 662 679 DO ji = 2, jpim1 663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) &680 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 664 681 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 682 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) &683 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 667 684 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 685 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) … … 674 691 ! 675 692 ! Add Coriolis trend: 676 ! zwz array below or triads normally depend on sea level with key_vvland should be updated693 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 677 694 ! at each time step. We however keep them constant here for optimization. 678 695 ! Recall that zwx and zwy arrays hold fluxes at this stage: … … 680 697 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 681 698 ! 682 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN 699 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 683 700 DO jj = 2, jpjm1 684 701 DO ji = fs_2, fs_jpim1 ! vector opt. … … 692 709 END DO 693 710 ! 694 ELSEIF ( ln_dynvor_ens ) THEN 711 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 695 712 DO jj = 2, jpjm1 696 713 DO ji = fs_2, fs_jpim1 ! vector opt. … … 704 721 END DO 705 722 ! 706 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!723 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 707 724 DO jj = 2, jpjm1 708 725 DO ji = fs_2, fs_jpim1 ! vector opt. … … 736 753 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 754 ! 755 ! Add top stresses: 756 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 757 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 758 ! 738 759 ! Surface pressure trend: 739 760 DO jj = 2, jpjm1 … … 748 769 ! 749 770 ! Set next velocities: 750 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN !Vector form771 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 751 772 DO jj = 2, jpjm1 752 773 DO ji = fs_2, fs_jpim1 ! vector opt. … … 755 776 & + zu_trd(ji,jj) & 756 777 & + zu_frc(ji,jj) ) & 757 & ) * umask(ji,jj,1)778 & ) * ssumask(ji,jj) 758 779 759 780 va_e(ji,jj) = ( vn_e(ji,jj) & … … 761 782 & + zv_trd(ji,jj) & 762 783 & + zv_frc(ji,jj) ) & 763 & ) * vmask(ji,jj,1)764 END DO 765 END DO 766 767 ELSE !Flux form784 & ) * ssvmask(ji,jj) 785 END DO 786 END DO 787 ! 788 ELSE !* Flux form 768 789 DO jj = 2, jpjm1 769 790 DO ji = fs_2, fs_jpim1 ! vector opt. 770 771 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 791 zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 792 zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 773 793 774 794 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 795 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 796 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 777 & + hu(ji,jj) * zu_frc(ji,jj) ) &797 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 778 798 & ) * zhura 779 799 … … 781 801 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 802 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 783 & + hv(ji,jj) * zv_frc(ji,jj) ) &803 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 784 804 & ) * zhvra 785 805 END DO … … 787 807 ENDIF 788 808 ! 789 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 790 ! ! ---------------------------------------------- 809 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 791 810 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 792 811 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 793 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )794 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )812 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 813 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 795 814 ! 796 815 ENDIF 797 ! !* domain lateral boundary 798 ! ! ----------------------- 799 ! 816 ! !* domain lateral boundary 800 817 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 801 818 ! 802 819 #if defined key_bdy 803 804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )820 ! ! open boundaries 821 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 822 #endif 806 823 #if defined key_agrif … … 824 841 ! ! ---------------------- 825 842 za1 = wgtbtp1(jn) 826 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities843 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 827 844 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 828 845 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 829 ELSE 846 ELSE ! Sum transports 830 847 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 831 848 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) … … 843 860 zwx(:,:) = un_adv(:,:) 844 861 zwy(:,:) = vn_adv(:,:) 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN846 un_adv(:,:) = zwx(:,:) *hur(:,:)847 vn_adv(:,:) = zwy(:,:) *hvr(:,:)862 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 863 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 864 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 848 865 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * hur(:,:)850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * hvr(:,:)866 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 867 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 851 868 END IF 852 869 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation870 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 854 871 ub2_b(:,:) = zwx(:,:) 855 872 vb2_b(:,:) = zwy(:,:) … … 857 874 ! 858 875 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN876 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 860 877 DO jk=1,jpkm1 861 878 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b … … 877 894 ! 878 895 DO jk=1,jpkm1 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b880 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b896 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 897 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 881 898 END DO 882 899 ! Save barotropic velocities not transport: 883 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )884 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )900 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 901 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 885 902 ENDIF 886 903 ! 887 904 DO jk = 1, jpkm1 888 905 ! Correct velocities: 889 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) *umask(:,:,jk)890 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) *vmask(:,:,jk)906 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 907 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 891 908 ! 892 909 END DO 910 ! 911 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 912 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current 893 913 ! 894 914 #if defined key_agrif … … 896 916 ! (used to update coarse grid transports at next time step) 897 917 ! 898 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw)) THEN899 IF ( Agrif_NbStepint().EQ.0 ) THEN900 ub2_i_b(:,:) = 0. e0901 vb2_i_b(:,:) = 0. e0918 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 919 IF( Agrif_NbStepint() == 0 ) THEN 920 ub2_i_b(:,:) = 0._wp 921 vb2_i_b(:,:) = 0._wp 902 922 END IF 903 923 ! … … 906 926 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 907 927 ENDIF 908 !909 !910 928 #endif 911 !912 929 ! !* write time-spliting arrays in the restart 913 IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' ) 914 ! 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 916 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 917 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 920 CALL wrk_dealloc( jpi, jpj, zhf ) 921 ! 930 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 931 ! 932 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 933 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 934 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 935 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 936 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 937 CALL wrk_dealloc( jpi,jpj, zhf ) 938 ! 939 IF ( ln_diatmb ) THEN 940 CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity 941 CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity 942 ENDIF 922 943 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 923 944 ! 924 945 END SUBROUTINE dyn_spg_ts 946 925 947 926 948 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) … … 1001 1023 END SUBROUTINE ts_wgt 1002 1024 1025 1003 1026 SUBROUTINE ts_rst( kt, cdrw ) 1004 1027 !!--------------------------------------------------------------------- … … 1054 1077 END SUBROUTINE ts_rst 1055 1078 1056 SUBROUTINE dyn_spg_ts_init( kt ) 1079 1080 SUBROUTINE dyn_spg_ts_init 1057 1081 !!--------------------------------------------------------------------- 1058 1082 !! *** ROUTINE dyn_spg_ts_init *** … … 1060 1084 !! ** Purpose : Set time splitting options 1061 1085 !!---------------------------------------------------------------------- 1062 INTEGER , INTENT(in) :: kt ! ocean time-step 1063 ! 1064 INTEGER :: ji ,jj 1065 REAL(wp) :: zxr2, zyr2, zcmax 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1067 !! 1086 INTEGER :: ji ,jj ! dummy loop indices 1087 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1088 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1068 1089 !!---------------------------------------------------------------------- 1069 1090 ! 1070 1091 ! Max courant number for ext. grav. waves 1071 1092 ! 1072 CALL wrk_alloc( jpi, jpj,zcu )1093 CALL wrk_alloc( jpi,jpj, zcu ) 1073 1094 ! 1074 1095 DO jj = 1, jpj … … 1084 1105 1085 1106 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1086 IF (ln_bt_auto)nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1107 IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1087 1108 1088 1109 rdtbt = rdt / REAL( nn_baro , wp ) … … 1114 1135 #if defined key_agrif 1115 1136 ! Restrict the use of Agrif to the forward case only 1116 IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root()))CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )1137 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1117 1138 #endif 1118 1139 ! 1119 1140 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1120 1141 SELECT CASE ( nn_bt_flt ) 1121 CASE( 0 ); IF(lwp) WRITE(numout,*) ' Dirac'1122 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'1123 CASE( 2 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro'1124 1142 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1143 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1144 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1145 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1125 1146 END SELECT 1126 1147 ! … … 1130 1151 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1131 1152 ! 1132 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN1153 IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 1133 1154 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1134 1155 ENDIF 1135 IF 1156 IF( zcmax>0.9_wp ) THEN 1136 1157 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1137 1158 ENDIF 1138 1159 ! 1139 CALL wrk_dealloc( jpi, jpj,zcu )1160 CALL wrk_dealloc( jpi,jpj, zcu ) 1140 1161 ! 1141 1162 END SUBROUTINE dyn_spg_ts_init -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5930 r6140 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 USE c1d ! 1D vertical configuration35 34 ! 36 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 76 75 77 76 !! * Substitutions 78 # include "domzgr_substitute.h90"79 77 # include "vectopt_loop_substitute.h90" 80 78 !!---------------------------------------------------------------------- … … 285 283 286 284 IF( ln_sco ) THEN 287 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)288 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)289 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)285 zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 286 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 287 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 290 288 ELSE 291 289 zwx(:,:) = e2u(:,:) * un(:,:,jk) … … 405 403 ! 406 404 IF( ln_sco ) THEN !== horizontal fluxes ==! 407 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)408 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)409 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)405 zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 406 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 407 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 410 408 ELSE 411 409 zwx(:,:) = e2u(:,:) * un(:,:,jk) … … 415 413 DO jj = 2, jpjm1 416 414 DO ji = fs_2, fs_jpim1 ! vector opt. 417 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1)&418 & + zwy(ji ,jj ) + zwy(ji+1,jj ))419 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1)&420 & + zwx(ji ,jj ) + zwx(ji ,jj+1))415 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 416 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 417 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 418 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 421 419 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 422 420 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 483 481 DO jj = 1, jpjm1 484 482 DO ji = 1, fs_jpim1 ! vector opt. 485 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &486 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk))487 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4. 0_wp / ze3488 ELSE ; z1_e3f(ji,jj) = 0. 0_wp483 ze3 = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 484 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 485 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3 486 ELSE ; z1_e3f(ji,jj) = 0._wp 489 487 ENDIF 490 488 END DO … … 493 491 DO jj = 1, jpjm1 494 492 DO ji = 1, fs_jpim1 ! vector opt. 495 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) &496 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk))497 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &498 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk))493 ze3 = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 494 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 495 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 496 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 499 497 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3 500 ELSE ; z1_e3f(ji,jj) = 0. 0_wp498 ELSE ; z1_e3f(ji,jj) = 0._wp 501 499 ENDIF 502 500 END DO … … 559 557 ! 560 558 ! !== horizontal fluxes ==! 561 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)562 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)559 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 560 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 563 561 564 562 ! !== compute and add the vorticity term trend =! … … 634 632 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 635 633 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 636 WRITE(numout,*) ' masked (= 1) or unmasked(=0) vorticity ln_dynvor_msk = ', ln_dynvor_msk634 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 637 635 ENDIF 638 636 … … 663 661 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 664 662 ! 665 IF( ( ioptio /= 1).AND.( .NOT.lk_c1d )) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )663 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 666 664 ! 667 665 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5836 r6140 32 32 33 33 !! * Substitutions 34 # include "domzgr_substitute.h90"35 34 # include "vectopt_loop_substitute.h90" 36 35 !!---------------------------------------------------------------------- … … 90 89 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 91 90 DO ji = fs_2, fs_jpim1 ! vector opt. 92 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) -un(ji,jj,jk) )93 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) -vn(ji,jj,jk) )91 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) 92 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) 94 93 END DO 95 94 END DO … … 121 120 DO ji = fs_2, fs_jpim1 ! vector opt. 122 121 ! ! vertical momentum advective trends 123 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk))124 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))122 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 123 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 125 124 ! ! add the trends to the general momentum trends 126 125 ua(ji,jj,jk) = ua(ji,jj,jk) + zua … … 252 251 DO ji = fs_2, fs_jpim1 ! vector opt. 253 252 ! ! vertical momentum advective trends 254 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk))255 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))253 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 254 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 256 255 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 257 256 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r6132 r6140 34 34 35 35 INTEGER :: nzdf = 0 ! type vertical diffusion algorithm used, defined from ln_zdf... namlist logicals 36 REAL(wp) :: r2dt ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=037 36 38 37 !! * Substitutions 39 # include "domzgr_substitute.h90"40 # include "zdfddm_substitute.h90"41 38 # include "vectopt_loop_substitute.h90" 42 39 !!---------------------------------------------------------------------- … … 63 60 ! 64 61 ! ! set time step 65 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt ra(restart with Euler time stepping)66 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt tra(leapfrog)62 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 63 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 67 64 ENDIF 68 65 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r5930 r6140 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 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5930 r6140 2 2 !!====================================================================== 3 3 !! *** MODULE dynzdf_imp *** 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend, implicit scheme 5 5 !!====================================================================== 6 6 !! History : OPA ! 1990-10 (B. Blanke) Original code … … 12 12 13 13 !!---------------------------------------------------------------------- 14 !! dyn_zdf_imp : update the momentum trend with the vertical diffusion using a implicit time-stepping 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE domvvl ! variable volume 19 USE sbc_oce ! surface boundary condition: ocean 20 USE zdf_oce ! ocean vertical physics 21 USE phycst ! physical constants 22 USE in_out_manager ! I/O manager 23 USE lib_mpp ! MPP library 24 USE zdfbfr ! Bottom friction setup 25 USE wrk_nemo ! Memory Allocation 26 USE timing ! Timing 27 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 14 !! dyn_zdf_imp : compute the vertical diffusion using a implicit scheme 15 !! together with the Leap-Frog time integration. 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE phycst ! physical constants 19 USE dom_oce ! ocean space and time domain 20 USE domvvl ! variable volume 21 USE sbc_oce ! surface boundary condition: ocean 22 USE dynadv , ONLY: ln_dynadv_vec ! Momentum advection form 23 USE zdf_oce ! ocean vertical physics 24 USE zdfbfr ! Bottom friction setup 25 ! 26 USE in_out_manager ! I/O manager 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 28 30 29 31 IMPLICIT NONE … … 32 34 PUBLIC dyn_zdf_imp ! called by step.F90 33 35 34 REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0otherwise36 REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 35 37 36 38 !! * Substitutions 37 # include "domzgr_substitute.h90"38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- … … 49 50 !! 50 51 !! ** Purpose : Compute the trend due to the vert. momentum diffusion 51 !! and the surface forcing, and add it to the general trend of52 !! the momentum equations.52 !! together with the Leap-Frog time stepping using an 53 !! implicit scheme. 53 54 !! 54 !! ** Method : The vertical momentum mixing trend is given by : 55 !! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 56 !! backward time stepping 57 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 58 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F) 59 !! Add this trend to the general trend ua : 60 !! ua = ua + dz( avmu dz(u) ) 55 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! ua = ub + 2*dt * ua vector form or linear free surf. 57 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 58 !! - update the after velocity with the implicit vertical mixing. 59 !! This requires to solver the following system: 60 !! ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ] 61 !! with the following surface/top/bottom boundary condition: 62 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 63 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfbfr.F) 61 64 !! 62 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.65 !! ** Action : (ua,va) after velocity 63 66 !!--------------------------------------------------------------------- 64 67 INTEGER , INTENT(in) :: kt ! ocean time-step index 65 68 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 66 ! !67 INTEGER :: ji, jj, jk ! dummy loop indices68 INTEGER :: ikbu, ikbv ! local integers69 REAL(wp) :: z 1_p2dt, zcoef, zzwi, zzws, zrhs! local scalars70 REAL(wp) :: z e3ua, ze3va69 ! 70 INTEGER :: ji, jj, jk ! dummy loop indices 71 INTEGER :: ikbu, ikbv ! local integers 72 REAL(wp) :: zzwi, ze3ua ! local scalars 73 REAL(wp) :: zzws, ze3va ! - - 71 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws 72 75 !!---------------------------------------------------------------------- … … 81 84 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 82 85 ! 83 I F( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator84 ELSE ; r_vvl = 0._wp86 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 87 ELSE ; r_vvl = 1._wp 85 88 ENDIF 86 89 ENDIF 87 88 ! 1. Time step dynamics 89 ! --------------------- 90 91 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 90 ! 91 ! !== Time step dynamics ==! 92 ! 93 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 92 94 DO jk = 1, jpkm1 93 95 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 94 96 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 95 97 END DO 96 ELSE 98 ELSE ! applied on thickness weighted velocity 97 99 DO jk = 1, jpkm1 98 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 99 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 100 & / fse3u_a(:,:,jk) * umask(:,:,jk) 101 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 102 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 103 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 104 END DO 105 ENDIF 106 107 ! 2. Apply semi-implicit bottom friction 108 ! -------------------------------------- 100 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 101 & + p2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 102 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 103 & + p2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 104 END DO 105 ENDIF 106 ! 107 ! !== Apply semi-implicit bottom friction ==! 108 ! 109 109 ! Only needed for semi-implicit bottom friction setup. The explicit 110 110 ! bottom friction has been included in "u(v)a" which act as the R.H.S … … 116 116 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 117 117 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 118 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)119 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)118 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1) 119 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1) 120 120 END DO 121 121 END DO … … 125 125 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 126 126 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 127 IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)128 IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)127 IF( ikbu >= 2 ) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu) 128 IF( ikbv >= 2 ) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv) 129 129 END DO 130 130 END DO 131 131 END IF 132 132 ENDIF 133 133 ! 134 134 ! With split-explicit free surface, barotropic stress is treated explicitly 135 135 ! Update velocities at the bottom. 136 136 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 137 137 ! not lead to the effective stress seen over the whole barotropic loop. 138 IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN 139 ! remove barotropic velocities: 140 DO jk = 1, jpkm1 141 ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 142 va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 143 END DO 144 ! Add bottom/top stress due to barotropic component only: 145 DO jj = 2, jpjm1 138 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 139 IF( ln_bfrimp .AND. ln_dynspg_ts ) THEN 140 DO jk = 1, jpkm1 ! remove barotropic velocities 141 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 142 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 143 END DO 144 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 146 145 DO ji = fs_2, fs_jpim1 ! vector opt. 147 146 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 148 147 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 149 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)150 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)148 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 149 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 151 150 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 152 151 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 153 152 END DO 154 153 END DO 155 IF ( ln_isfcav ) THEN154 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) 156 155 DO jj = 2, jpjm1 157 156 DO ji = fs_2, fs_jpim1 ! vector opt. 158 157 ikbu = miku(ji,jj) ! top ocean level at u- and v-points 159 158 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 160 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)161 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)159 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 160 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 162 161 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 163 162 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va … … 166 165 END IF 167 166 ENDIF 168 169 ! 3. Vertical diffusion on u170 ! ---------------------------167 ! 168 ! !== Vertical diffusion on u ==! 169 ! 171 170 ! Matrix and second member construction 172 171 ! bottom boundary condition: both zwi and zws must be masked as avmu can take … … 176 175 DO jj = 2, jpjm1 177 176 DO ji = fs_2, fs_jpim1 ! vector opt. 178 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point 179 zcoef = - p2dt / ze3ua 180 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 181 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 182 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 183 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 177 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 178 zzwi = - p2dt * avmu(ji,jj,jk ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) 179 zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 180 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 181 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 184 182 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 185 183 END DO … … 216 214 END DO 217 215 ! 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1)216 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 221 219 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 222 220 & / ( ze3ua * rau0 ) * umask(ji,jj,1) … … 226 224 DO jj = 2, jpjm1 227 225 DO ji = fs_2, fs_jpim1 228 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 229 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 230 END DO 231 END DO 232 END DO 233 ! 234 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 226 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 227 END DO 228 END DO 229 END DO 230 ! 231 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 235 232 DO ji = fs_2, fs_jpim1 ! vector opt. 236 233 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 244 241 END DO 245 242 END DO 246 247 ! 4. Vertical diffusion on v248 ! ---------------------------243 ! 244 ! !== Vertical diffusion on v ==! 245 ! 249 246 ! Matrix and second member construction 250 247 ! bottom boundary condition: both zwi and zws must be masked as avmv can take … … 254 251 DO jj = 2, jpjm1 255 252 DO ji = fs_2, fs_jpim1 ! vector opt. 256 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk) + r_vvl * fse3v_a(ji,jj,jk) ! after scale factor at T-point 257 zcoef = - p2dt / ze3va 258 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 259 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk) 260 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 261 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 253 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 254 zzwi = - p2dt * avmv (ji,jj,jk ) / ( ze3va * e3vw_n(ji,jj,jk ) ) 255 zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 256 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 257 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 262 258 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 263 259 END DO … … 294 290 END DO 295 291 ! 296 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 292 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 297 293 DO ji = fs_2, fs_jpim1 ! vector opt. 298 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1)294 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 299 295 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 300 296 & / ( ze3va * rau0 ) … … 304 300 DO jj = 2, jpjm1 305 301 DO ji = fs_2, fs_jpim1 ! vector opt. 306 zrhs = va(ji,jj,jk) ! zrhs=right hand side 307 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 308 END DO 309 END DO 310 END DO 311 ! 312 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 302 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 303 END DO 304 END DO 305 END DO 306 ! 307 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 313 308 DO ji = fs_2, fs_jpim1 ! vector opt. 314 309 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 322 317 END DO 323 318 END DO 324 319 325 320 ! J. Chanut: Lines below are useless ? 326 !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 321 !! restore bottom layer avmu(v) 322 !!gm I almost sure it is !!!! 327 323 IF( ln_bfrimp ) THEN 328 324 DO jj = 2, jpjm1 … … 330 326 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 331 327 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 332 avmu(ji,jj,ikbu+1) = 0. e0333 avmv(ji,jj,ikbv+1) = 0. e0328 avmu(ji,jj,ikbu+1) = 0._wp 329 avmv(ji,jj,ikbv+1) = 0._wp 334 330 END DO 335 331 END DO … … 339 335 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 340 336 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 341 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0342 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0337 IF( ikbu > 1 ) avmu(ji,jj,ikbu) = 0._wp 338 IF( ikbv > 1 ) avmv(ji,jj,ikbv) = 0._wp 343 339 END DO 344 340 END DO 345 END 346 ENDIF 347 ! 348 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)349 ! 350 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')341 ENDIF 342 ENDIF 343 ! 344 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 345 ! 346 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp') 351 347 ! 352 348 END SUBROUTINE dyn_zdf_imp -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5930 r6140 12 12 13 13 !!---------------------------------------------------------------------- 14 !! ssh_nxt 15 !! ssh_swp 16 !! wzv 17 !!---------------------------------------------------------------------- 18 USE oce 19 USE dom_oce 20 USE sbc_oce 21 USE domvvl 22 USE divhor 23 USE phycst 24 USE bdy_oce 25 USE bdy_par 26 USE bdydyn2d 14 !! ssh_nxt : after ssh 15 !! ssh_swp : filter ans swap the ssh arrays 16 !! wzv : compute now vertical velocity 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and tracers variables 19 USE dom_oce ! ocean space and time domain variables 20 USE sbc_oce ! surface boundary condition: ocean 21 USE domvvl ! Variable volume 22 USE divhor ! horizontal divergence 23 USE phycst ! physical constants 24 USE bdy_oce ! 25 USE bdy_par ! 26 USE bdydyn2d ! bdy_ssh routine 27 27 #if defined key_agrif 28 28 USE agrif_opa_interp 29 29 #endif 30 30 #if defined key_asminc 31 USE asminc ! Assimilation increment 32 #endif 33 USE in_out_manager ! I/O manager 34 USE restart ! only for lrst_oce 35 USE prtctl ! Print control 36 USE lbclnk ! ocean lateral boundary condition (or mpp link) 37 USE lib_mpp ! MPP library 38 USE wrk_nemo ! Memory Allocation 39 USE timing ! Timing 31 USE asminc ! Assimilation increment 32 #endif 33 ! 34 USE in_out_manager ! I/O manager 35 USE restart ! only for lrst_oce 36 USE prtctl ! Print control 37 USE lbclnk ! ocean lateral boundary condition (or mpp link) 38 USE lib_mpp ! MPP library 39 USE wrk_nemo ! Memory Allocation 40 USE timing ! Timing 40 41 41 42 IMPLICIT NONE … … 47 48 48 49 !! * Substitutions 49 # include "domzgr_substitute.h90"50 50 # include "vectopt_loop_substitute.h90" 51 51 !!---------------------------------------------------------------------- … … 97 97 zhdiv(:,:) = 0._wp 98 98 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 99 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)99 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 100 100 END DO 101 101 ! ! Sea surface elevation time stepping … … 194 194 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 195 195 ! computation of w 196 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)&197 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )) * tmask(:,:,jk)196 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) & 197 & + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 198 198 END DO 199 199 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 … … 202 202 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 203 203 ! computation of w 204 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk)&205 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )) * tmask(:,:,jk)204 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) & 205 & + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 206 206 END DO 207 207 ENDIF … … 240 240 !!---------------------------------------------------------------------- 241 241 INTEGER, INTENT(in) :: kt ! ocean time-step index 242 ! 243 REAL(wp) :: zcoef ! local scalar 242 244 !!---------------------------------------------------------------------- 243 245 ! … … 249 251 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 250 252 ENDIF 251 252 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN253 !** Euler time-stepping: no filter254 sshb(:,:) = sshn(:,:) ! before <-- now255 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)253 ! !== Euler time-stepping: no filter, just swap ==! 254 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. & 255 & ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 256 sshb(:,:) = sshn(:,:) ! before <-- now 257 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 256 258 ! 257 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 258 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 259 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) & 260 & - rnf_b(:,:) + rnf(:,:) & 261 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 262 sshn(:,:) = ssha(:,:) ! now <-- after 259 ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==! 260 ! ! before <-- now filtered 261 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 262 IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed 263 zcoef = atfp * rdt * r1_rau0 264 sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) & 265 & - rnf_b(:,:) + rnf (:,:) & 266 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 267 ENDIF 268 sshn(:,:) = ssha(:,:) ! now <-- after 263 269 ENDIF 264 270 ! 265 271 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 266 272 ! 267 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp')273 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp') 268 274 ! 269 275 END SUBROUTINE ssh_swp
Note: See TracChangeset
for help on using the changeset viewer.