Changeset 5930 for trunk/NEMOGCM/NEMO/OPA_SRC
- Timestamp:
- 2015-11-26T17:07:10+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 4 deleted
- 38 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
r5836 r5930 34 34 USE zdfmxl ! Mixed layer depth 35 35 USE dom_oce, ONLY : ndastp 36 USE sol_oce, ONLY : gcx ! Solver variables defined in memory37 36 USE in_out_manager ! I/O manager 38 37 USE iom ! I/O module … … 114 113 CALL iom_rstput( kt, nitbkg_r, inum, 'en' , en ) 115 114 #endif 116 CALL iom_rstput( kt, nitbkg_r, inum, 'gcx' , gcx )117 115 ! 118 116 CALL iom_close( inum ) -
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r5919 r5930 29 29 USE iom ! IOM library 30 30 USE in_out_manager ! I/O logical units 31 USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag32 31 #if defined key_lim2 33 32 USE ice_2 … … 388 387 END DO ! ib_bdy 389 388 390 ! bg jchanut tschanges391 389 #if defined key_tide 392 ! Add tides if not split-explicit free surface else this is done in ts loop 393 IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 394 #endif 395 ! end jchanut tschanges 390 IF (ln_dynspg_ts) THEN ! Fill temporary arrays with slow-varying bdy data 391 DO ib_bdy = 1, nb_bdy ! Tidal component added in ts loop 392 IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 393 nblen => idx_bdy(ib_bdy)%nblen 394 nblenrim => idx_bdy(ib_bdy)%nblenrim 395 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 396 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen1(1)) = dta_bdy(ib_bdy)%ssh(1:ilen1(1)) 397 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen1(2)) = dta_bdy(ib_bdy)%u2d(1:ilen1(2)) 398 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen1(3)) = dta_bdy(ib_bdy)%v2d(1:ilen1(3)) 399 ENDIF 400 END DO 401 ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 402 ! 403 CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 404 ENDIF 405 #endif 396 406 397 407 IF ( ln_apr_obc ) THEN … … 423 433 !! 424 434 !!---------------------------------------------------------------------- 425 USE dynspg_oce, ONLY: lk_dynspg_ts426 435 !! 427 436 INTEGER :: ib_bdy, jfld, jstart, jend, ierror ! local indices … … 670 679 ! sea ice 671 680 IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 672 ! Test for types of ice input (lim2 or lim3) 681 ! Test for types of ice input (lim2 or lim3) 673 682 ! Build file name to find dimensions 674 683 clname=TRIM(bn_a_i%clname) -
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
r4689 r5930 24 24 USE oce ! ocean dynamics and tracers 25 25 USE dom_oce ! ocean space and time domain 26 USE dynspg_oce27 26 USE bdy_oce ! ocean open boundary conditions 28 27 USE bdydyn2d ! open boundary conditions for barotropic solution … … 35 34 PRIVATE 36 35 37 PUBLIC bdy_dyn ! routine called in dynspg_flt (if lk_dynspg_flt) or 38 ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp) 36 PUBLIC bdy_dyn ! routine called in dyn_nxt 39 37 40 38 # include "domzgr_substitute.h90" -
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
r5215 r5930 23 23 USE bdy_oce ! ocean open boundary conditions 24 24 USE bdylib ! BDY library routines 25 USE dynspg_oce ! for barotropic variables26 25 USE phycst ! physical constants 27 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) -
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r5913 r5930 33 33 ! USE tide_mod ! Useless ?? 34 34 USE fldread 35 USE dynspg_oce, ONLY: lk_dynspg_ts36 35 37 36 IMPLICIT NONE … … 54 53 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 55 54 !$AGRIF_END_DO_NOT_TREAT 56 TYPE(OBC_DATA) , P RIVATE, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component)55 TYPE(OBC_DATA) , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component) 57 56 58 57 !!---------------------------------------------------------------------- … … 270 269 ENDIF 271 270 ! 272 IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 273 ! time splitting integration 274 ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 275 ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 276 ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 277 dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 278 dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 279 dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 280 ENDIF 271 ! Allocate slow varying data in the case of time splitting: 272 ! Do it anyway because at this stage knowledge of free surface scheme is unknown 273 ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 274 ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 275 ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 276 dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 277 dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 278 dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 281 279 ! 282 280 ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 … … 397 395 !! 398 396 LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step 399 INTEGER, 397 INTEGER, DIMENSION(jpbgrd) :: ilen0 400 398 INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim ! short cuts 401 399 INTEGER :: itide, ib_bdy, ib, igrd ! loop indices … … 456 454 zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 457 455 ! 458 ! If time splitting, save data at first barotropic iteration 459 IF ( PRESENT(kit) ) THEN 460 IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 461 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 462 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 463 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 464 465 ELSE ! Initialize arrays from slow varying open boundary data: 466 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 467 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 468 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 469 ENDIF 456 ! If time splitting, initialize arrays from slow varying open boundary data: 457 IF ( PRESENT(kit) ) THEN 458 IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 459 IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 460 IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 470 461 ENDIF 471 462 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
r5836 r5930 10 10 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 11 11 !!---------------------------------------------------------------------- 12 #if defined key_bdy && defined key_dynspg_flt12 #if defined key_bdy 13 13 !!---------------------------------------------------------------------- 14 !! 'key_bdy' AND unstructured open boundary conditions 15 !! 'key_dynspg_flt' filtered free surface 14 !! 'key_bdy' unstructured open boundary conditions 16 15 !!---------------------------------------------------------------------- 17 16 USE oce ! ocean dynamics and tracers … … 30 29 PRIVATE 31 30 32 PUBLIC bdy_vol ! routine called by dynspg_flt.h9031 PUBLIC bdy_vol 33 32 34 33 !! * Substitutions … … 45 44 !! *** ROUTINE bdyvol *** 46 45 !! 47 !! ** Purpose : This routine is called in dynspg_flt to control48 !! the volume of the system.A correction velocity is calculated46 !! ** Purpose : This routine controls the volume of the system. 47 !! A correction velocity is calculated 49 48 !! to correct the total transport through the unstructured OBC. 50 49 !! The total depth used is constant (H0) to be consistent with the -
trunk/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90
r5836 r5930 18 18 #endif 19 19 USE dyncor_c1d ! Coriolis term (c1d case) (dyn_cor_1d ) 20 USE dynnxt _c1d! time-stepping (dyn_nxt routine)20 USE dynnxt ! time-stepping (dyn_nxt routine) 21 21 USE dyndmp ! U & V momentum damping (dyn_dmp routine) 22 22 USE restart ! restart … … 139 139 CALL dyn_cor_c1d( kstp ) ! vorticity term including Coriolis 140 140 CALL dyn_zdf ( kstp ) ! vertical diffusion 141 CALL dyn_nxt _c1d( kstp ) ! lateral velocity at next time step141 CALL dyn_nxt ( kstp ) ! lateral velocity at next time step 142 142 143 143 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r5586 r5930 14 14 USE dom_oce ! ocean space and time domain 15 15 USE phycst 16 USE dynspg_oce17 USE dynspg_ts18 16 USE daymod 19 17 USE tide_mod -
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r5836 r5930 30 30 USE zdf_oce ! ocean vertical physics 31 31 USE ldftra ! lateral physics: eddy diffusivity coef. 32 USE sol_oce ! solver variables33 32 USE sbc_oce ! Surface boundary condition: ocean fields 34 33 USE sbc_ice ! Surface boundary condition: ice fields … … 46 45 USE iom 47 46 USE ioipsl 48 USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities49 47 50 48 #if defined key_lim2 … … 206 204 CALL iom_put( "sbu", z2d ) ! bottom i-current 207 205 ENDIF 208 #if defined key_dynspg_ts 209 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 210 #else 211 CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current 212 #endif 206 207 IF ( ln_dynspg_ts ) THEN 208 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 209 ELSE 210 CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current 211 ENDIF 213 212 214 213 CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current … … 223 222 CALL iom_put( "sbv", z2d ) ! bottom j-current 224 223 ENDIF 225 #if defined key_dynspg_ts 226 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current 227 #else 228 CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current 229 #endif 224 225 IF ( ln_dynspg_ts ) THEN 226 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current 227 ELSE 228 CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current 229 ENDIF 230 230 231 231 CALL iom_put( "woce", wn ) ! vertical velocity -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5836 r5930 45 45 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers 46 46 47 !! Free surface parameters 48 !! ======================= 49 LOGICAL, PUBLIC :: ln_dynspg_exp !: Explicit free surface flag 50 LOGICAL, PUBLIC :: ln_dynspg_ts !: Split-Explicit free surface flag 51 47 52 !! Time splitting parameters 48 53 !! ========================= 49 54 LOGICAL, PUBLIC :: ln_bt_fw !: Forward integration of barotropic sub-stepping 50 55 LOGICAL, PUBLIC :: ln_bt_av !: Time averaging of barotropic variables 51 LOGICAL, PUBLIC :: ln_bt_ nn_auto!: Set number of barotropic iterations automatically56 LOGICAL, PUBLIC :: ln_bt_auto !: Set number of barotropic iterations automatically 52 57 INTEGER, PUBLIC :: nn_bt_flt !: Filter choice 53 58 INTEGER, PUBLIC :: nn_baro !: Number of barotropic iterations during one baroclinic step (rdt) 54 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_ nn_auto=T)59 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_auto=T) 55 60 56 61 !! Horizontal grid parameters for domhgr -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5836 r5930 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE lib_mpp 30 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient31 30 USE wrk_nemo ! Memory allocation 32 31 USE timing ! Timing -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5836 r5930 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 36 USE domvvl ! varying vertical mesh 37 USE dynspg_oce ! pressure gradient schemes38 USE dynspg_flt ! filtered free surface39 USE sol_oce ! ocean solver variables40 37 ! 41 38 USE in_out_manager ! I/O manager … … 133 130 ENDIF 134 131 ! 135 IF( lk_agrif ) THEN ! read free surface arrays in restart file136 IF( ln_rstart ) THEN137 IF( lk_dynspg_flt ) THEN ! read or initialize the following fields138 ! ! gcx, gcxb for agrif_opa_init139 IF( sol_oce_alloc() > 0 ) CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')140 CALL flt_rst( nit000, 'READ' )141 ENDIF142 ENDIF ! explicit case not coded yet with AGRIF143 ENDIF144 !145 132 ! 146 133 ! Initialize "now" and "before" barotropic velocities: … … 445 432 !! p=integral [ rau*g dz ] 446 433 !!---------------------------------------------------------------------- 447 USE dynspg ! surface pressure gradient (dyn_spg routine)448 434 USE divhor ! hor. divergence (div_hor routine) 449 435 USE lbclnk ! ocean lateral boundary condition (or mpp link) 450 436 ! 451 437 INTEGER :: ji, jj, jk ! dummy loop indices 452 INTEGER :: indic ! ???453 438 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars 454 439 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn … … 517 502 vb(:,:,:) = vn(:,:,:) 518 503 519 ! WARNING !!!!!520 ! after initializing u and v, we need to calculate the initial streamfunction bsf.521 ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).522 ! to do that, we call dyn_spg with a special trick:523 ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the524 ! right value assuming the velocities have been set up in one time step.525 ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)526 ! sets up s false trend to calculate the barotropic streamfunction.527 528 ua(:,:,:) = ub(:,:,:) / rdt529 va(:,:,:) = vb(:,:,:) / rdt530 531 ! calls dyn_spg. we assume euler time step, starting from rest.532 indic = 0533 CALL dyn_spg( nit000, indic ) ! surface pressure gradient534 !535 ! the new velocity is ua*rdt536 !537 CALL lbc_lnk( ua, 'U', -1. )538 CALL lbc_lnk( va, 'V', -1. )539 540 ub(:,:,:) = ua(:,:,:) * rdt541 vb(:,:,:) = va(:,:,:) * rdt542 ua(:,:,:) = 0.e0543 va(:,:,:) = 0.e0544 un(:,:,:) = ub(:,:,:)545 vn(:,:,:) = vb(:,:,:)546 504 ! 547 505 !!gm Check here call to div_hor should not be necessary -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5836 r5930 36 36 USE trd_oce ! trends: ocean variables 37 37 USE trddyn ! trend manager: dynamics 38 !jc USE zpshde ! partial step: hor. derivative (zps_hde routine) 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 58 59 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 60 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag61 61 62 62 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) … … 132 132 !! 133 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf , ln_dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 135 135 !!---------------------------------------------------------------------- 136 136 ! … … 155 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 156 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 157 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp158 157 ENDIF 159 158 ! … … 293 292 ENDIF 294 293 294 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 295 !jc CALL zps_hde ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 295 296 296 297 ! Local constant initialization … … 491 492 ! iniitialised to 0. zhpi zhpi 492 493 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 495 ! Partial steps: top & bottom before horizontal gradient 496 !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 ) 493 499 494 500 !================================================================================== -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5643 r5930 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 20 !! 3.7 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 21 22 !!------------------------------------------------------------------------- 22 23 … … 28 29 USE sbc_oce ! Surface boundary condition: ocean fields 29 30 USE phycst ! physical constants 30 USE dynspg_oce ! type of surface pressure gradient31 31 USE dynadv ! dynamics: vector invariant versus flux form 32 32 USE domvvl ! variable volume … … 68 68 !! *** ROUTINE dyn_nxt *** 69 69 !! 70 !! ** Purpose : Compute the after horizontal velocity. Apply the boundary71 !! condition on the after velocity, achieve dthe time stepping70 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 71 !! condition on the after velocity, achieve the time stepping 72 72 !! by applying the Asselin filter on now fields and swapping 73 73 !! the fields. 74 74 !! 75 !! ** Method : * After velocity is compute using a leap-frog scheme: 76 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 77 !! Note that with flux form advection and variable volume layer 78 !! (lk_vvl=T), the leap-frog is applied on thickness weighted 79 !! velocity. 80 !! Note also that in filtered free surface (lk_dynspg_flt=T), 81 !! the time stepping has already been done in dynspg module 75 !! ** Method : * Ensure after velocities transport matches time splitting 76 !! estimate (ln_dynspg_ts=T) 82 77 !! 83 78 !! * Apply lateral boundary conditions on after velocity … … 101 96 INTEGER :: ji, jj, jk ! dummy loop indices 102 97 INTEGER :: iku, ikv ! local integers 103 #if ! defined key_dynspg_flt 104 REAL(wp) :: z2dt ! temporary scalar 105 #endif 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 107 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 108 100 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 112 104 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 113 105 ! 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 115 IF( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 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 ) 116 109 ! 117 110 IF( kt == nit000 ) THEN … … 121 114 ENDIF 122 115 123 #if defined key_dynspg_flt 124 ! 125 ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine 126 ! ------------- 127 128 ! Update after velocity on domain lateral boundaries (only local domain required) 129 ! -------------------------------------------------- 130 CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries 131 CALL lbc_lnk( va, 'V', -1. ) 132 ! 133 #else 134 135 # if defined key_dynspg_exp 136 ! Next velocity : Leap-frog time stepping 137 ! ------------- 138 z2dt = 2. * rdt ! Euler or leap-frog time step 139 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 140 ! 141 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 116 IF ( ln_dynspg_ts ) THEN 117 ! Ensure below that barotropic velocities match time splitting estimate 118 ! 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) 121 DO jk = 2, jpkm1 122 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 123 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 124 END DO 142 125 DO jk = 1, jpkm1 143 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 144 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 145 END DO 146 ELSE ! applied on thickness weighted velocity 147 DO jk = 1, jpkm1 148 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 149 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 150 & / fse3u_a(:,:,jk) * umask(:,:,jk) 151 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 152 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 153 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 154 END DO 155 ENDIF 156 # endif 157 158 # if defined key_dynspg_ts 159 !!gm IF ( lk_dynspg_ts ) THEN .... 160 ! Ensure below that barotropic velocities match time splitting estimate 161 ! Compute actual transport and replace it with ts estimate at "after" time step 162 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 163 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 164 DO jk = 2, jpkm1 165 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 167 END DO 168 DO jk = 1, jpkm1 169 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 171 END DO 172 173 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 174 ! Remove advective velocity from "now velocities" 175 ! prior to asselin filtering 176 ! In the forward case, this is done below after asselin filtering 177 ! so that asselin contribution is removed at the same time 178 DO jk = 1, jpkm1 179 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 180 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 181 END DO 182 ENDIF 183 !!gm ENDIF 184 # endif 126 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 127 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 128 END DO 129 130 IF ( .NOT.ln_bt_fw ) THEN 131 ! Remove advective velocity from "now velocities" 132 ! prior to asselin filtering 133 ! In the forward case, this is done below after asselin filtering 134 ! so that asselin contribution is removed at the same time 135 DO jk = 1, jpkm1 136 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 137 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 138 END DO 139 ENDIF 140 ENDIF 185 141 186 142 ! Update after velocity on domain lateral boundaries 187 143 ! -------------------------------------------------- 188 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries189 CALL lbc_lnk( va, 'V', -1. )190 !191 # if defined key_bdy192 ! !* BDY open boundaries193 IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )194 IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )195 196 !!$ Do we need a call to bdy_vol here??197 !198 # endif199 !200 144 # if defined key_agrif 201 145 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 202 146 # endif 203 #endif 204 147 ! 148 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries 149 CALL lbc_lnk( va, 'V', -1. ) 150 ! 151 # if defined key_bdy 152 ! !* 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. ) 155 156 !!$ Do we need a call to bdy_vol here?? 157 ! 158 # endif 159 ! 205 160 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 206 161 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step … … 259 214 ! (used as a now filtered scale factor until the swap) 260 215 ! ---------------------------------------------------- 261 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN216 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 262 217 ! No asselin filtering on thicknesses if forward time splitting 263 218 fse3t_b(:,:,:) = fse3t_n(:,:,:) 264 219 ELSE 265 220 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) … … 336 291 ENDIF 337 292 ! 338 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 339 294 ! Revert "before" velocities to time split estimate 340 295 ! Doing it here also means that asselin filter contribution is removed … … 400 355 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 401 356 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 402 ! 403 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 404 IF( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 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 ) 405 361 ! 406 362 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5836 r5930 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 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! dyn_spg : update the dynamics trend with the lateral diffusion 12 !! dyn_spg_ctl : initialization, namelist read, and parameters control 8 !! 3.7 ! 2015-11 (J. Chanut) Suppression of filtered free surface 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_spg : update the dynamics trend with surface pressure gradient 13 !! dyn_spg_init: initialization, namelist read, and parameters control 13 14 !!---------------------------------------------------------------------- 14 15 USE oce ! ocean dynamics and tracers variables … … 18 19 USE sbc_oce ! surface boundary condition: ocean 19 20 USE sbcapr ! surface boundary condition: atmospheric pressure 20 USE dynspg_oce ! surface pressure gradient variables21 21 USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine) 22 22 USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine) 23 USE dynspg_flt ! surface pressure gradient (dyn_spg_flt routine)24 USE dynadv ! dynamics: vector invariant versus flux form25 USE dynhpg, ONLY: ln_dynhpg_imp26 23 USE sbctide 27 24 USE updtide … … 32 29 USE in_out_manager ! I/O manager 33 30 USE lib_mpp ! MPP library 34 USE solver ! solver initialization35 31 USE wrk_nemo ! Memory Allocation 36 32 USE timing ! Timing … … 43 39 PUBLIC dyn_spg_init ! routine called by opa module 44 40 45 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from l k_dynspg_...41 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from ln_dynspg_... 46 42 47 43 !! * Substitutions … … 55 51 CONTAINS 56 52 57 SUBROUTINE dyn_spg( kt , kindic)53 SUBROUTINE dyn_spg( kt ) 58 54 !!---------------------------------------------------------------------- 59 55 !! *** ROUTINE dyn_spg *** … … 66 62 !!gm the after velocity, in the 2 other (ua,va) are still the trends 67 63 !! 68 !! ** Method : T hreeschemes:64 !! ** Method : Two schemes: 69 65 !! - explicit computation : the spg is evaluated at now 70 !! - filtered computation : the Roulet & madec (2000) technique is used71 66 !! - split-explicit computation: a time splitting technique is used 72 67 !! … … 80 75 ! 81 76 INTEGER, INTENT(in ) :: kt ! ocean time-step index 82 INTEGER, INTENT( out) :: kindic ! solver flag83 77 ! 84 78 INTEGER :: ji, jj, jk ! dummy loop indices … … 103 97 104 98 IF( ln_apr_dyn & ! atmos. pressure 105 .OR. ( .NOT.l k_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting)99 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting) 106 100 .OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice 107 101 ! … … 113 107 END DO 114 108 ! 115 IF( ln_apr_dyn .AND. (.NOT. l k_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==!109 IF( ln_apr_dyn .AND. (.NOT. ln_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==! 116 110 zg_2 = grav * 0.5 117 111 DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh … … 126 120 ! 127 121 ! !== tide potential forcing term ==! 128 IF( .NOT.l k_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case122 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case 129 123 ! 130 124 CALL upd_tide( kt ) ! update tide potential … … 171 165 CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit 172 166 CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 173 CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered174 167 ! 175 168 END SELECT 176 169 ! 177 170 IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics 178 SELECT CASE ( nspg ) 179 CASE ( 0, 1 ) 180 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 181 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 182 CASE( 2 ) 183 z2dt = 2. * rdt 184 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 185 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 186 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 187 END SELECT 171 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 188 173 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 189 174 ! … … 203 188 !! *** ROUTINE dyn_spg_init *** 204 189 !! 205 !! ** Purpose : Control the consistency between cppoptions for190 !! ** Purpose : Control the consistency between namelist options for 206 191 !! surface pressure gradient schemes 207 192 !!---------------------------------------------------------------------- 208 INTEGER :: ioptio 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_flt 209 198 !!---------------------------------------------------------------------- 210 199 ! 211 200 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_init') 212 201 ! 213 IF(lwp) THEN ! Control print 202 REWIND( numnam_ref ) ! Namelist namdyn_spg in reference namelist : Free surface 203 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 206 REWIND( numnam_cfg ) ! Namelist namdyn_spg in configuration namelist : Free surface 207 READ ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 208 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp ) 209 IF(lwm) WRITE ( numond, namdyn_spg ) 210 ! 211 IF(lwp) THEN ! Namelist print 214 212 WRITE(numout,*) 215 213 WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme' 216 214 WRITE(numout,*) '~~~~~~~~~~~' 217 WRITE(numout,*) ' Explicit free surface lk_dynspg_exp = ', lk_dynspg_exp 218 WRITE(numout,*) ' Free surface with time splitting lk_dynspg_ts = ', lk_dynspg_ts 219 WRITE(numout,*) ' Filtered free surface cst volume lk_dynspg_flt = ', lk_dynspg_flt 220 ENDIF 221 222 IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 223 ! (do it now, to set nn_baro, used to allocate some arrays later on) 224 ! ! allocate dyn_spg arrays 225 IF( lk_dynspg_ts ) THEN 226 IF( dynspg_oce_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays') 215 WRITE(numout,*) ' Explicit free surface ln_dynspg_exp = ', ln_dynspg_exp 216 WRITE(numout,*) ' Free surface with time splitting ln_dynspg_ts = ', ln_dynspg_ts 217 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 227 222 IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays') 228 223 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) … … 231 226 ! ! Control of surface pressure gradient scheme options 232 227 ioptio = 0 233 IF(lk_dynspg_exp) ioptio = ioptio + 1 234 IF(lk_dynspg_ts ) ioptio = ioptio + 1 235 IF(lk_dynspg_flt) ioptio = ioptio + 1 228 IF(ln_dynspg_exp) ioptio = ioptio + 1 229 IF(ln_dynspg_ts ) ioptio = ioptio + 1 236 230 ! 237 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 238 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 239 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 240 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 241 ! 242 IF( lk_dynspg_exp) nspg = 0 243 IF( lk_dynspg_ts ) nspg = 1 244 IF( lk_dynspg_flt) nspg = 2 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 245 238 ! 246 239 IF(lwp) THEN … … 248 241 IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface' 249 242 IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme' 250 IF( nspg == 2 ) WRITE(numout,*) ' filtered free surface' 251 ENDIF 252 253 #if defined key_dynspg_flt 254 CALL solver_init( nit000 ) ! Elliptic solver initialisation 255 #endif 256 ! ! Control of hydrostatic pressure choice 257 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 243 ENDIF 258 244 ! 259 245 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_init') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r5836 r5930 7 7 !! 3.2 ! 2009-06 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 8 8 !!---------------------------------------------------------------------- 9 #if defined key_dynspg_exp10 9 !!---------------------------------------------------------------------- 11 !! 'key_dynspg_exp'explicit free surface10 !! explicit free surface 12 11 !!---------------------------------------------------------------------- 13 12 !! dyn_spg_exp : update the momentum trend with the surface … … 31 30 PRIVATE 32 31 33 PUBLIC dyn_spg_exp ! routine called by step.F9032 PUBLIC dyn_spg_exp ! routine called by dyn_spg 34 33 35 34 !! * Substitutions … … 101 100 END SUBROUTINE dyn_spg_exp 102 101 103 #else104 !!----------------------------------------------------------------------105 !! Default case : Empty module No standart explicit free surface106 !!----------------------------------------------------------------------107 CONTAINS108 SUBROUTINE dyn_spg_exp( kt ) ! Empty routine109 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt110 END SUBROUTINE dyn_spg_exp111 #endif112 113 102 !!====================================================================== 114 103 END MODULE dynspg_exp -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5913 r5930 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 14 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts15 15 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts'split explicit free surface16 !! split explicit free surface 17 17 !!---------------------------------------------------------------------- 18 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 23 23 USE sbc_oce ! surface boundary condition: ocean 24 24 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables26 25 USE phycst ! physical constants 27 26 USE dynvor ! vorticity term 28 27 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 28 USE bdytides ! open boundary condition data 30 29 USE bdydyn2d ! open boundary conditions on barotropic variables 31 30 USE sbctide ! tides … … 68 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 68 70 ! Arrays below are saved to allow testing of the "no time averaging" option71 ! If this option is not retained, these could be replaced by temporary arrays72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays73 ubb_e, ub_e, &74 vbb_e, vb_e75 76 69 !! * Substitutions 77 70 # include "domzgr_substitute.h90" … … 88 81 !! *** routine dyn_spg_ts_alloc *** 89 82 !!---------------------------------------------------------------------- 90 INTEGER :: ierr( 3)83 INTEGER :: ierr(4) 91 84 !!---------------------------------------------------------------------- 92 85 ierr(:) = 0 93 86 94 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 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) ) 97 91 98 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) … … 101 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 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 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 104 105 105 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 107 107 ! 108 108 END FUNCTION dyn_spg_ts_alloc … … 146 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 150 REAL(wp) :: zu_spg, zv_spg ! - -151 REAL(wp) :: zhura, zhvra 152 REAL(wp) :: za0, za1, za2, za3 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: z un_e, zvn_e, zsshp2_e148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 150 REAL(wp) :: zu_spg, zv_spg ! - - 151 REAL(wp) :: zhura, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 155 155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 156 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv156 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 158 158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 164 164 ! !* Allocate temporary arrays 165 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)167 CALL wrk_alloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)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 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 485 485 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 486 sshn_e(:,:) = sshn (:,:) 487 zun_e(:,:) = un_b (:,:)488 zvn_e(:,:) = vn_b (:,:)487 un_e (:,:) = un_b (:,:) 488 vn_e (:,:) = vn_b (:,:) 489 489 ! 490 490 hu_e (:,:) = hu (:,:) … … 494 494 ELSE ! CENTRED integration: start from BEFORE fields 495 495 sshn_e(:,:) = sshb (:,:) 496 zun_e(:,:) = ub_b (:,:)497 zvn_e(:,:) = vb_b (:,:)496 un_e (:,:) = ub_b (:,:) 497 vn_e (:,:) = vb_b (:,:) 498 498 ! 499 499 hu_e (:,:) = hu_b (:,:) … … 509 509 va_b (:,:) = 0._wp 510 510 ssha (:,:) = 0._wp ! Sum for after averaged sea level 511 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop512 zv_sum(:,:) = 0._wp511 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 512 vn_adv(:,:) = 0._wp 513 513 ! ! ==================== ! 514 514 DO jn = 1, icycle ! sub-time-step loop ! … … 534 534 535 535 ! Extrapolate barotropic velocities at step jit+0.5: 536 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)537 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)536 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 537 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 538 539 539 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 597 597 ! Sum over sub-time-steps to compute advective velocities 598 598 za2 = wgtbtp2(jn) 599 zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)600 zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)599 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 601 601 ! 602 602 ! Set next sea level: … … 733 733 ! 734 734 ! Add bottom stresses: 735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 737 ! 738 738 ! Surface pressure trend: … … 751 751 DO jj = 2, jpjm1 752 752 DO ji = fs_2, fs_jpim1 ! vector opt. 753 ua_e(ji,jj) = ( zun_e(ji,jj) &753 ua_e(ji,jj) = ( un_e(ji,jj) & 754 754 & + rdtbt * ( zwx(ji,jj) & 755 755 & + zu_trd(ji,jj) & … … 757 757 & ) * umask(ji,jj,1) 758 758 759 va_e(ji,jj) = ( zvn_e(ji,jj) &759 va_e(ji,jj) = ( vn_e(ji,jj) & 760 760 & + rdtbt * ( zwy(ji,jj) & 761 761 & + zv_trd(ji,jj) & … … 772 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 773 773 774 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &774 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 778 778 & ) * zhura 779 779 780 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &780 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 781 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 802 802 #if defined key_bdy 803 803 ! open boundaries 804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 805 #endif 806 806 #if defined key_agrif … … 810 810 ! ! ---- 811 811 ubb_e (:,:) = ub_e (:,:) 812 ub_e (:,:) = zun_e(:,:)813 zun_e(:,:) = ua_e (:,:)812 ub_e (:,:) = un_e (:,:) 813 un_e (:,:) = ua_e (:,:) 814 814 ! 815 815 vbb_e (:,:) = vb_e (:,:) 816 vb_e (:,:) = zvn_e(:,:)817 zvn_e(:,:) = va_e (:,:)816 vb_e (:,:) = vn_e (:,:) 817 vn_e (:,:) = va_e (:,:) 818 818 ! 819 819 sshbb_e(:,:) = sshb_e(:,:) … … 840 840 ! ----------------------------------------------------------------------------- 841 841 ! 842 ! At this stage ssha holds a time averaged value 843 ! ! Sea Surface Height at u-,v- and f-points 844 IF( lk_vvl ) THEN ! (required only in key_vvl case) 842 ! Set advection velocity correction: 843 zwx(:,:) = un_adv(:,:) 844 zwy(:,:) = vn_adv(:,:) 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 846 un_adv(:,:) = zwx(:,:)*hur(:,:) 847 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 848 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 851 END IF 852 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 854 ub2_b(:,:) = zwx(:,:) 855 vb2_b(:,:) = zwy(:,:) 856 ENDIF 857 ! 858 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 860 DO jk=1,jpkm1 861 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 862 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 863 END DO 864 ELSE 865 ! At this stage, ssha has been corrected: compute new depths at velocity points 845 866 DO jj = 1, jpjm1 846 867 DO ji = 1, jpim1 ! NO Vector Opt. … … 854 875 END DO 855 876 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 856 ENDIF 857 ! 858 ! Set advection velocity correction: 859 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 860 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 861 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 862 ELSE 863 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 864 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 865 END IF 866 867 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 868 ub2_b(:,:) = zu_sum(:,:) 869 vb2_b(:,:) = zv_sum(:,:) 870 ENDIF 871 ! 872 ! Update barotropic trend: 873 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 874 DO jk=1,jpkm1 875 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 876 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 877 END DO 878 ELSE 877 ! 879 878 DO jk=1,jpkm1 880 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 915 914 ! 916 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 917 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)918 CALL wrk_dealloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )916 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 917 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 919 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 920 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1064 1063 ! 1065 1064 INTEGER :: ji ,jj 1066 INTEGER :: ios ! Local integer output status for namelist read1067 1065 REAL(wp) :: zxr2, zyr2, zcmax 1068 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1069 1067 !! 1070 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1071 & nn_baro, rn_bt_cmax, nn_bt_flt1072 1068 !!---------------------------------------------------------------------- 1073 1069 ! 1074 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1075 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1076 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1077 1078 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1079 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1080 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1081 IF(lwm) WRITE ( numond, namsplit ) 1082 ! 1083 ! ! Max courant number for ext. grav. waves 1070 ! Max courant number for ext. grav. waves 1084 1071 ! 1085 1072 CALL wrk_alloc( jpi, jpj, zcu ) 1086 1073 ! 1087 IF (lk_vvl) THEN 1088 DO jj = 1, jpj 1089 DO ji =1, jpi 1090 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1091 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1092 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1093 END DO 1094 END DO 1095 ELSE 1096 DO jj = 1, jpj 1097 DO ji =1, jpi 1098 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1099 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1100 zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) ) 1101 END DO 1102 END DO 1103 ENDIF 1104 1074 DO jj = 1, jpj 1075 DO ji =1, jpi 1076 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1077 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1078 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1079 END DO 1080 END DO 1081 ! 1105 1082 zcmax = MAXVAL( zcu(:,:) ) 1106 1083 IF( lk_mpp ) CALL mpp_max( zcmax ) 1107 1084 1108 1085 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1109 IF (ln_bt_ nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1086 IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1110 1087 1111 1088 rdtbt = rdt / REAL( nn_baro , wp ) … … 1115 1092 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1116 1093 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1117 IF( ln_bt_ nn_auto ) THEN1118 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1094 IF( ln_bt_auto ) THEN 1095 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1119 1096 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1120 1097 ELSE 1121 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1098 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1122 1099 ENDIF 1123 1100 … … 1164 1141 END SUBROUTINE dyn_spg_ts_init 1165 1142 1166 #else1167 !!---------------------------------------------------------------------------1168 !! Default case : Empty module No split explicit free surface1169 !!---------------------------------------------------------------------------1170 CONTAINS1171 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1172 dyn_spg_ts_alloc = 01173 END FUNCTION dyn_spg_ts_alloc1174 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1175 INTEGER, INTENT(in) :: kt1176 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1177 END SUBROUTINE dyn_spg_ts1178 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1179 INTEGER , INTENT(in) :: kt ! ocean time-step1180 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1181 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1182 END SUBROUTINE ts_rst1183 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1184 INTEGER , INTENT(in) :: kt ! ocean time-step1185 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1186 END SUBROUTINE dyn_spg_ts_init1187 #endif1188 1189 1143 !!====================================================================== 1190 1144 END MODULE dynspg_ts -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5907 r5930 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 USE c1d ! 1D vertical configuration 34 35 ! 35 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 662 663 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 663 664 ! 664 IF( ioptio /= 1) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )665 IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 665 666 ! 666 667 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r3625 r5930 8 8 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 !! 3.7 ! 2015-11 (J. Chanut) output velocities instead of trends 10 11 !!---------------------------------------------------------------------- 11 12 … … 18 19 USE phycst ! physical constants 19 20 USE zdf_oce ! ocean vertical physics 21 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 20 22 USE sbc_oce ! surface boundary condition: ocean 21 23 USE lib_mpp ! MPP library … … 118 120 ! 119 121 END DO ! End of time splitting 122 123 ! Time step momentum here to be compliant with what is done in the implicit case 124 ! 125 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 126 DO jk = 1, jpkm1 127 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 128 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 129 END DO 130 ELSE ! applied on thickness weighted velocity 131 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) 138 END DO 139 ENDIF 120 140 ! 121 141 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5836 r5930 25 25 USE wrk_nemo ! Memory Allocation 26 26 USE timing ! Timing 27 USE dynadv ! dynamics: vector invariant versus flux form 28 USE dynspg_oce, ONLY: lk_dynspg_ts 27 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 29 28 30 29 IMPLICIT NONE … … 87 86 ENDIF 88 87 89 ! 0. Local constant initialization 90 ! -------------------------------- 91 z1_p2dt = 1._wp / p2dt ! inverse of the timestep 92 93 ! 1. Apply semi-implicit bottom friction 88 ! 1. Time step dynamics 89 ! --------------------- 90 91 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 92 DO jk = 1, jpkm1 93 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 94 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 95 END DO 96 ELSE ! applied on thickness weighted velocity 97 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 94 108 ! -------------------------------------- 95 109 ! Only needed for semi-implicit bottom friction setup. The explicit … … 97 111 ! column vector of the tri-diagonal matrix equation 98 112 ! 99 100 113 IF( ln_bfrimp ) THEN 101 114 DO jj = 2, jpjm1 … … 119 132 ENDIF 120 133 121 #if defined key_dynspg_ts 122 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 123 DO jk = 1, jpkm1 124 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 125 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 126 END DO 127 ELSE ! applied on thickness weighted velocity 128 DO jk = 1, jpkm1 129 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 130 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 131 & / fse3u_a(:,:,jk) * umask(:,:,jk) 132 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 133 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 134 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 135 END DO 136 ENDIF 137 138 IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 134 ! With split-explicit free surface, barotropic stress is treated explicitly 135 ! Update velocities at the bottom. 136 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 137 ! not lead to the effective stress seen over the whole barotropic loop. 138 IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN 139 139 ! remove barotropic velocities: 140 140 DO jk = 1, jpkm1 … … 166 166 END IF 167 167 ENDIF 168 #endif 169 170 ! 2. Vertical diffusion on u 168 169 ! 3. Vertical diffusion on u 171 170 ! --------------------------- 172 171 ! Matrix and second member construction … … 219 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 220 219 DO ji = fs_2, fs_jpim1 ! vector opt. 221 #if defined key_dynspg_ts222 220 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 223 221 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 224 222 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 225 #else226 ua(ji,jj,1) = ub(ji,jj,1) &227 & + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &228 & / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) )229 #endif230 223 END DO 231 224 END DO … … 233 226 DO jj = 2, jpjm1 234 227 DO ji = fs_2, fs_jpim1 235 #if defined key_dynspg_ts236 228 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 237 #else238 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)239 #endif240 229 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 241 230 END DO … … 256 245 END DO 257 246 258 #if ! defined key_dynspg_ts 259 ! Normalization to obtain the general momentum trend ua 260 DO jk = 1, jpkm1 261 DO jj = 2, jpjm1 262 DO ji = fs_2, fs_jpim1 ! vector opt. 263 ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 264 END DO 265 END DO 266 END DO 267 #endif 268 269 ! 3. Vertical diffusion on v 247 ! 4. Vertical diffusion on v 270 248 ! --------------------------- 271 249 ! Matrix and second member construction … … 317 295 ! 318 296 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 319 DO ji = fs_2, fs_jpim1 ! vector opt. 320 #if defined key_dynspg_ts 297 DO ji = fs_2, fs_jpim1 ! vector opt. 321 298 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 322 299 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 323 300 & / ( ze3va * rau0 ) 324 #else325 va(ji,jj,1) = vb(ji,jj,1) &326 & + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &327 & / ( fse3v(ji,jj,1) * rau0 ) )328 #endif329 301 END DO 330 302 END DO … … 332 304 DO jj = 2, jpjm1 333 305 DO ji = fs_2, fs_jpim1 ! vector opt. 334 #if defined key_dynspg_ts335 306 zrhs = va(ji,jj,jk) ! zrhs=right hand side 336 #else337 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)338 #endif339 307 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 340 308 END DO … … 355 323 END DO 356 324 357 ! Normalization to obtain the general momentum trend va358 #if ! defined key_dynspg_ts359 DO jk = 1, jpkm1360 DO jj = 2, jpjm1361 DO ji = fs_2, fs_jpim1 ! vector opt.362 va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt363 END DO364 END DO365 END DO366 #endif367 368 325 ! J. Chanut: Lines below are useless ? 369 !! restore bottom layer avmu(v)326 !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 370 327 IF( ln_bfrimp ) THEN 371 328 DO jj = 2, jpjm1 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5836 r5930 106 106 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 107 107 108 #if ! defined key_dynspg_ts 109 ! These lines are not necessary with time splitting since110 ! boundary condition on sea level is set during ts loop108 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 111 # if defined key_agrif 112 CALL agrif_ssh( kt )112 CALL agrif_ssh( kt ) 113 113 # endif 114 114 # if defined key_bdy 115 IF( lk_bdy ) THEN116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries118 ENDIF115 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 119 # endif 120 #endif 120 ENDIF 121 121 122 122 #if defined key_asminc … … 250 250 ENDIF 251 251 252 # if defined key_dynspg_ts 253 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 254 # else 255 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 256 #endif 252 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 253 !** Euler time-stepping: no filter 257 254 sshb(:,:) = sshn(:,:) ! before <-- now 258 255 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
r5836 r5930 12 12 USE dom_oce ! ocean space and time domain 13 13 USE sbc_oce ! surface boundary condition 14 USE dynspg_oce ! surface pressure gradient variables15 14 USE phycst ! physical constants 16 15 USE fldread ! read input fields … … 112 111 IF(lwp) WRITE(numout,*) ' Inverse barometer added to OBC ssh data' 113 112 ENDIF 114 IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts ) & 115 CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' ) 113 !jc: stop below should rather be a warning 116 114 IF( ( ln_apr_obc ) .AND. .NOT. ln_apr_dyn ) & 117 115 CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY requires ln_apr_dyn=T' ) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
r5215 r5930 10 10 USE phycst 11 11 USE daymod 12 USE dynspg_oce13 12 USE tideini 14 13 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
r5215 r5930 10 10 USE phycst 11 11 USE daymod 12 USE dynspg_oce13 12 USE tide_mod 14 13 ! … … 96 95 & CALL ctl_stop('rdttideramp must be positive') 97 96 ! 98 IF( .NOT. lk_dynspg_ts ) CALL ctl_warn( 'sbc_tide : use of time splitting is recommended' )99 97 ! 100 98 ALLOCATE( ntide(nb_harmo) ) -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r5836 r5930 26 26 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 27 27 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 28 USE c1d ! 1D vertical configuration 28 29 ! 29 30 USE in_out_manager ! I/O manager … … 214 215 CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 215 216 ENDIF 216 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 217 IF( (ioptio /= 1).AND. (.NOT. lk_c1d ) ) & 218 CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 217 219 ! 218 220 IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 ) & ! Centered -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r5836 r5930 19 19 USE trd_oce ! trends: ocean variables 20 20 USE trdtra ! tracers trends 21 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient22 21 USE diaptr ! poleward transport diagnostics 23 22 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
r5836 r5930 21 21 USE trd_oce ! trends: ocean variables 22 22 USE trdtra ! tracers trends manager 23 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient24 23 USE sbcrnf ! river runoffs 25 24 USE diaptr ! poleward transport diagnostics -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r5836 r5930 20 20 USE trd_oce ! trends: ocean variables 21 21 USE trdtra ! trends manager: tracers 22 USE dynspg_oce ! surface pressure gradient variables23 22 USE diaptr ! poleward transport diagnostics 24 23 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5836 r5930 18 18 USE traadv_fct ! acces to routine interp_4th_cpt 19 19 USE trdtra ! trends manager: tracers 20 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient21 20 USE diaptr ! poleward transport diagnostics 22 21 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r5836 r5930 31 31 USE zdf_oce ! ocean vertical mixing 32 32 USE domvvl ! variable volume 33 USE dynspg_oce ! surface pressure gradient variables34 USE dynhpg ! hydrostatic pressure gradient35 33 USE trd_oce ! trends: ocean variables 36 34 USE trdtra ! trends manager: tracers … … 58 56 PUBLIC tra_nxt_vvl ! to be used in trcnxt 59 57 60 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg61 58 62 59 !! * Substitutions … … 90 87 !! 91 88 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 92 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)89 !! 93 90 !!---------------------------------------------------------------------- 94 91 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 105 102 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 106 103 IF(lwp) WRITE(numout,*) '~~~~~~~' 107 !108 rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp) ! Brown & Campana parameter for semi-implicit hpg109 104 ENDIF 110 105 111 106 ! Update after tracer on domain lateral boundaries 112 107 ! 108 ! 113 109 #if defined key_agrif 114 110 CALL Agrif_tra ! AGRIF zoom boundaries … … 181 177 !! 182 178 !! ** Method : - Apply a Asselin time filter on now fields. 183 !! - save in (ta,sa) an average over the three time levels184 !! which will be used to compute rdn and thus the semi-implicit185 !! hydrostatic pressure gradient (ln_dynhpg_imp = T)186 179 !! - swap tracer fields to prepare the next time_step. 187 !! This can be summurized for tempearture as:188 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T189 !! ztm = 0 otherwise190 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp)191 !! tb = tn + atfp*[ tb - 2 tn + ta ]192 !! tn = ta193 !! ta = ztm (NB: reset to 0 after eos_bn2 call)194 180 !! 195 181 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 196 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)182 !! 197 183 !!---------------------------------------------------------------------- 198 184 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 205 191 ! 206 192 INTEGER :: ji, jj, jk, jn ! dummy loop indices 207 LOGICAL :: ll_tra_hpg ! local logical208 193 REAL(wp) :: ztn, ztd ! local scalars 209 194 !!---------------------------------------------------------------------- … … 215 200 ENDIF 216 201 ! 217 IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg218 ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg219 ENDIF220 202 ! 221 203 DO jn = 1, kjpt … … 230 212 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 231 213 ! 232 IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average233 214 END DO 234 215 END DO … … 248 229 !! 249 230 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 250 !! - save in (ta,sa) a thickness weighted average over the three251 !! time levels which will be used to compute rdn and thus the semi-252 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)253 231 !! - swap tracer fields to prepare the next time_step. 254 232 !! This can be summurized for tempearture as: 255 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T256 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )257 !! ztm = 0 otherwise258 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )259 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] )260 !! tn = ta261 !! ta = zt (NB: reset to 0 after eos_bn2 call)262 233 !! 263 234 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 264 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T)235 !! 265 236 !!---------------------------------------------------------------------- 266 237 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 276 247 277 248 !! 278 LOGICAL :: ll_tra _hpg, ll_traqsr, ll_rnf, ll_isf ! local logical249 LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical 279 250 INTEGER :: ji, jj, jk, jn ! dummy loop indices 280 251 REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar … … 289 260 ! 290 261 IF( cdtype == 'TRA' ) THEN 291 ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg292 262 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 293 263 ll_rnf = ln_rnf ! active tracers case and river runoffs … … 298 268 END IF 299 269 ELSE 300 ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg301 270 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 302 271 ll_rnf = .FALSE. ! passive tracers or NO river runoffs … … 356 325 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 357 326 ! 358 IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only)359 ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )360 pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average361 ENDIF362 327 END DO 363 328 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r5836 r5930 18 18 USE zdf_oce ! ocean vertical physics variables 19 19 USE sbc_oce ! surface boundary condition: ocean 20 USE dynspg_oce21 20 ! 22 21 USE ldftra ! lateral diffusion: eddy diffusivity -
trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90
r5215 r5930 71 71 INTEGER, PUBLIC, PARAMETER :: jpdyn_bfri = 12 !: implicit bottom friction (ln_bfrimp=.TRUE.) 72 72 INTEGER, PUBLIC, PARAMETER :: jpdyn_ken = 13 !: use for calculation of KE 73 INTEGER, PUBLIC, PARAMETER :: jpdyn_spgflt = 14 !: filter contribution to surface pressure gradient (spg_flt)74 INTEGER, PUBLIC, PARAMETER :: jpdyn_spgexp = 15 !: explicit contribution to surface pressure gradient (spg_flt)75 73 ! 76 74 !!---------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90
r5215 r5930 113 113 CASE( jpdyn_spg ) ; CALL iom_put( "utrd_spg", putrd ) ! surface pressure gradient 114 114 CALL iom_put( "vtrd_spg", pvtrd ) 115 CASE( jpdyn_spgexp ); CALL iom_put( "utrd_spgexp", putrd ) ! surface pressure gradient (explicit)116 CALL iom_put( "vtrd_spgexp", pvtrd )117 CASE( jpdyn_spgflt ); CALL iom_put( "utrd_spgflt", putrd ) ! surface pressure gradient (filtered)118 CALL iom_put( "vtrd_spgflt", pvtrd )119 115 CASE( jpdyn_pvo ) ; CALL iom_put( "utrd_pvo", putrd ) ! planetary vorticity 120 116 CALL iom_put( "vtrd_pvo", pvtrd ) -
trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90
r5836 r5930 120 120 CASE( jpdyn_hpg ) ; CALL iom_put( "ketrd_hpg", zke ) ! hydrostatic pressure gradient 121 121 CASE( jpdyn_spg ) ; CALL iom_put( "ketrd_spg", zke ) ! surface pressure gradient 122 CASE( jpdyn_spgexp ); CALL iom_put( "ketrd_spgexp", zke ) ! surface pressure gradient (explicit)123 CASE( jpdyn_spgflt ); CALL iom_put( "ketrd_spgflt", zke ) ! surface pressure gradient (filter)124 122 CASE( jpdyn_pvo ) ; CALL iom_put( "ketrd_pvo", zke ) ! planetary vorticity 125 123 CASE( jpdyn_rvo ) ; CALL iom_put( "ketrd_rvo", zke ) ! relative vorticity (or metric term) -
trunk/NEMOGCM/NEMO/OPA_SRC/oce.F90
r5836 r5930 21 21 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ub , un , ua !: i-horizontal velocity [m/s] 22 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vb , vn , va !: j-horizontal velocity [m/s] 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ua_sv, va_sv !: Saved trends (time spliting) [m/s2]24 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wn !: vertical velocity [m/s] 25 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn !: horizontal divergence [s-1] … … 38 37 ! 39 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: spgu, spgv !: horizontal surface pressure gradient 39 40 !! Specific to split explicit free surface (allocated in dynspg_ts module): 41 ! 42 !! Time filtered arrays at baroclinic time step: 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv, vn_adv ! Advection vel. at "now" barocl. step 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b ! Half step fluxes (ln_bt_fw=T) 45 #if defined key_agrif 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b ! Half step time integrated fluxes 47 #endif 48 49 !! Arrays at barotropic time step: ! bef before ! before ! now ! after ! 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubb_e , ub_e , un_e , ua_e 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vbb_e , vb_e , vn_e , va_e 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e , sshb_e , sshn_e , ssha_e 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_e 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e 55 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e 40 57 41 58 !! interpolated gradient (only used in zps case) … … 85 102 ! 86 103 ALLOCATE( ub (jpi,jpj,jpk) , un (jpi,jpj,jpk) , ua(jpi,jpj,jpk) , & 87 & vb (jpi,jpj,jpk) , vn (jpi,jpj,jpk) , va(jpi,jpj,jpk) , & 88 & ua_sv(jpi,jpj,jpk) , va_sv(jpi,jpj,jpk) , & 104 & vb (jpi,jpj,jpk) , vn (jpi,jpj,jpk) , va(jpi,jpj,jpk) , & 89 105 & wn (jpi,jpj,jpk) , hdivn(jpi,jpj,jpk) , & 90 106 & tsb (jpi,jpj,jpk,jpts) , tsn (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) , & -
trunk/NEMOGCM/NEMO/OPA_SRC/step.F90
r5836 r5930 28 28 !! 3.7 ! 2014-10 (G. Madec) LDF simplication 29 29 !! - ! 2014-12 (G. Madec) remove KPP scheme 30 !! - ! 2015-11 (J. Chanut) free surface simplification 30 31 !!---------------------------------------------------------------------- 31 32 … … 179 180 180 181 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 181 ! Ocean dynamics : hdiv, ssh, e3, wn 182 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 182 ! Ocean dynamics : hdiv, ssh, e3, u, v, w 183 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 184 183 185 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 184 186 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 185 187 CALL wzv ( kstp ) ! now cross-level velocity 186 188 187 IF( lk_dynspg_ts ) THEN188 ! In case the time splitting case, update almost all momentum trends here:189 ! Note that the computation of vertical velocity above, hence "after" sea level190 ! is necessary to compute momentum advection for the rhs of barotropic loop:191 189 !!gm : why also here ???? 192 IF(ln_sto_eos )CALL sto_pts( tsn ) ! Random T/S fluctuations190 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 193 191 !!gm 194 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 195 192 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 193 194 !!jc: fs simplification 195 !!jc: lines below are useless if lk_vvl=T. Keep them here (which maintains a bug if lk_vvl=F and ln_zps=T, cf ticket #1636) 196 !! but ensures reproductible results 197 !! with previous versions using split-explicit free surface 196 198 IF( ln_zps .AND. .NOT. ln_isfcav) & ! Partial steps: bottom before horizontal gradient 197 199 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! of t, s, rd at the last ocean level … … 201 203 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 202 204 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 203 204 ua(:,:,:) = 0._wp ! set dynamics trends to zero 205 va(:,:,:) = 0._wp 206 IF( lk_asminc .AND. ln_asmiau .AND. & 207 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 208 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 209 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 210 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 211 CALL dyn_ldf ( kstp ) ! lateral mixing 212 #if defined key_agrif 213 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 214 #endif 215 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 216 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 217 218 ua_sv(:,:,:) = ua(:,:,:) ! Save trends (barotropic trend has been fully updated at this stage) 219 va_sv(:,:,:) = va(:,:,:) 220 221 CALL div_hor( kstp ) ! Horizontal divergence (2nd call in time-split case) 222 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 223 CALL wzv ( kstp ) ! now cross-level velocity 224 ENDIF 225 226 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 227 ! diagnostics and outputs (ua, va, tsa used as workspace) 228 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 229 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 230 IF( lk_diahth ) CALL dia_hth( kstp ) ! Thermocline depth (20 degres isotherm depth) 231 IF(.NOT.ln_cpl ) CALL dia_fwb( kstp ) ! Fresh water budget diagnostics 232 IF( lk_diadct ) CALL dia_dct( kstp ) ! Transports 233 IF( lk_diaar5 ) CALL dia_ar5( kstp ) ! ar5 diag 234 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 235 CALL dia_wri( kstp ) ! ocean model: outputs 236 ! 237 IF( ln_crs ) CALL crs_fld( kstp ) ! ocean model: online field coarsening & output 205 !!jc: fs simplification 206 207 ua(:,:,:) = 0._wp ! set dynamics trends to zero 208 va(:,:,:) = 0._wp 209 210 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) & 211 CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 212 IF( lk_bdy ) CALL bdy_dyn3d_dmp ( kstp ) ! bdy damping trends 213 #if defined key_agrif 214 IF(.NOT. Agrif_Root()) & 215 & CALL Agrif_Sponge_dyn ! momentum sponge 216 #endif 217 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 218 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 219 CALL dyn_ldf ( kstp ) ! lateral mixing 220 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 221 CALL dyn_spg ( kstp ) ! surface pressure gradient 222 223 ! With split-explicit free surface, since now transports have been updated and ssha as well 224 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 225 CALL div_hor ( kstp ) ! Horizontal divergence (2nd call in time-split case) 226 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 227 CALL wzv ( kstp ) ! now cross-level velocity 228 ENDIF 229 230 CALL dyn_bfr ( kstp ) ! bottom friction 231 CALL dyn_zdf ( kstp ) ! vertical diffusion 232 233 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 234 ! diagnostics and outputs 235 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 236 IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats 237 IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 238 IF(.NOT.ln_cpl ) CALL dia_fwb ( kstp ) ! Fresh water budget diagnostics 239 IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports 240 IF( lk_diaar5 ) CALL dia_ar5 ( kstp ) ! ar5 diag 241 IF( lk_diaharm ) CALL dia_harm ( kstp ) ! Tidal harmonic analysis 242 CALL dia_wri ( kstp ) ! ocean model: outputs 243 ! 244 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output 238 245 239 246 #if defined key_top … … 241 248 ! Passive Tracer Model 242 249 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 243 CALL trc_stp ( kstp )! time-stepping244 #endif 245 246 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 247 ! Active tracers (ua, va used as workspace)248 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 249 tsa(:,:,:,:) = 0._wp! set tracer trends to zero250 CALL trc_stp ( kstp ) ! time-stepping 251 #endif 252 253 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 254 ! Active tracers 255 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 256 tsa(:,:,:,:) = 0._wp ! set tracer trends to zero 250 257 251 258 IF( lk_asminc .AND. ln_asmiau .AND. & 252 & ln_trainc ) CALL tra_asm_inc( kstp ) ! apply tracer assimilation increment 253 CALL tra_sbc ( kstp ) ! surface boundary condition 254 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 255 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 256 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 257 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 258 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 259 CALL tra_adv ( kstp ) ! horizontal & vertical advection 260 CALL tra_ldf ( kstp ) ! lateral mixing 259 & ln_trainc ) CALL tra_asm_inc ( kstp ) ! apply tracer assimilation increment 260 CALL tra_sbc ( kstp ) ! surface boundary condition 261 IF( ln_traqsr ) CALL tra_qsr ( kstp ) ! penetrative solar radiation qsr 262 IF( ln_trabbc ) CALL tra_bbc ( kstp ) ! bottom heat flux 263 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 264 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 265 IF( lk_bdy ) CALL bdy_tra_dmp ( kstp ) ! bdy damping trends 266 #if defined key_agrif 267 IF(.NOT. Agrif_Root()) & 268 & CALL Agrif_Sponge_tra ! tracers sponge 269 #endif 270 CALL tra_adv ( kstp ) ! horizontal & vertical advection 271 CALL tra_ldf ( kstp ) ! lateral mixing 261 272 262 273 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 263 IF( ln_diaptr ) CALL dia_ptr! Poleward adv/ldf TRansports diagnostics274 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 264 275 !!gm 265 266 #if defined key_agrif 267 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra ! tracers sponge 268 #endif 269 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 270 271 IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg (time stepping then eos) 272 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 273 CALL tra_nxt( kstp ) ! tracer fields at next time step 274 !!gm : why again a call to sto_pts ??? 275 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 276 !!gm 277 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 278 IF( ln_zps .AND. .NOT. ln_isfcav) & 279 & CALL zps_hde ( kstp, jpts, tsa, gtsu, gtsv, & ! Partial steps: before horizontal gradient 280 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 281 IF( ln_zps .AND. ln_isfcav) & 282 & CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top/bottom cells 283 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 284 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 285 ELSE ! centered hpg (eos then time stepping) 286 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 287 !!gm : why again a call to sto_pts ??? 288 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 289 !!gm 290 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 291 IF( ln_zps .AND. .NOT. ln_isfcav) & 292 & CALL zps_hde ( kstp, jpts, tsn, gtsu, gtsv, & ! Partial steps: bottom before horizontal gradient 293 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 294 IF( ln_zps .AND. ln_isfcav) & 295 & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top/bottom cells 296 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 297 & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 298 ENDIF 299 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 300 CALL tra_nxt( kstp ) ! tracer fields at next time step 301 ENDIF 302 303 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 304 ! Dynamics (tsa used as workspace) 305 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 306 IF( lk_dynspg_ts ) THEN 307 ! revert to previously computed momentum tendencies 308 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 309 ua(:,:,:) = ua_sv(:,:,:) 310 va(:,:,:) = va_sv(:,:,:) 311 312 CALL dyn_bfr( kstp ) ! bottom friction 313 CALL dyn_zdf( kstp ) ! vertical diffusion 314 ELSE 315 ua(:,:,:) = 0._wp ! set dynamics trends to zero 316 va(:,:,:) = 0._wp 317 318 IF( lk_asminc .AND. ln_asmiau .AND. & 319 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 320 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 321 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 322 CALL dyn_adv( kstp ) ! advection (vector or flux form) 323 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 324 CALL dyn_ldf( kstp ) ! lateral mixing 325 #if defined key_agrif 326 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 327 #endif 328 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 329 CALL dyn_bfr( kstp ) ! bottom friction 330 CALL dyn_zdf( kstp ) ! vertical diffusion 331 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 332 ENDIF 333 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 334 335 CALL ssh_swp( kstp ) ! swap of sea surface height 336 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 276 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields 277 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 278 279 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 280 ! Set boundary conditions and Swap 281 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 282 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 283 !! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 284 !! If so: 285 !! (i) no need to call agrif update at initialization time 286 !! (ii) no need to update "before" fields 287 !! 288 !! Apart from creating new tra_swp/dyn_swp routines, this however: 289 !! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 290 !! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 291 !! e.g. a shift of the feedback interface inside child domain. 292 !! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 293 !! place. 294 !! 295 !!jc2: dynnxt must be the latest call. fse3t_b are indeed updated in that routine 296 CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap 297 CALL dyn_nxt ( kstp ) ! finalize (bcs) velocities at next time step and swap 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 337 300 ! 338 301 339 302 !!gm : This does not only concern the dynamics ==>>> add a new title 340 303 !!gm2: why ouput restart before AGRIF update? 341 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 304 !! 305 !!jc: That would be better, but see comment above 306 !! 307 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 342 308 343 309 #if defined key_agrif … … 345 311 ! AGRIF 346 312 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 347 CALL Agrif_Integrate_ChildGrids( stp ) 348 349 IF ( Agrif_NbStepint().EQ.0 ) THEN 350 CALL Agrif_Update_Tra() ! Update active tracers 351 CALL Agrif_Update_Dyn() ! Update momentum 352 ENDIF 353 #endif 354 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 355 IF( lk_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 313 CALL Agrif_Integrate_ChildGrids( stp ) 314 315 IF ( Agrif_NbStepint().EQ.0 ) THEN ! AGRIF Update 316 !!jc in fact update i useless at last time step, but do it for global diagnostics 317 CALL Agrif_Update_Tra() ! Update active tracers 318 CALL Agrif_Update_Dyn() ! Update momentum 319 ENDIF 320 #endif 321 IF( ln_diahsb ) CALL dia_hsb ( kstp ) ! - ML - global conservation diagnostics 322 IF( lk_diaobs ) CALL dia_obs ( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 356 323 357 324 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 358 325 ! Control 359 326 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 360 CALL stp_ctl( kstp, indic )361 IF( indic < 0 362 363 364 ENDIF 365 IF( kstp == nit000 327 CALL stp_ctl ( kstp, indic ) 328 IF( indic < 0 ) THEN 329 CALL ctl_stop( 'step: indic < 0' ) 330 CALL dia_wri_state( 'output.abort', kstp ) 331 ENDIF 332 IF( kstp == nit000 ) THEN 366 333 CALL iom_close( numror ) ! close input ocean restart file 367 334 IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce … … 374 341 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 375 342 !!gm why lk_oasis and not lk_cpl ???? 376 IF( lk_oasis 343 IF( lk_oasis ) CALL sbc_cpl_snd( kstp ) ! coupled mode : field exchanges 377 344 ! 378 345 #if defined key_iomput -
trunk/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r5836 r5930 40 40 USE dynldf ! lateral momentum diffusion (dyn_ldf routine) 41 41 USE dynzdf ! vertical diffusion (dyn_zdf routine) 42 USE dynspg_oce ! surface pressure gradient (dyn_spg routine)43 42 USE dynspg ! surface pressure gradient (dyn_spg routine) 44 43 -
trunk/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
r3294 r5930 16 16 USE oce ! ocean dynamics and tracers variables 17 17 USE dom_oce ! ocean space and time domain variables 18 USE sol_oce ! ocean space and time domain variables19 18 USE in_out_manager ! I/O manager 20 19 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 21 20 USE lib_mpp ! distributed memory computing 22 USE dynspg_oce ! pressure gradient schemes23 21 USE c1d ! 1D vertical configuration 24 22 … … 43 41 !! ** Method : - Save the time step in numstp 44 42 !! - Print it each 50 time steps 45 !! - Print solver statistics in numsol 46 !! - Stop the run IF problem for the solver ( indec < 0 ) 43 !! - Stop the run IF problem ( indic < 0 ) 47 44 !! 48 45 !! ** Actions : 'time.step' file containing the last ocean time-step … … 50 47 !!---------------------------------------------------------------------- 51 48 INTEGER, INTENT( in ) :: kt ! ocean time-step index 52 INTEGER, INTENT( inout ) :: kindic ! indicator of solver convergence49 INTEGER, INTENT( inout ) :: kindic ! error indicator 53 50 !! 54 51 INTEGER :: ji, jj, jk ! dummy loop indices … … 143 140 IF( lk_c1d ) RETURN ! No log file in case of 1D vertical configuration 144 141 145 ! log file (solver or ssh statistics) 146 ! -------- 147 IF( lk_dynspg_flt ) THEN ! elliptic solver statistics (if required) 148 ! 149 IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps ! Solver 150 ! 151 IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN ! create a abort file if problem found 152 IF(lwp) THEN 153 WRITE(numout,*) ' stpctl: the elliptic solver DO not converge or explode' 154 WRITE(numout,*) ' ====== ' 155 WRITE(numout,9200) kt, niter, res, sqrt(epsr)/eps 156 WRITE(numout,*) 157 WRITE(numout,*) ' stpctl: output of last fields' 158 WRITE(numout,*) ' ====== ' 159 ENDIF 160 ENDIF 161 ! 162 ELSE !* ssh statistics (and others...) 163 IF( kt == nit000 .AND. lwp ) THEN ! open ssh statistics file (put in solver.stat file) 164 CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 165 ENDIF 166 ! 167 zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) ) 168 IF( lk_mpp ) CALL mpp_sum( zssh2 ) ! sum over the global domain 169 ! 170 IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin ! ssh statistics 171 ! 142 ! log file (ssh statistics) 143 ! -------- !* ssh statistics (and others...) 144 IF( kt == nit000 .AND. lwp ) THEN ! open ssh statistics file (put in solver.stat file) 145 CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 172 146 ENDIF 147 ! 148 zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) ) 149 IF( lk_mpp ) CALL mpp_sum( zssh2 ) ! sum over the global domain 150 ! 151 IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin ! ssh statistics 152 ! 173 153 174 154 9200 FORMAT('it:', i8, ' iter:', i4, ' r: ',e16.10, ' b: ',e16.10 )
Note: See TracChangeset
for help on using the changeset viewer.