- Timestamp:
- 2013-11-20T17:28:04+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90
r4230 r4292 95 95 ! Update data, open boundaries, surface boundary condition (including sea-ice) 96 96 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 97 IF( lk_tide ) CALL sbc_tide( kstp ) 98 IF( lk_obc ) CALL obc_dta ( kstp ) ! update dynamic and tracer data at open boundaries 99 IF( lk_obc ) CALL obc_rad ( kstp ) ! compute phase velocities at open boundaries 100 IF( lk_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 101 97 102 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 98 99 IF( lk_tide.AND.(kstp /= nit000 )) CALL tide_init ( kstp ) 100 IF( lk_tide ) CALL sbc_tide( kstp ) 101 IF( lk_obc ) CALL obc_dta( kstp ) ! update dynamic and tracer data at open boundaries 102 IF( lk_obc ) CALL obc_rad( kstp ) ! compute phase velocities at open boundaries 103 IF( lk_bdy ) CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries 104 105 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 106 ! Ocean dynamics : ssh, wn, hdiv, rot ! 107 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 108 CALL ssh_wzv( kstp ) ! after ssh & vertical velocity 103 ! clem: moved here for bdy ice purpose 104 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 105 ! Ocean dynamics : hdiv, rot, ssh, e3, wn 106 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 107 CALL zdf_bfr( kstp ) ! bottom friction (if quadratic) 108 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_cur) 109 IF( lk_dynspg_ts ) THEN 110 CALL wzv ( kstp ) ! now cross-level velocity 111 ! In case the time splitting case, update almost all momentum trends here: 112 ! Note that the computation of vertical velocity above, hence "after" sea level 113 ! is necessary to compute momentum advection for the rhs of barotropic loop: 114 CALL eos ( tsn, rhd, rhop ) ! now in situ density for hpg computation 115 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 116 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 117 118 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 119 va(:,:,:) = 0.e0 120 IF( ln_asmiau .AND. & 121 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 122 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! subtract Neptune velocities (simplified) 123 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 124 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 125 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 126 CALL dyn_ldf ( kstp ) ! lateral mixing 127 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! add Neptune velocities (simplified) 128 #if defined key_agrif 129 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 130 #endif 131 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 132 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 133 134 hdivb(:,:,:) = hdivn(:,:,:) ! Store now divergence and rot temporarly, revert to these below 135 rotb(:,:,:) = rotn(:,:,:) 136 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated) 137 va_sv(:,:,:) = va(:,:,:) 138 139 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 140 ENDIF 141 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 142 CALL wzv ( kstp ) ! now cross-level velocity (original) 109 143 110 144 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 115 149 ! 116 150 ! VERTICAL PHYSICS 117 CALL zdf_bfr( kstp ) ! bottom friction118 119 151 ! ! Vertical eddy viscosity and diffusivity coefficients 120 152 IF( lk_zdfric ) CALL zdf_ric( kstp ) ! Richardson number dependent Kz … … 122 154 IF( lk_zdfgls ) CALL zdf_gls( kstp ) ! GLS closure scheme for Kz 123 155 IF( lk_zdfkpp ) CALL zdf_kpp( kstp ) ! KPP closure scheme for Kz 124 IF( lk_zdfcst ) THEN! Constant Kz (reset avt, avm[uv] to the background value)156 IF( lk_zdfcst ) THEN ! Constant Kz (reset avt, avm[uv] to the background value) 125 157 avt (:,:,:) = rn_avt0 * tmask(:,:,:) 126 158 avmu(:,:,:) = rn_avm0 * umask(:,:,:) … … 146 178 ! 147 179 IF( lk_ldfslp ) THEN ! slope of lateral mixing 148 CALL eos( tsb, rhd )! before in situ density180 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 149 181 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 150 182 & rhd, gru , grv ) ! of t, s, rd at the last ocean level … … 193 225 tsa(:,:,:,:) = 0.e0 ! set tracer trends to zero 194 226 195 !write(numout,*) "MAV kt",kstp196 !write(numout,'(a5,3(1x,f21.18))') "INIn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)197 !write(numout,'(a5,3(1x,f21.18))') "INIa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)198 227 IF( ln_asmiau .AND. & 199 228 & ln_trainc ) CALL tra_asm_inc( kstp ) ! apply tracer assimilation increment … … 205 234 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 206 235 CALL tra_adv ( kstp ) ! horizontal & vertical advection 207 !write(numout,'(a5,3(1x,f21.18))') "ADVn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)208 !write(numout,'(a5,3(1x,f21.18))') "ADVa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)209 236 IF( lk_zdfkpp ) CALL tra_kpp ( kstp ) ! KPP non-local tracer fluxes 210 237 CALL tra_ldf ( kstp ) ! lateral mixing 211 !write(numout,'(a5,3(1x,f21.18))') "LDFn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11)212 !write(numout,'(a5,3(1x,f21.18))') "LDFa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)213 238 #if defined key_agrif 214 239 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 215 240 #endif 216 241 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 217 !do jk=1,jpk218 !write(numout,'(a5,3(1x,f21.18))') "ZDFn:",tsn(5,10,jk,jp_tem),tsn(5,10,jk,jp_sal),tmask(5,10,jk)219 !write(numout,'(a5,3(1x,f21.18))') "ZDFa:",tsa(5,10,jk,jp_tem),tsa(5,10,jk,jp_sal),ssha(5,10)220 !end do221 242 222 243 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 223 244 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 224 245 CALL tra_nxt( kstp ) ! tracer fields at next time step 225 CALL eos ( tsa, rhd, rhop )! Time-filtered in situ density for hpg computation246 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 226 247 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv, & ! zps: time filtered hor. derivative 227 248 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 228 249 229 250 ELSE ! centered hpg (eos then time stepping) 230 CALL eos ( tsn, rhd, rhop ) ! now in situ density for hpg computation 231 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 251 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 252 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 253 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 232 254 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 233 !write(numout,'(a5,3(1x,f21.18))') "ZPSn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(24,11) 234 !write(numout,'(a5,3(1x,f21.18))') "ZPSa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11) 255 ENDIF 235 256 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 236 257 CALL tra_nxt( kstp ) ! tracer fields at next time step 237 !write(numout,'(a5,3(1x,f21.18))') "NXTn:",tsn(24,11,1,jp_tem),tsn(24,11,1,jp_sal),sshn(25,11)238 !write(numout,'(a5,3(1x,f21.18))') "NXTa:",tsa(24,11,1,jp_tem),tsa(24,11,1,jp_sal),ssha(24,11)239 258 ENDIF 240 259 … … 242 261 ! Dynamics (tsa used as workspace) 243 262 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 263 IF( lk_dynspg_ts ) THEN 264 ! revert to previously computed momentum tendencies 265 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 266 ua(:,:,:) = ua_sv(:,:,:) 267 va(:,:,:) = va_sv(:,:,:) 268 ! Revert now divergence and rotational to previously computed ones 269 !(needed because of the time swap in div_cur, at the beginning of each time step) 270 hdivn(:,:,:) = hdivb(:,:,:) 271 rotn(:,:,:) = rotb(:,:,:) 272 273 CALL dyn_bfr( kstp ) ! bottom friction 274 CALL dyn_zdf( kstp ) ! vertical diffusion 275 ELSE 244 276 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 245 277 va(:,:,:) = 0.e0 246 278 247 IF( ln_asmiau .AND. &248 & ln_dyninc )CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment249 IF( ln_bkgwri )CALL asm_bkg_wri( kstp ) ! output background fields250 IF( ln_neptsimp )CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified)251 IF( lk_bdy )CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends279 IF( ln_asmiau .AND. & 280 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 281 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 282 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 283 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 252 284 CALL dyn_adv( kstp ) ! advection (vector or flux form) 253 285 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 254 286 CALL dyn_ldf( kstp ) ! lateral mixing 255 IF( ln_neptsimp )CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified)256 #if defined key_agrif 257 IF(.NOT. Agrif_Root())CALL Agrif_Sponge_dyn ! momemtum sponge287 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 288 #if defined key_agrif 289 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 258 290 #endif 259 291 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure … … 261 293 CALL dyn_zdf( kstp ) ! vertical diffusion 262 294 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 295 ENDIF 263 296 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 264 297 265 CALL ssh_nxt( kstp ) ! sea surface height at next time step 298 CALL ssh_swp( kstp ) ! swap of sea surface height 299 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 266 300 267 301 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics
Note: See TracChangeset
for help on using the changeset viewer.
