- Timestamp:
- 2011-07-19T18:35:40+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r2779 r2807 24 24 PRIVATE 25 25 26 PUBLIC dom_vvl ! called by domain.F9027 PUBLIC dom_vvl_ alloc ! called by nemogcm.F9028 29 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ee_t, ee_u, ee_v, ee_f !: ??? 30 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: ???26 PUBLIC dom_vvl ! called by domain.F90 27 PUBLIC dom_vvl_2 ! called by domain.F90 28 PUBLIC dom_vvl_alloc ! called by nemogcm.F90 29 30 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: 1/H_0 at t-,u-,v-,f-points 31 31 32 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra … … 49 49 ! 50 50 ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) , & 51 & ee_t(jpi,jpj) , ee_u(jpi,jpj) , ee_v(jpi,jpj) , ee_f(jpi,jpj) , &52 51 & r2dt (jpk) , STAT=dom_vvl_alloc ) 53 52 ! … … 62 61 !! *** ROUTINE dom_vvl *** 63 62 !! 64 !! ** Purpose : compute coefficients muX at T-U-V-F points to spread 65 !! ssh over the whole water column (scale factors) 63 !! ** Purpose : compute mu coefficients at t-, u-, v- and f-points to 64 !! spread ssh over the whole water column (scale factors) 65 !! set the before and now ssh at u- and v-points 66 !! (also f-point in now case) 66 67 !!---------------------------------------------------------------------- 67 68 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 68 USE wrk_nemo, ONLY: z s_t => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3! 2D workspace69 USE wrk_nemo, ONLY: zee_t => wrk_2d_1, zee_u => wrk_2d_2, zee_v => wrk_2d_3, zee_f => wrk_2d_4 ! 2D workspace 69 70 ! 70 71 INTEGER :: ji, jj, jk ! dummy loop indices 71 REAL(wp) :: zcoefu , zcoefv , zcoeff! local scalars72 REAL(wp) :: zv _t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1 ! - -73 !!---------------------------------------------------------------------- 74 75 IF( wrk_in_use(2, 1,2,3 ) ) THEN72 REAL(wp) :: zcoefu, zcoefv , zcoeff ! local scalars 73 REAL(wp) :: zvt , zvt_ip1, zvt_jp1, zvt_ip1jp1 ! - - 74 !!---------------------------------------------------------------------- 75 76 IF( wrk_in_use(2, 1,2,3,4) ) THEN 76 77 CALL ctl_stop('dom_vvl: requested workspace arrays unavailable') ; RETURN 77 78 ENDIF … … 97 98 98 99 ! !== mu computation ==! 99 ee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level100 ee_u(:,:) = fse3u_0(:,:,1)101 ee_v(:,:) = fse3v_0(:,:,1)102 ee_f(:,:) = fse3f_0(:,:,1)100 zee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 101 zee_u(:,:) = fse3u_0(:,:,1) 102 zee_v(:,:) = fse3v_0(:,:,1) 103 zee_f(:,:) = fse3f_0(:,:,1) 103 104 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 104 ee_t(:,:) =ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)105 ee_u(:,:) =ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)106 ee_v(:,:) =ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)105 zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 106 zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 107 zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 107 108 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 108 ee_f(:,jj) =ee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)109 zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 109 110 END DO 110 111 END DO 111 112 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 112 ee_t(:,:) = 1. /ee_t(:,:) * tmask(:,:,1)113 ee_u(:,:) = 1. /ee_u(:,:) * umask(:,:,1)114 ee_v(:,:) = 1. /ee_v(:,:) * vmask(:,:,1)113 zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1) 114 zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1) 115 zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1) 115 116 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 116 ee_f(:,jj) = 1. /ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)117 END DO 118 CALL lbc_lnk( ee_f, 'F', 1. )! lateral boundary condition on ee_f117 zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 118 END DO 119 CALL lbc_lnk( zee_f, 'F', 1. ) ! lateral boundary condition on ee_f 119 120 ! 120 121 DO jk = 1, jpk ! mu coefficients 121 mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk) ! T-point at T levels122 muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk) ! U-point at T levels123 muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk) ! V-point at T levels122 mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk) ! T-point at T levels 123 muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk) ! U-point at T levels 124 muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk) ! V-point at T levels 124 125 END DO 125 126 DO jk = 1, jpk ! F-point : fmask=shlat at coasts, use the product of umask 126 127 DO jj = 1, jpjm1 127 muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels128 END DO 129 muf(:,jpj,jk) = 0. e0128 muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 129 END DO 130 muf(:,jpj,jk) = 0._wp 130 131 END DO 131 132 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition … … 139 140 END DO 140 141 141 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations142 ! for ssh and scale factors143 zs_t (:,:) = e1t(:,:) * e2t(:,:)144 zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )145 zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )146 147 142 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 148 143 DO ji = 1, jpim1 ! NO vector opt. 149 zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)150 zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)151 zcoeff = 0. 5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj))152 ! before fields153 zv _t_ij = zs_t(ji ,jj ) * sshb(ji ,jj )154 zv _t_ip1j = zs_t(ji+1,jj ) * sshb(ji+1,jj )155 zv _t_ijp1 = zs_t(ji ,jj+1) * sshb(ji ,jj+1)156 sshu_b(ji,jj) = zcoefu * ( zv _t_ij + zv_t_ip1j)157 sshv_b(ji,jj) = zcoefv * ( zv _t_ij + zv_t_ijp1 )158 ! now fields159 zv _t_ij = zs_t(ji ,jj ) * sshn(ji ,jj )160 zv _t_ip1j = zs_t(ji+1,jj ) * sshn(ji+1,jj )161 zv _t_ijp1 = zs_t(ji ,jj+1) * sshn(ji ,jj+1)162 zv _t_ip1jp1 = zs_t(ji ,jj+1) * sshn(ji,jj+1)163 sshu_n(ji,jj) = zcoefu * ( zv _t_ij + zv_t_ip1j)164 sshv_n(ji,jj) = zcoefv * ( zv _t_ij + zv_t_ijp1 )165 sshf_n(ji,jj) = zcoeff * ( zv _t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )144 zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1) 145 zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1) 146 zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1) 147 ! 148 zvt = e1e2t(ji ,jj ) * sshb(ji ,jj ) ! before fields 149 zvt_ip1 = e1e2t(ji+1,jj ) * sshb(ji+1,jj ) 150 zvt_jp1 = e1e2t(ji ,jj+1) * sshb(ji ,jj+1) 151 sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 152 sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 153 ! 154 zvt = e1e2t(ji ,jj ) * sshn(ji ,jj ) ! now fields 155 zvt_ip1 = e1e2t(ji+1,jj ) * sshn(ji+1,jj ) 156 zvt_jp1 = e1e2t(ji ,jj+1) * sshn(ji ,jj+1) 157 zvt_ip1jp1 = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) 158 sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 159 sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 160 sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 ) 166 161 END DO 167 162 END DO … … 169 164 CALL lbc_lnk( sshv_n, 'V', 1. ) ; CALL lbc_lnk( sshv_b, 'V', 1. ) 170 165 CALL lbc_lnk( sshf_n, 'F', 1. ) 171 172 ! initialise before scale factors at (u/v)-points 173 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 174 DO jk = 1, jpkm1 175 DO jj = 1, jpjm1 176 DO ji = 1, jpim1 177 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 178 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 179 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 180 fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 181 fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 166 ! 167 IF( wrk_not_released(2, 1,2,3,4) ) CALL ctl_stop('dom_vvl: failed to release workspace arrays') 168 ! 169 END SUBROUTINE dom_vvl 170 171 172 SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b ) 173 !!---------------------------------------------------------------------- 174 !! *** ROUTINE dom_vvl_2 *** 175 !! 176 !! ** Purpose : compute the vertical scale factors at u- and v-points 177 !! in variable volume case. 178 !! 179 !! ** Method : In variable volume case (non linear sea surface) the 180 !! the vertical scale factor at velocity points is computed 181 !! as the average of the cell surface weighted e3t. 182 !! It uses the sea surface heigth so it have to be initialized 183 !! after ssh is read/set 184 !!---------------------------------------------------------------------- 185 INTEGER , INTENT(in ) :: kt ! ocean time-step index 186 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3u_b, pe3v_b ! before vertical scale factor at u- & v-pts 187 ! 188 INTEGER :: ji, jj, jk ! dummy loop indices 189 INTEGER :: iku, ikv ! local integers 190 REAL(wp) :: zvt ! local scalars 191 !!---------------------------------------------------------------------- 192 193 IF( lwp .AND. kt == nit000 ) THEN 194 WRITE(numout,*) 195 WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 196 WRITE(numout,*) '~~~~~~~~~ ' 197 pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 198 pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 199 ENDIF 200 201 DO jk = 1, jpkm1 ! set the before scale factors at u- & v-points 202 DO jj = 2, jpjm1 203 DO ji = fs_2, fs_jpim1 204 zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 205 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) ) 206 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) ) 182 207 END DO 183 208 END DO 184 209 END DO 185 CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 186 CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 187 ! Add initial scale factor to scale factor anomaly 188 fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 189 fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 190 ! 191 IF( wrk_not_released(2, 1,2,3) ) CALL ctl_stop('dom_vvl: failed to release workspace arrays') 192 ! 193 END SUBROUTINE dom_vvl 194 210 IF( ln_zps ) THEN ! minimum of the e3t at partial cell level 211 DO jj = 2, jpjm1 212 DO ji = fs_2, fs_jpim1 213 iku = mbku(ji,jj) 214 ikv = mbkv(ji,jj) 215 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 216 pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji ,jj+1,ikv) ) 217 END DO 218 END DO 219 ENDIF 220 pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:) ! anomaly to avoid zero along closed boundary/extra halos 221 pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 222 CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 223 CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 224 pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:) ! recover the full scale factor 225 pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 226 ! 227 END SUBROUTINE dom_vvl_2 228 195 229 #else 196 230 !!---------------------------------------------------------------------- … … 200 234 SUBROUTINE dom_vvl 201 235 END SUBROUTINE dom_vvl 236 SUBROUTINE dom_vvl_2 237 END SUBROUTINE dom_vvl_2 202 238 #endif 203 239
Note: See TracChangeset
for help on using the changeset viewer.
