Changeset 5989 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
- Timestamp:
- 2015-12-03T09:10:32+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5260 r5989 20 20 USE sbc_oce ! surface boundary condition: ocean 21 21 USE domvvl ! Variable volume 22 USE divcur ! hor. divergence and curl (div & cur routines) 23 USE iom ! I/O library 24 USE restart ! only for lrst_oce 25 USE in_out_manager ! I/O manager 26 USE prtctl ! Print control 27 USE phycst 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) 29 USE lib_mpp ! MPP library 22 USE divhor ! horizontal divergence 23 USE phycst ! physical constants 30 24 USE bdy_oce 31 25 USE bdy_par 32 26 USE bdydyn2d ! bdy_ssh routine 33 USE iom34 27 #if defined key_agrif 35 USE agrif_opa_update36 28 USE agrif_opa_interp 37 29 #endif … … 39 31 USE asminc ! Assimilation increment 40 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 41 38 USE wrk_nemo ! Memory Allocation 42 39 USE timing ! Timing … … 69 66 !! by the time step. 70 67 !! 71 !! ** action : ssha :after sea surface height68 !! ** action : ssha, after sea surface height 72 69 !! 73 70 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 74 71 !!---------------------------------------------------------------------- 75 ! 76 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 77 INTEGER, INTENT(in) :: kt ! time step 72 INTEGER, INTENT(in) :: kt ! time step 78 73 ! 79 INTEGER :: jk ! dummy loop indice 80 REAL(wp) :: z2dt, z1_rau0 ! local scalars 81 !!---------------------------------------------------------------------- 82 ! 83 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 84 ! 85 CALL wrk_alloc( jpi, jpj, zhdiv ) 74 INTEGER :: jk ! dummy loop indice 75 REAL(wp) :: z2dt, zcoef ! local scalars 76 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv ! 2D workspace 77 !!---------------------------------------------------------------------- 78 ! 79 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 80 ! 81 CALL wrk_alloc( jpi,jpj, zhdiv ) 86 82 ! 87 83 IF( kt == nit000 ) THEN 88 !89 84 IF(lwp) WRITE(numout,*) 90 85 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 91 86 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 92 ! 93 ENDIF 94 ! 95 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 87 ENDIF 88 ! 89 CALL div_hor( kt ) ! Horizontal divergence 96 90 ! 97 91 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) … … 109 103 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 110 104 ! 111 z 1_rau0= 0.5_wp * r1_rau0112 ssha(:,:) = ( sshb(:,:) - z2dt * ( z 1_rau0* ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)113 114 #if ! defined key_dynspg_ts 115 ! These lines are not necessary with time splitting since116 ! boundary condition on sea level is set during ts loop117 # if defined key_agrif118 CALL agrif_ssh( kt )119 # endif120 # if defined key_bdy121 IF (lk_bdy) THEN122 CALL lbc_lnk( ssha, 'T', 1. )! Not sure that's necessary123 CALL bdy_ssh( ssha )! Duplicate sea level across open boundaries124 ENDIF125 # endif126 #endif 105 zcoef = 0.5_wp * r1_rau0 106 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 107 108 IF ( .NOT.ln_dynspg_ts ) THEN 109 ! These lines are not necessary with time splitting since 110 ! boundary condition on sea level is set during ts loop 111 # if defined key_agrif 112 CALL agrif_ssh( kt ) 113 # endif 114 # if defined key_bdy 115 IF( lk_bdy ) THEN 116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 118 ENDIF 119 # endif 120 ENDIF 127 121 128 122 #if defined key_asminc 129 ! ! Include the IAU weighted SSH increment 130 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 123 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN ! Include the IAU weighted SSH increment 131 124 CALL ssh_asm_inc( kt ) 132 125 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 133 126 ENDIF 134 127 #endif 135 136 128 ! !------------------------------! 137 129 ! ! outputs ! 138 130 ! !------------------------------! 139 CALL iom_put( "ssh" , sshn ) ! sea surface height140 if( iom_use('ssh2') ) CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height141 131 ! 142 132 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) … … 164 154 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 165 155 !!---------------------------------------------------------------------- 166 ! 167 INTEGER, INTENT(in) :: kt ! time step 156 INTEGER, INTENT(in) :: kt ! time step 157 ! 158 INTEGER :: ji, jj, jk ! dummy loop indices 159 REAL(wp) :: z1_2dt ! local scalars 168 160 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 169 161 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 170 ! 171 INTEGER :: ji, jj, jk ! dummy loop indices 172 REAL(wp) :: z1_2dt ! local scalars 173 !!---------------------------------------------------------------------- 174 175 IF( nn_timing == 1 ) CALL timing_start('wzv') 162 !!---------------------------------------------------------------------- 163 ! 164 IF( nn_timing == 1 ) CALL timing_start('wzv') 176 165 ! 177 166 IF( kt == nit000 ) THEN 178 !179 167 IF(lwp) WRITE(numout,*) 180 168 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' … … 182 170 ! 183 171 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 184 !185 172 ENDIF 186 173 ! !------------------------------! … … 198 185 DO jj = 2, jpjm1 199 186 DO ji = fs_2, fs_jpim1 ! vector opt. 200 zhdiv(ji,jj,jk) = r1_e1 2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )187 zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 201 188 END DO 202 189 END DO … … 221 208 222 209 #if defined key_bdy 223 IF (lk_bdy) THEN210 IF( lk_bdy ) THEN 224 211 DO jk = 1, jpkm1 225 212 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) … … 228 215 #endif 229 216 ! 230 ! !------------------------------!231 ! ! outputs !232 ! !------------------------------!233 CALL iom_put( "woce", wn ) ! vertical velocity234 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value235 CALL wrk_alloc( jpi, jpj, z2d )236 CALL wrk_alloc( jpi, jpj, jpk, z3d )237 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.238 z2d(:,:) = rau0 * e12t(:,:)239 DO jk = 1, jpk240 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)241 END DO242 CALL iom_put( "w_masstr" , z3d )243 IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )244 CALL wrk_dealloc( jpi, jpj, z2d )245 CALL wrk_dealloc( jpi, jpj, jpk, z3d )246 ENDIF247 !248 217 IF( nn_timing == 1 ) CALL timing_stop('wzv') 249 250 218 ! 251 219 END SUBROUTINE wzv 220 252 221 253 222 SUBROUTINE ssh_swp( kt ) … … 281 250 ENDIF 282 251 283 # if defined key_dynspg_ts 284 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 285 # else 286 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 287 #endif 252 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 253 !** Euler time-stepping: no filter 288 254 sshb(:,:) = sshn(:,:) ! before <-- now 289 255 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 256 ! 290 257 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 291 258 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 292 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * ssmask(:,:) 259 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) & 260 & - rnf_b(:,:) + rnf(:,:) & 261 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 293 262 sshn(:,:) = ssha(:,:) ! now <-- after 294 263 ENDIF 295 !296 ! Update velocity at AGRIF zoom boundaries297 #if defined key_agrif298 IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )299 #endif300 264 ! 301 265 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
Note: See TracChangeset
for help on using the changeset viewer.
