- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- Location:
- trunk/NEMO
- Files:
-
- 24 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/NST_SRC/agrif_opa_update.F90
r1300 r1438 80 80 81 81 82 CALL Agrif_ChildGrid_To_ParentGrid() 83 CALL recompute_diags( kt ) 84 CALL Agrif_ParentGrid_To_ChildGrid() 82 !Done in step 83 ! CALL Agrif_ChildGrid_To_ParentGrid() 84 ! CALL recompute_diags( kt ) 85 ! CALL Agrif_ParentGrid_To_ChildGrid() 85 86 86 87 #endif -
trunk/NEMO/OPA_SRC/DOM/dom_oce.F90
r1241 r1438 114 114 #if defined key_vvl 115 115 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .TRUE. !: variable grid flag 116 117 116 !! All coordinates 118 117 !! --------------- 119 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: 120 gdep3w_1 , & !: depth of T-points (sum of e3w) (m) 121 gdept_1, gdepw_1, & !: analytical depth at T-W points (m) 122 e3v_1 , e3f_1 , & !: analytical vertical scale factors at V--F 123 e3t_1 , e3u_1 , & !: T--U points (m) 124 e3vw_1 , & !: analytical vertical scale factors at VW-- 125 e3w_1 , e3uw_1 !: W--UW points (m) 118 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: gdep3w_1 !: depth of T-points (sum of e3w) (m) 119 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: gdept_1, gdepw_1 !: analytical depth at T-W points (m) 120 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3v_1 , e3f_1 !: analytical vertical scale factors at V--F 121 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3t_1 , e3u_1 !: T--U points (m) 122 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3vw_1 !: analytical vertical scale factors at VW-- 123 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e3w_1 , e3uw_1 !: W--UW points (m) 126 124 #else 127 125 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .FALSE. !: fixed grid flag … … 129 127 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 130 128 hur, hvr, & !: inverse of u and v-points ocean depth (1/m) 131 hu , hv !: depth at u- and v-points (meters) 129 hu , hv, & !: depth at u- and v-points (meters) 130 hu_0 , hv_0 !: refernce depth at u- and v-points (meters) 132 131 133 132 !! z-coordinate with full steps (also used in the other cases as reference z-coordinate) -
trunk/NEMO/OPA_SRC/DOM/domain.F90
r1335 r1438 4 4 !! Ocean initialization : domain initialization 5 5 !!============================================================================== 6 6 !! History : OPA ! 1990-10 (C. Levy - G. Madec) Original code 7 !! ! 1991-11 (G. Madec) 8 !! ! 1992-01 (M. Imbard) insert time step initialization 9 !! ! 1996-06 (G. Madec) generalized vertical coordinate 10 !! ! 1997-02 (G. Madec) creation of domwri.F 11 !! ! 2001-05 (E.Durand - G. Madec) insert closed sea 12 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 13 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 14 !!---------------------------------------------------------------------- 15 7 16 !!---------------------------------------------------------------------- 8 17 !! dom_init : initialize the space and time domain … … 10 19 !! dom_ctl : control print for the ocean domain 11 20 !!---------------------------------------------------------------------- 12 !! * Modules used13 21 USE oce ! 14 22 USE dom_oce ! ocean space and time domain … … 30 38 PRIVATE 31 39 32 !! * Routine accessibility 33 PUBLIC dom_init ! called by opa.F90 40 PUBLIC dom_init ! called by opa.F90 34 41 35 42 !! * Substitutions 36 43 # include "domzgr_substitute.h90" 37 !!---------------------------------------------------------------------- 38 !! OPA 9.0 , LOCEAN-IPSL (2005)44 !!------------------------------------------------------------------------- 45 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 39 46 !! $Id$ 40 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt41 !!---------------------------------------------------------------------- 47 !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 48 !!------------------------------------------------------------------------- 42 49 43 50 CONTAINS … … 58 65 !! - dom_stp: defined the model time step 59 66 !! - dom_wri: create the meshmask file if nmsh=1 60 !! 61 !! History : 62 !! ! 90-10 (C. Levy - G. Madec) Original code 63 !! ! 91-11 (G. Madec) 64 !! ! 92-01 (M. Imbard) insert time step initialization 65 !! ! 96-06 (G. Madec) generalized vertical coordinate 66 !! ! 97-02 (G. Madec) creation of domwri.F 67 !! ! 01-05 (E.Durand - G. Madec) insert closed sea 68 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 69 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 70 !!---------------------------------------------------------------------- 71 !! * Local declarations 67 !!---------------------------------------------------------------------- 72 68 INTEGER :: jk ! dummy loop argument 73 69 INTEGER :: iconf = 0 ! temporary integers … … 90 86 CALL dom_msk ! Masks 91 87 92 IF( lk_vvl ) CALL dom_vvl _ini! Vertical variable mesh88 IF( lk_vvl ) CALL dom_vvl ! Vertical variable mesh 93 89 94 90 ! Local depth or Inverse of the local depth of the water column at u- and v-points 95 91 ! ------------------------------ 96 92 ! Ocean depth at U- and V-points 97 hu(:,:) = 0. 98 hv(:,:) = 0. 99 93 hu(:,:) = 0.e0 94 hv(:,:) = 0.e0 100 95 DO jk = 1, jpk 101 96 hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) … … 105 100 hur(:,:) = fse3u(:,:,1) ! Lower bound : thickness of the first model level 106 101 hvr(:,:) = fse3v(:,:,1) 107 108 102 DO jk = 2, jpk ! Sum of the vertical scale factors 109 103 hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 110 104 hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 111 105 END DO 112 113 106 ! Compute and mask the inverse of the local depth 114 107 hur(:,:) = 1. / hur(:,:) * umask(:,:,1) 115 108 hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) 116 109 117 118 110 CALL dom_stp ! Time step 119 111 … … 121 113 122 114 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 123 115 ! 124 116 END SUBROUTINE dom_init 125 117 … … 134 126 !! - namdom namelist 135 127 !! - namcla namelist 136 !! 137 !! History : 138 !! 9.0 ! 03-08 (G. Madec) Original code 139 !!---------------------------------------------------------------------- 140 !! * Modules used 128 !!---------------------------------------------------------------------- 141 129 USE ioipsl 142 130 NAMELIST/namrun/ no , cexper, cn_ocerst_in, cn_ocerst_out, ln_rstart, nrstdt, & … … 156 144 ENDIF 157 145 158 ! Namelist namrun : parameters of the run 159 REWIND( numnam ) 146 REWIND( numnam ) ! Namelist namrun : parameters of the run 160 147 READ ( numnam, namrun ) 161 162 148 IF(lwp) THEN 163 149 WRITE(numout,*) ' Namelist namrun' … … 228 214 ENDIF 229 215 230 ! Namelist namdom : space/time domain (bathymetry, mesh, timestep) 231 REWIND( numnam ) 216 REWIND( numnam ) ! Namelist namdom : space/time domain (bathymetry, mesh, timestep) 232 217 READ ( numnam, namdom ) 233 218 … … 252 237 ENDIF 253 238 254 ! Default values255 239 n_cla = 0 256 257 ! Namelist cross land advection 258 REWIND( numnam ) 240 REWIND( numnam ) ! Namelist cross land advection 259 241 READ ( numnam, namcla ) 260 242 IF(lwp) THEN … … 264 246 ENDIF 265 247 266 IF( nbit_cmp == 1 .AND. n_cla /= 0 ) THEN 267 CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' ) 268 END IF 269 248 IF( nbit_cmp == 1 .AND. n_cla /= 0 ) CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' ) 249 ! 270 250 END SUBROUTINE dom_nam 271 251 … … 278 258 !! 279 259 !! ** Method : compute and print extrema of masked scale factors 280 !! 281 !! History : 282 !! 8.5 ! 02-08 (G. Madec) Original code 283 !!---------------------------------------------------------------------- 284 !! * Local declarations 260 !!---------------------------------------------------------------------- 285 261 INTEGER :: iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 286 262 INTEGER, DIMENSION(2) :: iloc ! 287 263 REAL(wp) :: ze1min, ze1max, ze2min, ze2max 288 264 !!---------------------------------------------------------------------- 289 290 ! Extrema of the scale factors291 265 292 266 IF(lwp)WRITE(numout,*) … … 325 299 WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 326 300 ENDIF 327 301 ! 328 302 END SUBROUTINE dom_ctl 329 303 -
trunk/NEMO/OPA_SRC/DOM/domvvl.F90
r1146 r1438 4 4 !! Ocean : 5 5 !!====================================================================== 6 !! History : 9.0 !06-06 (B. Levier, L. Marie) original code7 !! " ! 07-07 (D. Storkey) Bug fixes and code for BDY option.6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 8 !!---------------------------------------------------------------------- 9 9 #if defined key_vvl 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_vvl' variable volume 12 12 !!---------------------------------------------------------------------- 13 !! dom_vvl : defined coefficients to distribute ssh on each layers 13 14 !!---------------------------------------------------------------------- 14 !! dom_vvl : defined scale factors & depths at each time step15 !! dom_vvl_ini : defined coefficients to distribute ssh on each layers16 !! dom_vvl_ssh : defined the ocean sea level at each time step17 !!----------------------------------------------------------------------18 !! * Modules used19 15 USE oce ! ocean dynamics and tracers 20 16 USE dom_oce ! ocean space and time domain 21 17 USE sbc_oce ! surface boundary condition: ocean 22 USE dynspg_oce ! surface pressure gradient variables23 18 USE phycst ! physical constants 24 19 USE in_out_manager ! I/O manager 25 20 USE lib_mpp ! distributed memory computing library 26 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 USE bdy_oce ! unstructured open boundary conditions28 22 29 23 IMPLICIT NONE 30 24 PRIVATE 31 25 32 !! * Routine accessibility 33 PUBLIC dom_vvl_ini ! called by dom_init.F90 34 PUBLIC dom_vvl_ssh ! called by trazdf.F90 35 PUBLIC dom_vvl ! called by istate.F90 and step.F90 36 PUBLIC dom_vvl_sf_ini ! 37 PUBLIC dom_vvl_sf ! 26 PUBLIC dom_vvl ! called by domain.F90 38 27 39 !! * Module variables40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: 41 mut, muu, muv, muf !:28 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ee_t, ee_u, ee_v, ee_f !: ??? 29 30 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: mut, muu, muv, muf !: ??? 42 31 43 32 REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time-step, = 2 rdttra … … 48 37 # include "vectopt_loop_substitute.h90" 49 38 !!---------------------------------------------------------------------- 50 !! OPA 9.0 , LOCEAN-IPSL (2005)39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 51 40 !! $Id$ 52 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 55 44 CONTAINS 56 45 57 #if defined key_vvl 58 59 SUBROUTINE dom_vvl_ini 46 SUBROUTINE dom_vvl 60 47 !!---------------------------------------------------------------------- 61 !! *** ROUTINE dom_vvl _ini***48 !! *** ROUTINE dom_vvl *** 62 49 !! 63 50 !! ** Purpose : compute coefficients muX at T-U-V-F points to spread … … 65 52 !! 66 53 !!---------------------------------------------------------------------- 67 INTEGER :: ji, jj, jk , zpk68 REAL(wp) , DIMENSION(jpi,jpj) :: zmut, zmuu, zmuv, zmuf ! 2D workspace54 INTEGER :: ji, jj, jk 55 REAL(wp) :: zcoefu, zcoefv, zcoeff 69 56 !!---------------------------------------------------------------------- 70 57 71 58 IF(lwp) THEN 72 59 WRITE(numout,*) 73 WRITE(numout,*) 'dom_vvl _ini: Variable volume activated'74 WRITE(numout,*) '~~~~~~~~ ~~~compute coef. used to spread ssh over each layers'60 WRITE(numout,*) 'dom_vvl : Variable volume activated' 61 WRITE(numout,*) '~~~~~~~~ compute coef. used to spread ssh over each layers' 75 62 ENDIF 76 77 IF( ln_zps ) CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl')78 63 79 64 #if defined key_zco || defined key_dynspg_rl … … 81 66 #endif 82 67 83 fs vdept(:,:,:) = gdept (:,:,:)84 fs vdepw(:,:,:) = gdepw (:,:,:)85 fs vde3w(:,:,:) = gdep3w(:,:,:)86 fs ve3t(:,:,:) = e3t (:,:,:)87 fs ve3u(:,:,:) = e3u (:,:,:)88 fs ve3v(:,:,:) = e3v (:,:,:)89 fs ve3f(:,:,:) = e3f (:,:,:)90 fs ve3w(:,:,:) = e3w (:,:,:)91 fs ve3uw(:,:,:) = e3uw (:,:,:)92 fs ve3vw(:,:,:) = e3vw (:,:,:)68 fsdept(:,:,:) = gdept (:,:,:) 69 fsdepw(:,:,:) = gdepw (:,:,:) 70 fsde3w(:,:,:) = gdep3w(:,:,:) 71 fse3t (:,:,:) = e3t (:,:,:) 72 fse3u (:,:,:) = e3u (:,:,:) 73 fse3v (:,:,:) = e3v (:,:,:) 74 fse3f (:,:,:) = e3f (:,:,:) 75 fse3w (:,:,:) = e3w (:,:,:) 76 fse3uw(:,:,:) = e3uw (:,:,:) 77 fse3vw(:,:,:) = e3vw (:,:,:) 93 78 94 79 ! mu computation 95 zmut(:,:) = 0.e0 96 zmuu(:,:) = 0.e0 97 zmuv(:,:) = 0.e0 98 zmuf(:,:) = 0.e0 99 mut (:,:,:) = 0.e0 100 muu (:,:,:) = 0.e0 101 muv (:,:,:) = 0.e0 102 muf (:,:,:) = 0.e0 80 ! -------------- 81 ! define ee_t, u, v and f as in sigma coordinate (ee_t = 1/ht, ...) 82 ee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 83 ee_u(:,:) = fse3u_0(:,:,1) 84 ee_v(:,:) = fse3v_0(:,:,1) 85 ee_f(:,:) = fse3f_0(:,:,1) 86 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 87 ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 88 ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 89 ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 90 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 91 ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 92 END DO 93 END DO 94 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 95 ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 96 ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 97 ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 98 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 99 ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 100 END DO 101 CALL lbc_lnk( ee_f, 'F', 1. ) ! lateral boundary condition on ee_f 102 ! 103 DO jk = 1, jpk 104 mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk) ! at T levels 105 muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk) ! at T levels 106 muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk) ! at T levels 107 END DO 108 DO jk = 1, jpk 109 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 110 muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 111 END DO 112 muf(:,jpj,jk) = 0.e0 113 END DO 114 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition on ee_f 103 115 104 DO jj = 1, jpj 105 DO ji = 1, jpi 106 zpk = mbathy(ji,jj) - 1 107 DO jk = 1, zpk 108 zmut(ji,jj) = zmut(ji,jj) + fsve3t(ji,jj,jk) * SUM( fsve3t(ji,jj,jk:zpk) ) 109 zmuu(ji,jj) = zmuu(ji,jj) + fsve3u(ji,jj,jk) * SUM( fsve3u(ji,jj,jk:zpk) ) 110 zmuv(ji,jj) = zmuv(ji,jj) + fsve3v(ji,jj,jk) * SUM( fsve3v(ji,jj,jk:zpk) ) 111 zmuf(ji,jj) = zmuf(ji,jj) + fsve3f(ji,jj,jk) * SUM( fsve3f(ji,jj,jk:zpk) ) 112 END DO 116 117 ! Reference ocean depth at U- and V-points 118 hu_0(:,:) = 0.e0 119 hv_0(:,:) = 0.e0 120 DO jk = 1, jpk 121 hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 122 hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 123 END DO 124 125 ! before and now Sea Surface Height at u-, v-, f-points 126 DO jj = 1, jpjm1 127 DO ji = 1, jpim1 128 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 129 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 130 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 131 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 132 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 133 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 134 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 135 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) & 136 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 137 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 138 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 139 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 140 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 141 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 142 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 113 143 END DO 114 144 END DO 115 116 DO jj = 1, jpj117 DO ji = 1, jpi118 zpk = mbathy(ji,jj) - 1119 DO jk = 1, zpk120 #if defined key_sigma_vvl121 mut(ji,jj,jk) = 1./SUM( fsve3t(ji,jj,1:zpk) )122 muu(ji,jj,jk) = 1./SUM( fsve3u(ji,jj,1:zpk) )123 muv(ji,jj,jk) = 1./SUM( fsve3v(ji,jj,1:zpk) )124 muf(ji,jj,jk) = 1./SUM( fsve3f(ji,jj,1:zpk) )125 #else126 mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj)127 muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj)128 muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj)129 muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj)130 #endif131 END DO132 END DO133 END DO134 135 136 END SUBROUTINE dom_vvl_ini137 138 139 140 SUBROUTINE dom_vvl141 !!----------------------------------------------------------------------142 !! *** ROUTINE dom_vvl ***143 !!144 !! ** Purpose : compute ssh at U-V-F points, T-W scale factors and local145 !! depths at each time step.146 !!----------------------------------------------------------------------147 !! * Local declarations148 INTEGER :: ji, jj, jk ! dummy loop indices149 REAL(wp), DIMENSION(jpi,jpj) :: zsshf ! 2D workspace150 !!----------------------------------------------------------------------151 152 ! Sea Surface Height at u-v-fpoints153 DO jj = 1, jpjm1154 DO ji = 1,jpim1155 sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &156 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &157 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )158 !159 sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &160 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &161 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )162 !163 zsshf(ji,jj) = 0.25 * fmask(ji,jj,1) &164 & * ( sshn(ji ,jj) + sshn(ji ,jj+1) &165 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) )166 END DO167 END DO168 169 145 ! Boundaries conditions 170 CALL lbc_lnk( sshu, 'U', 1. ) 171 CALL lbc_lnk( sshv, 'V', 1. ) 172 CALL lbc_lnk( zsshf, 'F', 1. ) 173 174 ! Scale factors at T levels 175 DO jk = 1, jpkm1 176 fse3t(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshn(:,:) * mut(:,:,jk) ) 177 fse3u(:,:,jk) = fsve3u(:,:,jk) * ( 1 + sshu(:,:) * muu(:,:,jk) ) 178 fse3v(:,:,jk) = fsve3v(:,:,jk) * ( 1 + sshv(:,:) * muv(:,:,jk) ) 179 fse3f(:,:,jk) = fsve3f(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) ) 180 END DO 181 182 ! Scale factors at W levels 183 fse3w (:,:,1) = fse3t(:,:,1) 184 fse3uw(:,:,1) = fse3u(:,:,1) 185 fse3vw(:,:,1) = fse3v(:,:,1) 186 DO jk = 2, jpk 187 fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) ) 188 fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) ) 189 fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) ) 190 END DO 191 192 ! T and W points depth 193 fsdept(:,:,1) = 0.5 * fse3w(:,:,1) 194 fsdepw(:,:,1) = 0.e0 195 fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:) 196 DO jk = 2, jpk 197 fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk) 198 fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1) 199 fsde3w(:,:,jk) = fsdept(:,:,jk ) - sshn(:,:) 200 END DO 201 202 IF( MINVAL(fse3t(:,:,:)) < 0.0 ) THEN 203 CALL ctl_stop('E R R O R : dom_vvl','Level thickness fse3t less than zero.') 204 nstop = nstop+1 205 ENDIF 206 207 ! Local depth or Inverse of the local depth of the water column at u- and v-points 208 ! ------------------------------ 209 210 ! Ocean depth at U- and V-points 211 hu(:,:) = 0.e0 ; hv(:,:) = 0.e0 212 213 DO jk = 1, jpk 214 hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 215 hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 216 END DO 217 218 219 ! Inverse of the local depth 220 hur(:,:) = fse3u(:,:,1) ! Lower bound : thickness of the first model level 221 hvr(:,:) = fse3v(:,:,1) 222 223 DO jk = 2, jpk ! Sum of the vertical scale factors 224 hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk) 225 hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) 226 END DO 227 228 ! Compute and mask the inverse of the local depth 229 hur(:,:) = 1. / hur(:,:) * umask(:,:,1) 230 hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1) 231 146 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 147 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 148 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 149 ! 232 150 END SUBROUTINE dom_vvl 233 151 234 235 236 SUBROUTINE dom_vvl_ssh( kt )237 !!----------------------------------------------------------------------238 !! *** ROUTINE dom_vvl_ssh ***239 !!240 !! ** Purpose : compute the ssha field used in tra_zdf, dynnxt, tranxt241 !! and dynspg routines242 !!----------------------------------------------------------------------243 !! * Arguments244 INTEGER, INTENT( in ) :: kt ! time step245 246 !! * Local declarations247 INTEGER :: ji, jj, jk ! dummy loop indices248 REAL(wp) :: z2dt, zraur ! temporary scalars249 REAL(wp), DIMENSION(jpi,jpj) :: zun, zvn, zhdiv ! 2D workspace250 !!----------------------------------------------------------------------251 252 !! Sea surface height after (ssha) in variable volume case253 !! --------------------------====-----===============-----254 !! ssha is needed for the tracer time-stepping (trazdf_imp or255 !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the256 !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv).257 !! un, vn (or un_b and vn_b) and emp are needed, so it must be258 !! done after the call of oce_sbc.259 260 IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000261 r2dt(:) = rdttra(:) ! = rdtra (restarting with Euler time stepping)262 ELSEIF( kt <= nit000 + 1) THEN ! at nit000 or nit000+1263 r2dt(:) = 2. * rdttra(:) ! = 2 rdttra (leapfrog)264 ENDIF265 266 z2dt = r2dt(1)267 zraur = 1. / rauw268 269 ! Vertically integrated quantities270 ! --------------------------------271 IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN272 zun(:,:) = 0.e0 ; zvn(:,:) = 0.e0273 !274 DO jk = 1, jpkm1275 ! ! Vertically integrated transports (now)276 zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)277 zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)278 END DO279 ELSE ! lk_dynspg_ts is true280 zun (:,:) = un_b(:,:) ; zvn (:,:) = vn_b(:,:)281 ENDIF282 283 ! Horizontal divergence of barotropic transports284 !--------------------------------------------------285 DO jj = 2, jpjm1286 DO ji = 2, jpim1 ! vector opt.287 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) &288 & - e2u(ji-1,jj ) * zun(ji-1,jj ) &289 & + e1v(ji ,jj ) * zvn(ji ,jj ) &290 & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) &291 & / ( e1t(ji,jj) * e2t(ji,jj) )292 END DO293 END DO294 295 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts )296 ! open boundaries (div must be zero behind the open boundary)297 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column298 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east299 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west300 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north301 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south302 #endif303 304 #if defined key_bdy305 zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:)306 #endif307 308 CALL lbc_lnk( zhdiv, 'T', 1. )309 310 ! Sea surface elevation time stepping311 ! -----------------------------------312 IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN313 ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)314 ELSE ! lk_dynspg_ts is true315 ssha(:,:) = ( sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)316 ENDIF317 318 319 END SUBROUTINE dom_vvl_ssh320 321 SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 )322 !!----------------------------------------------------------------------323 !! *** ROUTINE dom_vvl_sf ***324 !!325 !! ** Purpose : compute vertical scale factor at each time step326 !!----------------------------------------------------------------------327 !! * Arguments328 CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type329 REAL(wp), DIMENSION(jpi,jpj) , INTENT( in ) :: zssh ! 2D workspace330 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: sfe3 ! 3D workspace331 332 !! * Local declarations333 INTEGER :: jk ! dummy loop indices334 !!----------------------------------------------------------------------335 336 SELECT CASE (gridp)337 338 CASE ('U')339 !340 DO jk = 1, jpk341 sfe3(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zssh(:,:) * muu(:,:,jk) )342 ENDDO343 344 CASE ('V')345 !346 DO jk = 1, jpk347 sfe3(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zssh(:,:) * muv(:,:,jk) )348 ENDDO349 350 CASE ('T')351 !352 DO jk = 1, jpk353 sfe3(:,:,jk) = fsve3t(:,:,jk) * ( 1 + zssh(:,:) * mut(:,:,jk) )354 ENDDO355 356 END SELECT357 358 END SUBROUTINE dom_vvl_sf359 360 SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini )361 !!----------------------------------------------------------------------362 !! *** ROUTINE dom_vvl_sf_ini ***363 !!364 !! ** Purpose : affect the appropriate vertical scale factor. It is done365 !! to avoid compiling error when using key_zco cpp key366 !!----------------------------------------------------------------------367 !! * Arguments368 CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type369 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini ! 3D workspace370 !!----------------------------------------------------------------------371 372 SELECT CASE (gridp)373 374 CASE ('U')375 !376 sfe3ini(:,:,:) = fse3u(:,:,:)377 378 CASE ('V')379 !380 sfe3ini(:,:,:) = fse3v(:,:,:)381 382 CASE ('T')383 !384 sfe3ini(:,:,:) = fse3t(:,:,:)385 386 END SELECT387 388 END SUBROUTINE dom_vvl_sf_ini389 152 #else 390 153 !!---------------------------------------------------------------------- 391 154 !! Default option : Empty routine 392 155 !!---------------------------------------------------------------------- 393 SUBROUTINE dom_vvl_ini 394 END SUBROUTINE dom_vvl_ini 156 CONTAINS 395 157 SUBROUTINE dom_vvl 396 158 END SUBROUTINE dom_vvl 397 SUBROUTINE dom_vvl_ssh( kt )398 !! * Arguments399 INTEGER, INTENT( in ) :: kt ! time step400 WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt401 END SUBROUTINE dom_vvl_ssh402 SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 )403 !! * Arguments404 CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type405 REAL(wp), DIMENSION(jpi,jpj) , INTENT( in ) :: zssh ! 2D workspace406 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3 ! 3D workspace407 sfe3(:,:,:) = 0.e0408 WRITE(*,*) 'sfe3: You should not have seen this print! error?', gridp409 WRITE(*,*) 'sfe3: You should not have seen this print! error?', zssh(1,1)410 END SUBROUTINE dom_vvl_sf411 SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini )412 !! * Arguments413 CHARACTER(LEN=1) , INTENT( in ) :: gridp ! grid point type414 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini ! 3D workspace415 sfe3ini(:,:,:) = 0.e0416 WRITE(*,*) 'sfe3ini: You should not have seen this print! error?', gridp417 END SUBROUTINE dom_vvl_sf_ini418 159 #endif 419 160 -
trunk/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
r1156 r1438 5 5 !! factors depending on the vertical coord. used, using CPP macro. 6 6 !!---------------------------------------------------------------------- 7 !! History : 1.0 ! 2005-10 (A. Beckmann, G. Madec) generalisation to all coord. 8 !! 3.1 ! 2009-02 (G. Madec, M. Leclair) pure z* coordinate 7 9 !!---------------------------------------------------------------------- 8 !! OPA 9.0 , LOCEAN-IPSL (2005) 10 #if defined key_zco 11 ! reference for pure z-coordinate (1D - no i,j and time dependency) 12 # define fsdept_0(i,j,k) gdept_0(k) 13 # define fsdepw_0(i,j,k) gdepw_0(k) 14 # define fsde3w_0(i,j,k) gdepw_0(k) 15 # define fse3t_0(i,j,k) e3t_0(k) 16 # define fse3u_0(i,j,k) e3t_0(k) 17 # define fse3v_0(i,j,k) e3t_0(k) 18 # define fse3f_0(i,j,k) e3t_0(k) 19 # define fse3w_0(i,j,k) e3w_0(k) 20 # define fse3uw_0(i,j,k) e3w_0(k) 21 # define fse3vw_0(i,j,k) e3w_0(k) 22 #else 23 ! reference for s- or zps-coordinate (3D no time dependency) 24 # define fsdept_0(i,j,k) gdept(i,j,k) 25 # define fsdepw_0(i,j,k) gdepw(i,j,k) 26 # define fsde3w_0(i,j,k) gdep3w(i,j,k) 27 # define fse3t_0(i,j,k) e3t(i,j,k) 28 # define fse3u_0(i,j,k) e3u(i,j,k) 29 # define fse3v_0(i,j,k) e3v(i,j,k) 30 # define fse3f_0(i,j,k) e3f(i,j,k) 31 # define fse3w_0(i,j,k) e3w(i,j,k) 32 # define fse3uw_0(i,j,k) e3uw(i,j,k) 33 # define fse3vw_0(i,j,k) e3vw(i,j,k) 34 #endif 35 #if defined key_vvl 36 ! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1) 37 # define fsdept(i,j,k) gdept_1(i,j,k) 38 # define fsdepw(i,j,k) gdepw_1(i,j,k) 39 # define fsde3w(i,j,k) gdep3w_1(i,j,k) 40 # define fse3t(i,j,k) e3t_1(i,j,k) 41 # define fse3u(i,j,k) e3u_1(i,j,k) 42 # define fse3v(i,j,k) e3v_1(i,j,k) 43 # define fse3f(i,j,k) e3f_1(i,j,k) 44 # define fse3w(i,j,k) e3w_1(i,j,k) 45 # define fse3uw(i,j,k) e3uw_1(i,j,k) 46 # define fse3vw(i,j,k) e3vw_1(i,j,k) 47 48 # define fsdept_b(i,j,k) (fsdept_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 49 # define fsdepw_b(i,j,k) (fsdepw_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 50 # define fsde3w_b(i,j,k) (fsde3w_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))-sshb(i,j)) 51 # define fse3t_b(i,j,k) (fse3t_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 52 # define fse3u_b(i,j,k) (fse3u_0(i,j,k)*(1+sshu_b(i,j)*muu(i,j,k))) 53 # define fse3v_b(i,j,k) (fse3v_0(i,j,k)*(1+sshv_b(i,j)*muv(i,j,k))) 54 # define fse3f_b(i,j,k) (fse3f_0(i,j,k)*(1+sshf_b(i,j)*muf(i,j,k))) 55 # define fse3w_b(i,j,k) (fse3w_0(i,j,k)*(1+sshb(i,j)*mut(i,j,k))) 56 # define fse3uw_b(i,j,k) (fse3uw_0(i,j,k)*(1+sshu_b(i,j)*muu(i,j,k))) 57 # define fse3vw_b(i,j,k) (fse3vw_0(i,j,k)*(1+sshv_b(i,j)*muv(i,j,k))) 58 59 # define fsdept_n(i,j,k) (fsdept_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 60 # define fsdepw_n(i,j,k) (fsdepw_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 61 # define fsde3w_n(i,j,k) (fsde3w_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))-sshn(i,j)) 62 # define fse3t_n(i,j,k) (fse3t_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 63 # define fse3u_n(i,j,k) (fse3u_0(i,j,k)*(1+sshu_n(i,j)*muu(i,j,k))) 64 # define fse3v_n(i,j,k) (fse3v_0(i,j,k)*(1+sshv_n(i,j)*muv(i,j,k))) 65 # define fse3f_n(i,j,k) (fse3f_0(i,j,k)*(1+sshf_n(i,j)*muf(i,j,k))) 66 # define fse3w_n(i,j,k) (fse3w_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k))) 67 # define fse3uw_n(i,j,k) (fse3uw_0(i,j,k)*(1+sshu_n(i,j)*muu(i,j,k))) 68 # define fse3vw_n(i,j,k) (fse3vw_0(i,j,k)*(1+sshv_n(i,j)*muv(i,j,k))) 69 70 # define fsdept_a(i,j,k) (fsdept_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 71 # define fsdepw_a(i,j,k) (fsdepw_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 72 # define fsde3w_a(i,j,k) (fsde3w_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))-ssha(i,j)) 73 # define fse3t_a(i,j,k) (fse3t_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 74 # define fse3u_a(i,j,k) (fse3u_0(i,j,k)*(1+sshu_a(i,j)*muu(i,j,k))) 75 # define fse3v_a(i,j,k) (fse3v_0(i,j,k)*(1+sshv_a(i,j)*muv(i,j,k))) 76 # define fse3f_a(i,j,k) (fse3f_0(i,j,k)*(1+sshf_a(i,j)*muf(i,j,k))) 77 # define fse3w_a(i,j,k) (fse3w_0(i,j,k)*(1+ssha(i,j)*mut(i,j,k))) 78 # define fse3uw_a(i,j,k) (fse3uw_0(i,j,k)*(1+sshu_a(i,j)*muu(i,j,k))) 79 # define fse3vw_a(i,j,k) (fse3vw_0(i,j,k)*(1+sshv_a(i,j)*muv(i,j,k))) 80 81 #else 82 ! z- or s-coordinate (1D or 3D + no time dependency) use reference in all cases 83 # define fsdept(i,j,k) fsdept_0(i,j,k) 84 # define fsdepw(i,j,k) fsdepw_0(i,j,k) 85 # define fsde3w(i,j,k) fsde3w_0(i,j,k) 86 # define fse3t(i,j,k) fse3t_0(i,j,k) 87 # define fse3u(i,j,k) fse3u_0(i,j,k) 88 # define fse3v(i,j,k) fse3v_0(i,j,k) 89 # define fse3f(i,j,k) fse3f_0(i,j,k) 90 # define fse3w(i,j,k) fse3w_0(i,j,k) 91 # define fse3uw(i,j,k) fse3uw_0(i,j,k) 92 # define fse3vw(i,j,k) fse3vw_0(i,j,k) 93 94 # define fsdept_b(i,j,k) fsdept_0(i,j,k) 95 # define fsdepw_b(i,j,k) fsdepw_0(i,j,k) 96 # define fsde3w_b(i,j,k) fsde3w_0(i,j,k) 97 # define fse3t_b(i,j,k) fse3t_0(i,j,k) 98 # define fse3u_b(i,j,k) fse3u_0(i,j,k) 99 # define fse3v_b(i,j,k) fse3v_0(i,j,k) 100 # define fse3f_b(i,j,k) fse3f_0(i,j,k) 101 # define fse3w_b(i,j,k) fse3w_0(i,j,k) 102 # define fse3uw_b(i,j,k) fse3uw_0(i,j,k) 103 # define fse3vw_b(i,j,k) fse3vw_0(i,j,k) 104 105 # define fsdept_n(i,j,k) fsdept_0(i,j,k) 106 # define fsdepw_n(i,j,k) fsdepw_0(i,j,k) 107 # define fsde3w_n(i,j,k) fsde3w_0(i,j,k) 108 # define fse3t_n(i,j,k) fse3t_0(i,j,k) 109 # define fse3u_n(i,j,k) fse3u_0(i,j,k) 110 # define fse3v_n(i,j,k) fse3v_0(i,j,k) 111 # define fse3f_n(i,j,k) fse3f_0(i,j,k) 112 # define fse3w_n(i,j,k) fse3w_0(i,j,k) 113 # define fse3uw_n(i,j,k) fse3uw_0(i,j,k) 114 # define fse3vw_n(i,j,k) fse3vw_0(i,j,k) 115 116 # define fsdept_a(i,j,k) fsdept_0(i,j,k) 117 # define fsdepw_a(i,j,k) fsdepw_0(i,j,k) 118 # define fsde3w_a(i,j,k) fsde3w_0(i,j,k) 119 # define fse3t_a(i,j,k) fse3t_0(i,j,k) 120 # define fse3u_a(i,j,k) fse3u_0(i,j,k) 121 # define fse3v_a(i,j,k) fse3v_0(i,j,k) 122 # define fse3f_a(i,j,k) fse3f_0(i,j,k) 123 # define fse3w_a(i,j,k) fse3w_0(i,j,k) 124 # define fse3uw_a(i,j,k) fse3uw_0(i,j,k) 125 # define fse3vw_a(i,j,k) fse3vw_0(i,j,k) 126 #endif 127 !!---------------------------------------------------------------------- 128 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 9 129 !! $Id$ 10 130 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 11 !!12 !! History :13 !! 9.0 ! 05-10 (A. Beckmann, G. Madec) generalisation to all coord.14 131 !!---------------------------------------------------------------------- 15 #if defined key_zco16 # define fsdept(i,j,k) gdept_0(k)17 18 # define fsdepw(i,j,k) gdepw_0(k)19 # define fsde3w(i,j,k) gdepw_0(k)20 21 # define fse3t(i,j,k) e3t_0(k)22 # define fse3u(i,j,k) e3t_0(k)23 # define fse3v(i,j,k) e3t_0(k)24 # define fse3f(i,j,k) e3t_0(k)25 26 # define fse3w(i,j,k) e3w_0(k)27 # define fse3uw(i,j,k) e3w_0(k)28 # define fse3vw(i,j,k) e3w_0(k)29 #else30 # define fsdept(i,j,k) gdept(i,j,k)31 32 # define fsdepw(i,j,k) gdepw(i,j,k)33 # define fsde3w(i,j,k) gdep3w(i,j,k)34 35 # define fse3t(i,j,k) e3t(i,j,k)36 # define fse3u(i,j,k) e3u(i,j,k)37 # define fse3v(i,j,k) e3v(i,j,k)38 # define fse3f(i,j,k) e3f(i,j,k)39 40 # define fse3w(i,j,k) e3w(i,j,k)41 # define fse3uw(i,j,k) e3uw(i,j,k)42 # define fse3vw(i,j,k) e3vw(i,j,k)43 #endif44 45 #if defined key_vvl46 # define fsvdept(i,j,k) gdept_1(i,j,k)47 48 # define fsvdepw(i,j,k) gdepw_1(i,j,k)49 # define fsvde3w(i,j,k) gdep3w_1(i,j,k)50 51 # define fsve3t(i,j,k) e3t_1(i,j,k)52 # define fsve3u(i,j,k) e3u_1(i,j,k)53 # define fsve3v(i,j,k) e3v_1(i,j,k)54 # define fsve3f(i,j,k) e3f_1(i,j,k)55 56 # define fsve3w(i,j,k) e3w_1(i,j,k)57 # define fsve3uw(i,j,k) e3uw_1(i,j,k)58 # define fsve3vw(i,j,k) e3vw_1(i,j,k)59 #else60 # define fsvdept(i,j,k) fsdept(i,j,k)61 62 # define fsvdepw(i,j,k) fsdepw(i,j,k)63 # define fsvde3w(i,j,k) fsde3w(i,j,k)64 65 # define fsve3t(i,j,k) fse3t(i,j,k)66 # define fsve3u(i,j,k) fse3u(i,j,k)67 # define fsve3v(i,j,k) fse3v(i,j,k)68 # define fsve3f(i,j,k) fse3f(i,j,k)69 70 # define fsve3w(i,j,k) fse3w(i,j,k)71 # define fsve3uw(i,j,k) fse3uw(i,j,k)72 # define fsve3vw(i,j,k) fse3vw(i,j,k)73 #endif74 -
trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
r1152 r1438 4 4 !! Ocean dynamics: time stepping 5 5 !!====================================================================== 6 6 !!====================================================================== 7 !! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code 8 !! ! 1990-10 (C. Levy, G. Madec) 9 !! 7.0 ! 1993-03 (M. Guyon) symetrical conditions 10 !! 8.0 ! 1997-02 (G. Madec & M. Imbard) opa, release 8.0 11 !! 8.2 ! 1997-04 (A. Weaver) Euler forward step 12 !! - ! 1997-06 (G. Madec) lateral boudary cond., lbc routine 13 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 14 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 15 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 16 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 17 !! 3.2 ! 2009-04 (G. Madec, R.Benshila)) re-introduce the vvl option 18 !!---------------------------------------------------------------------- 19 7 20 !!---------------------------------------------------------------------- 8 21 !! dyn_nxt : update the horizontal velocity from the momentum trend 9 22 !!---------------------------------------------------------------------- 10 !! * Modules used11 23 USE oce ! ocean dynamics and tracers 12 24 USE dom_oce ! ocean space and time domain … … 29 41 PRIVATE 30 42 31 !! * Accessibility32 PUBLIC dyn_nxt ! routine called by step.F90 43 PUBLIC dyn_nxt ! routine called by step.F90 44 33 45 !! * Substitutions 34 46 # include "domzgr_substitute.h90" 35 !!---------------------------------------------------------------------- 47 !!------------------------------------------------------------------------- 48 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 49 !! $Id$ 50 !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 51 !!------------------------------------------------------------------------- 36 52 37 53 CONTAINS … … 59 75 !! - Update un,vn arrays, the now horizontal velocity 60 76 !! 61 !! History :62 !! ! 87-02 (P. Andrich, D. L Hostis) Original code63 !! ! 90-10 (C. Levy, G. Madec)64 !! ! 93-03 (M. Guyon) symetrical conditions65 !! ! 97-02 (G. Madec & M. Imbard) opa, release 8.066 !! ! 97-04 (A. Weaver) Euler forward step67 !! ! 97-06 (G. Madec) lateral boudary cond., lbc routine68 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module69 !! ! 02-10 (C. Talandier, A-M. Treguier) Open boundary cond.70 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization71 !! " ! 07-07 (D. Storkey) Calls to BDY routines.72 77 !!---------------------------------------------------------------------- 73 !! * Arguments74 78 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 76 !! * Local declarations 79 !! 77 80 INTEGER :: ji, jj, jk ! dummy loop indices 78 81 REAL(wp) :: z2dt ! temporary scalar 79 REAL(wp) :: zsshun1, zsshvn1 ! temporary scalar 80 !! Variable volume 81 REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace 82 zsshub, zsshun, zsshua, & 83 zsshvb, zsshvn, zsshva 84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 85 zfse3ub, zfse3un, zfse3ua, & ! 3D workspace 86 zfse3vb, zfse3vn, zfse3va 87 !!---------------------------------------------------------------------- 88 !! OPA 9.0 , LOCEAN-IPSL (2005) 89 !! $Id$ 90 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 91 !!---------------------------------------------------------------------- 82 REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar 83 REAL(wp) :: zve3a , zve3n , zve3b ! - - 84 REAL(wp) :: ze3u_b, ze3u_n, ze3u_a ! - - 85 REAL(wp) :: ze3v_b, ze3v_n, ze3v_a ! - - 86 REAL(wp) :: zuf , zvf ! - - 87 !!---------------------------------------------------------------------- 92 88 93 89 IF( kt == nit000 ) THEN … … 102 98 103 99 !! Explicit physics with thickness weighted updates 104 IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN 105 106 ! Sea surface elevation time stepping 107 ! ----------------------------------- 108 ! 109 DO jj = 1, jpjm1 110 DO ji = 1,jpim1 111 112 ! Sea Surface Height at u-point before 113 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 114 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) & 115 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * sshbb(ji+1,jj ) ) 116 117 ! Sea Surface Height at v-point before 118 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 119 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) & 120 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * sshbb(ji ,jj+1) ) 121 122 ! Sea Surface Height at u-point after 123 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 124 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) & 125 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * ssha(ji+1,jj ) ) 126 127 ! Sea Surface Height at v-point after 128 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 129 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) & 130 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * ssha(ji ,jj+1) ) 131 132 END DO 133 END DO 134 135 ! Boundaries conditions 136 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. ) 137 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. ) 138 139 ! Scale factors at before and after time step 140 ! ------------------------------------------- 141 CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua ) 142 CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va ) 143 144 ! Asselin filtered scale factor at now time step 145 ! ---------------------------------------------- 146 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 147 CALL dom_vvl_sf_ini( 'U', zfse3un ) ; CALL dom_vvl_sf_ini( 'V', zfse3vn ) 148 ELSE 149 zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:) 150 zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:) 151 CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ; CALL dom_vvl_sf( zsshvn, 'V', zfse3vn ) 152 ENDIF 153 154 ! Thickness weighting 155 ! ------------------- 100 101 ! Lateral boundary conditions on ( ua, va ) 102 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) 103 104 ! Next velocity 105 ! ------------- 106 #if defined key_dynspg_flt 107 ! Leap-frog time stepping already done in dynspg_flt.F routine 108 #else 109 IF( lk_vvl ) THEN ! Varying levels 156 110 DO jk = 1, jpkm1 157 DO jj = 1, jpj 158 DO ji = 1, jpi 159 ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk) 160 va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk) 161 162 un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk) 163 vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk) 164 165 ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk) 166 vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk) 167 END DO 168 END DO 169 END DO 170 171 ENDIF 172 173 ! Lateral boundary conditions on ( ua, va ) 174 CALL lbc_lnk( ua, 'U', -1. ) 175 CALL lbc_lnk( va, 'V', -1. ) 176 177 ! ! =============== 178 DO jk = 1, jpkm1 ! Horizontal slab 179 ! ! =============== 180 ! Next velocity 181 ! ------------- 182 #if defined key_dynspg_flt 183 ! Leap-frog time stepping already done in dynspg.F routine 184 #else 185 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 186 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 187 ! Leap-frog time stepping 188 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 189 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 190 END DO 191 END DO 192 193 IF( lk_vvl ) THEN 194 ! Unweight velocities prior to updating open boundaries. 195 196 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 197 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 198 ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk) 199 va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk) 200 201 un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk) 202 vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk) 203 204 ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk) 205 vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk) 206 END DO 207 END DO 208 209 ENDIF 111 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 112 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 113 & / fse3u_a(:,:,jk) * umask(:,:,jk) 114 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 115 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 116 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 117 END DO 118 ELSE 119 DO jk = 1, jpkm1 120 ! Leap-frog time stepping 121 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 122 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 123 END DO 124 ENDIF 210 125 211 126 # if defined key_obc 212 ! ! ===============213 END DO ! End of slab214 ! ! ===============215 127 ! Update (ua,va) along open boundaries (only in the rigid-lid case) 216 128 CALL obc_dyn( kt ) 217 129 218 130 IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 219 ! Flather boundary condition :131 ! Flather boundary condition : 220 132 ! - Update sea surface height on each open boundary 221 133 ! sshn (= after ssh) for explicit case … … 223 135 ! - Correct the barotropic velocities 224 136 CALL obc_dyn_bt( kt ) 225 226 ! Boundary conditions on sshn ( after ssh)137 ! 138 ! Boundary conditions on sshn ( after ssh) 227 139 CALL lbc_lnk( sshn, 'T', 1. ) 228 229 IF(ln_ctl) THEN ! print sum trends (used for debugging) 230 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask) 231 ENDIF 232 233 IF ( ln_vol_cst ) CALL obc_vol( kt ) 234 235 ENDIF 236 237 ! ! =============== 238 DO jk = 1, jpkm1 ! Horizontal slab 239 ! ! =============== 140 ! 141 IF( ln_vol_cst ) CALL obc_vol( kt ) 142 ! 143 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) 144 ENDIF 145 240 146 # elif defined key_bdy 241 ! ! ===============242 END DO ! End of slab243 ! ! ===============244 147 ! Update (ua,va) along open boundaries (for exp or ts options). 245 IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN246 148 IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 149 ! 247 150 CALL bdy_dyn_frs( kt ) 248 249 IF ( ln_bdy_fla ) THEN 250 251 ua_e(:,:)=0.0 252 va_e(:,:)=0.0 253 151 ! 152 IF( ln_bdy_fla ) THEN 153 ua_e(:,:) = 0.e0 154 va_e(:,:) = 0.e0 254 155 ! Set these variables for use in bdy_dyn_fla 255 156 hu_e(:,:) = hu(:,:) 256 157 hv_e(:,:) = hv(:,:) 257 258 DO jk = 1, jpkm1 259 !! Vertically integrated momentum trends 158 DO jk = 1, jpkm1 !! Vertically integrated momentum trends 260 159 ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 261 160 va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 262 161 END DO 263 264 162 DO jk = 1 , jpkm1 265 163 ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:) 266 164 va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:) 267 165 END DO 268 269 166 CALL bdy_dta_bt( kt+1, 0) 270 167 CALL bdy_dyn_fla 271 272 168 ENDIF 273 169 ! 274 170 DO jk = 1 , jpkm1 275 171 ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:) 276 172 va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:) 277 173 END DO 278 279 ENDIF 280 281 ! ! =============== 282 DO jk = 1, jpkm1 ! Horizontal slab 283 ! ! =============== 174 ENDIF 284 175 # endif 176 285 177 # if defined key_agrif 286 ! ! ===============287 END DO ! End of slab288 ! ! ===============289 ! Update (ua,va) along open boundaries (only in the rigid-lid case)290 178 CALL Agrif_dyn( kt ) 291 ! ! ===============292 DO jk = 1, jpkm1 ! Horizontal slab293 ! ! ===============294 179 # endif 295 180 #endif 296 181 297 ! Time filter and swap of dynamics arrays 298 ! ------------------------------------------ 299 IF( neuler == 0 .AND. kt == nit000 ) THEN 300 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 301 ! caution: don't use (:,:) for this loop 302 ! it causes optimization problems on NEC in auto-tasking 182 ! Time filter and swap of dynamics arrays 183 ! ------------------------------------------ 184 IF( neuler == 0 .AND. kt == nit000 ) THEN 185 DO jk = 1, jpkm1 186 un(:,:,jk) = ua(:,:,jk) 187 vn(:,:,jk) = va(:,:,jk) 188 END DO 189 ELSE 190 IF( lk_vvl ) THEN ! Varying levels 191 DO jk = 1, jpkm1 ! filter applied on thickness weighted velocities 303 192 DO jj = 1, jpj 304 193 DO ji = 1, jpi 305 zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk) 306 zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk) 307 ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 308 vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 309 zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 310 zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 311 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 312 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 194 ze3u_a = fse3u_a(ji,jj,jk) 195 ze3v_a = fse3v_a(ji,jj,jk) 196 ze3u_n = fse3u_n(ji,jj,jk) 197 ze3v_n = fse3v_n(ji,jj,jk) 198 ze3u_b = fse3u_b(ji,jj,jk) 199 ze3v_b = fse3v_b(ji,jj,jk) 200 ! 201 zue3a = ua(ji,jj,jk) * ze3u_a 202 zve3a = va(ji,jj,jk) * ze3v_a 203 zue3n = un(ji,jj,jk) * ze3u_n 204 zve3n = vn(ji,jj,jk) * ze3v_n 205 zue3b = ub(ji,jj,jk) * ze3u_b 206 zve3b = vb(ji,jj,jk) * ze3v_b 207 ! 208 ub(ji,jj,jk) = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 209 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 210 vb(ji,jj,jk) = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 211 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 212 un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk) 213 vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk) 313 214 END DO 314 215 END DO 315 ELSE ! Fixed levels 216 END DO 217 ELSE ! Fixed levels 218 DO jk = 1, jpkm1 ! filter applied on velocities 316 219 DO jj = 1, jpj 317 220 DO ji = 1, jpi 318 ! Euler (forward) time stepping 319 ub(ji,jj,jk) = un(ji,jj,jk) 320 vb(ji,jj,jk) = vn(ji,jj,jk) 221 zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 222 zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 223 ub(ji,jj,jk) = zuf 224 vb(ji,jj,jk) = zvf 321 225 un(ji,jj,jk) = ua(ji,jj,jk) 322 226 vn(ji,jj,jk) = va(ji,jj,jk) 323 227 END DO 324 228 END DO 325 ENDIF 326 ELSE 327 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 328 ! caution: don't use (:,:) for this loop 329 ! it causes optimization problems on NEC in auto-tasking 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk) 333 zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk) 334 ub(ji,jj,jk) = ( atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) & 335 & + atfp1 * un(ji,jj,jk) ) * zsshun1 336 vb(ji,jj,jk) = ( atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) & 337 & + atfp1 * vn(ji,jj,jk) ) * zsshvn1 338 zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 339 zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 340 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 341 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 342 END DO 343 END DO 344 ELSE ! Fixed levels 345 DO jj = 1, jpj 346 DO ji = 1, jpi 347 ! Leap-frog time stepping 348 ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 349 vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 350 un(ji,jj,jk) = ua(ji,jj,jk) 351 vn(ji,jj,jk) = va(ji,jj,jk) 352 END DO 353 END DO 354 ENDIF 229 END DO 355 230 ENDIF 356 ! ! ===============357 END DO ! End of slab358 ! ! ===============359 360 IF(ln_ctl) THEN361 CALL prt_ctl(tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &362 & tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)363 231 ENDIF 364 232 … … 367 235 #endif 368 236 237 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 238 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 239 ! 369 240 END SUBROUTINE dyn_nxt 370 241 -
trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r1146 r1438 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 9.0 ! 2005-11 (V. Garnier, G. Madec, L. Bessieres) Original code 7 !!---------------------------------------------------------------------- 6 8 #if defined key_dynspg_exp || defined key_esopa 7 9 !!---------------------------------------------------------------------- … … 11 13 !! pressure gradient in the free surface constant 12 14 !! volume case with vector optimization 13 !! exp_rst : read/write the explicit restart fields in the ocean restart file14 15 !!---------------------------------------------------------------------- 15 !! * Modules used16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain … … 31 31 PRIVATE 32 32 33 !! * Accessibility 34 PUBLIC dyn_spg_exp ! routine called by step.F90 35 PUBLIC exp_rst ! routine called by istate.F90 33 PUBLIC dyn_spg_exp ! routine called by step.F90 36 34 37 35 !! * Substitutions … … 39 37 # include "vectopt_loop_substitute.h90" 40 38 !!---------------------------------------------------------------------- 41 !! OPA 9.0 , LOCEAN-IPSL (2005)39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 42 40 !! $Id$ 43 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 42 !!---------------------------------------------------------------------- 45 46 43 47 44 CONTAINS … … 62 59 !! -1- Compute the now surface pressure gradient 63 60 !! -2- Add it to the general trend 64 !! -3- Compute the horizontal divergence of velocities65 !! - the now divergence is given by :66 !! zhdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )67 !! - integrate the horizontal divergence from the bottom68 !! to the surface69 !! - apply lateral boundary conditions on zhdivn70 !! -4- Estimate the after sea surface elevation from the kinematic71 !! surface boundary condition:72 !! zssha = sshb - 2 rdt ( zhdiv + emp )73 !! - Time filter applied on now sea surface elevation to avoid74 !! the divergence of two consecutive time-steps and swap of free75 !! surface arrays to start the next time step:76 !! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]77 !! sshn = zssha78 !! - apply lateral boundary conditions on sshn79 61 !! 80 62 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 63 !!--------------------------------------------------------------------- 64 INTEGER, INTENT( in ) :: kt ! ocean time-step index 81 65 !! 82 !! References :83 !!84 !! History :85 !! 9.0 ! 05-11 (V. Garnier, G. Madec, L. Bessieres) Original code86 !!87 !!---------------------------------------------------------------------88 !! * Arguments89 INTEGER, INTENT( in ) :: kt ! ocean time-step index90 91 !! * Local declarations92 66 INTEGER :: ji, jj, jk ! dummy loop indices 93 REAL(wp) :: z2dt, zraur, zssha ! temporary scalars94 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary arrays95 & zhdiv96 67 !!---------------------------------------------------------------------- 97 68 … … 104 75 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 105 76 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 106 107 CALL exp_rst( nit000, 'READ' ) ! read or initialize the following fields:108 ! ! sshb, sshn109 110 77 ENDIF 111 78 112 IF( .NOT. lk_vvl ) THEN 79 ! read or estimate sea surface height and vertically integrated velocities 80 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 113 81 114 ! 0. Initialization115 ! -----------------116 ! read or estimate sea surface height and vertically integrated velocities117 IF( lk_obc ) CALL obc_dta_bt( kt, 0)118 z2dt = 2. * rdt ! time step: leap-frog119 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest120 zraur = 1.e0 / rauw82 ! Surface pressure gradient (now) 83 DO jj = 2, jpjm1 84 DO ji = fs_2, fs_jpim1 ! vector opt. 85 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 86 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 87 END DO 88 END DO 121 89 122 ! 1. Surface pressure gradient (now)123 ! ----------------------------90 ! Add the surface pressure trend to the general trend 91 DO jk = 1, jpkm1 124 92 DO jj = 2, jpjm1 125 93 DO ji = fs_2, fs_jpim1 ! vector opt. 126 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)127 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)94 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 95 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 128 96 END DO 129 97 END DO 98 END DO 130 99 131 ! 2. Add the surface pressure trend to the general trend132 ! ------------------------------------------------------133 DO jk = 1, jpkm1134 DO jj = 2, jpjm1135 DO ji = fs_2, fs_jpim1 ! vector opt.136 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)137 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)138 END DO139 END DO140 END DO141 142 ! 3. Vertical integration of horizontal divergence of velocities143 ! --------------------------------144 zhdiv(:,:) = 0.e0145 DO jk = jpkm1, 1, -1146 DO jj = 2, jpjm1147 DO ji = fs_2, fs_jpim1 ! vector opt.148 zhdiv(ji,jj) = zhdiv(ji,jj) + ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * un(ji ,jj ,jk) &149 & - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) &150 & + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * vn(ji ,jj ,jk) &151 & - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) &152 & / ( e1t(ji,jj) * e2t(ji,jj) )153 END DO154 END DO155 END DO156 157 #if defined key_obc158 ! open boundaries (div must be zero behind the open boundary)159 ! mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column160 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east161 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west162 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north163 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south164 #endif165 166 ! 4. Sea surface elevation time stepping167 ! --------------------------------------168 ! Euler (forward) time stepping, no time filter169 IF( neuler == 0 .AND. kt == nit000 ) THEN170 DO jj = 1, jpj171 DO ji = 1, jpi172 ! after free surface elevation173 zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)174 ! swap of arrays175 sshb(ji,jj) = sshn(ji,jj)176 sshn(ji,jj) = zssha177 END DO178 END DO179 ELSE180 ! Leap-frog time stepping and time filter181 DO jj = 1, jpj182 DO ji = 1, jpi183 ! after free surface elevation184 zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)185 ! time filter and array swap186 sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)187 sshn(ji,jj) = zssha188 END DO189 END DO190 ENDIF191 192 ELSE !! Variable volume, ssh time-stepping already done193 194 ! read or estimate sea surface height and vertically integrated velocities195 IF( lk_obc ) CALL obc_dta_bt( kt, 0 )196 197 ! Sea surface elevation swap198 ! -----------------------------199 !200 sshbb(:,:) = sshb(:,:) ! Needed for the dynamics time-stepping201 !202 IF( neuler == 0 .AND. kt == nit000 ) THEN203 ! No time filter204 sshb(:,:) = sshn(:,:)205 sshn(:,:) = ssha(:,:)206 ELSE207 ! Time filter208 sshb(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:)209 sshn(:,:) = ssha(:,:)210 ENDIF211 212 ENDIF213 214 ! Boundary conditions on sshn215 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )216 217 ! write filtered free surface arrays in restart file218 ! --------------------------------------------------219 IF( lrst_oce ) CALL exp_rst( kt, 'WRITE' )220 221 IF(ln_ctl) THEN ! print sum trends (used for debugging)222 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask)223 ENDIF224 225 100 END SUBROUTINE dyn_spg_exp 226 101 227 SUBROUTINE exp_rst( kt, cdrw )228 !!---------------------------------------------------------------------229 !! *** ROUTINE exp_rst ***230 !!231 !! ** Purpose : Read or write explicit arrays in restart file232 !!----------------------------------------------------------------------233 INTEGER , INTENT(in) :: kt ! ocean time-step234 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag235 !236 !!----------------------------------------------------------------------237 !238 IF( TRIM(cdrw) == 'READ' ) THEN239 IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN240 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) )241 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) )242 IF( neuler == 0 ) sshb(:,:) = sshn(:,:)243 ELSE244 IF( nn_rstssh == 1 ) THEN245 sshb(:,:) = 0.e0246 sshn(:,:) = 0.e0247 ENDIF248 ENDIF249 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN250 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) )251 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) )252 ENDIF253 !254 END SUBROUTINE exp_rst255 102 #else 256 103 !!---------------------------------------------------------------------- … … 261 108 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt 262 109 END SUBROUTINE dyn_spg_exp 263 SUBROUTINE exp_rst( kt, cdrw ) ! Empty routine264 INTEGER , INTENT(in) :: kt ! ocean time-step265 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag266 WRITE(*,*) 'exp_rst: You should not have seen this print! error?', kt, cdrw267 END SUBROUTINE exp_rst268 110 #endif 269 111 270 112 !!====================================================================== 271 113 END MODULE dynspg_exp -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r1200 r1438 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History 8.0 ! 98-05 (G. Roullet) free surface 7 !! ! 98-10 (G. Madec, M. Imbard) release 8.2 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 9 !! " " ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 10 !! 9.0 ! 04-08 (C. Talandier) New trends organization 11 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 12 !! " " ! 06-07 (S. Masson) distributed restart using iom 13 !! " " ! 05-01 (J.Chanut, A.Sellar) Calls to BDY routines. 6 !! History OPA ! 1998-05 (G. Roullet) free surface 7 !! ! 1998-10 (G. Madec, M. Imbard) release 8.2 8 !! NEMO O.1 ! 2002-08 (G. Madec) F90: Free form and module 9 !! - ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries 10 !! 1.0 ! 2004-08 (C. Talandier) New trends organization 11 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 12 !! 2.0 ! 2006-07 (S. Masson) distributed restart using iom 13 !! - ! 2006-08 (J.Chanut, A.Sellar) Calls to BDY routines. 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 14 15 !!---------------------------------------------------------------------- 15 16 #if defined key_dynspg_flt || defined key_esopa 16 17 !!---------------------------------------------------------------------- 17 18 !! 'key_dynspg_flt' filtered free surface 18 !!----------------------------------------------------------------------19 19 !!---------------------------------------------------------------------- 20 20 !! dyn_spg_flt : update the momentum trend with the surface pressure … … 53 53 PRIVATE 54 54 55 PUBLIC dyn_spg_flt ! routine called by step.F9056 PUBLIC flt_rst ! routine called by istate.F9055 PUBLIC dyn_spg_flt ! routine called by step.F90 56 PUBLIC flt_rst ! routine called by istate.F90 57 57 58 58 !! * Substitutions … … 60 60 # include "vectopt_loop_substitute.h90" 61 61 !!---------------------------------------------------------------------- 62 !! OPA 9.0 , LOCEAN-IPSL (2005)62 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 63 63 !! $Id$ 64 64 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 80 80 !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + rnu btda ) 81 81 !! where sshn is the free surface elevation and btda is the after 82 !! of the free surface elevation 83 !! -1- compute the after sea surface elevation from the kinematic 84 !! surface boundary condition: 85 !! zssha = sshb + 2 rdt ( wn - emp ) 86 !! Time filter applied on now sea surface elevation to avoid 87 !! the divergence of two consecutive time-steps and swap of free 88 !! surface arrays to start the next time step: 89 !! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ] 90 !! sshn = zssha 91 !! -2- evaluate the surface presure trend (including the addi- 82 !! time derivative of the free surface elevation 83 !! -1- evaluate the surface presure trend (including the addi- 92 84 !! tional force) in three steps: 93 85 !! a- compute the right hand side of the elliptic equation: … … 105 97 !! iterative solver 106 98 !! c- apply the solver by a call to sol... routine 107 !! - 3- compute and add the free surface pressure gradient inclu-99 !! -2- compute and add the free surface pressure gradient inclu- 108 100 !! ding the additional force used to stabilize the equation. 109 101 !! … … 112 104 !! References : Roullet and Madec 1999, JGR. 113 105 !!--------------------------------------------------------------------- 114 !! * Modules used 115 USE oce , ONLY : zub => ta, & ! ta used as workspace 116 zvb => sa ! sa " " 117 118 INTEGER, INTENT( in ) :: kt ! ocean time-step index 119 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) 106 USE oce, ONLY : zub => ta ! ta used as workspace 107 USE oce, ONLY : zvb => sa ! ta used as workspace 108 !! 109 INTEGER, INTENT( in ) :: kt ! ocean time-step index 110 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) 120 111 !! 121 112 INTEGER :: ji, jj, jk ! dummy loop indices 122 REAL(wp) :: z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 123 & znurau, zgcb, zbtd, & ! " " 124 & ztdgu, ztdgv ! " " 125 !! Variable volume 126 REAL(wp), DIMENSION(jpi,jpj) :: & 127 zsshub, zsshua, zsshvb, zsshva, zssha ! 2D workspace 128 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace 129 zfse3ub, zfse3ua, zfse3vb, zfse3va ! 3D workspace 113 REAL(wp) :: z2dt, z2dtg, zraur, znugdt ! temporary scalars 114 REAL(wp) :: znurau, zgcb, zbtd ! " " 115 REAL(wp) :: ztdgu, ztdgv ! " " 130 116 !!---------------------------------------------------------------------- 131 117 ! … … 143 129 ! when using agrif, sshn, gcx have to be read in istate 144 130 IF (.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 145 ! ! gcx, gcxb , sshb, sshn131 ! ! gcx, gcxb 146 132 ENDIF 147 133 … … 158 144 IF( lk_vvl ) THEN ! variable volume 159 145 160 DO jj = 1, jpjm1161 DO ji = 1,jpim1162 163 ! Sea Surface Height at u-point before164 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &165 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &166 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )167 168 ! Sea Surface Height at v-point before169 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &170 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &171 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )172 173 ! Sea Surface Height at u-point after174 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &175 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &176 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )177 178 ! Sea Surface Height at v-point after179 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &180 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &181 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )182 183 END DO184 END DO185 186 ! Boundaries conditions187 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. )188 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. )189 190 ! Scale factors at before and after time step191 ! -------------------------------------------192 CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua )193 CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va )194 195 196 ! Thickness weighting197 ! -------------------198 DO jk = 1, jpkm1199 DO jj = 1, jpj200 DO ji = 1, jpi201 ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)202 va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)203 204 zub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)205 zvb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)206 END DO207 END DO208 END DO209 210 146 ! Evaluate the masked next velocity (effect of the additional force not included) 147 ! ------------------- (thickness weighted velocity, surface pressure gradient already included in dyn_hpg) 211 148 DO jk = 1, jpkm1 212 149 DO jj = 2, jpjm1 213 150 DO ji = fs_2, fs_jpim1 ! vector opt. 214 ua(ji,jj,jk) = ( zub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) /zfse3ua(ji,jj,jk) * umask(ji,jj,jk) 215 va(ji,jj,jk) = ( zvb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) /zfse3va(ji,jj,jk) * vmask(ji,jj,jk) 151 ua(ji,jj,jk) = ( ub(ji,jj,jk) * fse3u_b(ji,jj,jk) & 152 & + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk) ) & 153 & / fse3u_a(ji,jj,jk) * umask(ji,jj,jk) 154 va(ji,jj,jk) = ( vb(ji,jj,jk) * fse3v_b(ji,jj,jk) & 155 & + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk) ) & 156 & / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk) 216 157 END DO 217 158 END DO … … 227 168 END DO 228 169 END DO 229 230 ! Add the surface pressure trend to the general trend 170 ! 171 ! add the surface pressure trend to the general trend and 172 ! evaluate the masked next velocity (effect of the additional force not included) 231 173 DO jk = 1, jpkm1 232 174 DO jj = 2, jpjm1 233 175 DO ji = fs_2, fs_jpim1 ! vector opt. 234 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)235 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)176 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) ) ) * umask(ji,jj,jk) 177 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) ) ) * vmask(ji,jj,jk) 236 178 END DO 237 179 END DO 238 180 END DO 239 240 ! Evaluate the masked next velocity (effect of the additional force not included) 241 DO jk = 1, jpkm1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 245 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 181 ! 250 182 ENDIF 251 183 … … 311 243 CALL lbc_lnk( spgv, 'V', -1. ) 312 244 313 IF( lk_vvl ) CALL sol_mat( kt ) 245 IF( lk_vvl ) CALL sol_mat( kt ) ! build the matrix at kt (vvl case only) 314 246 315 247 ! Right hand side of the elliptic equation and first guess … … 347 279 ! Relative precision (computation on one processor) 348 280 ! ------------------ 349 rnorme =0. 281 rnorme =0.e0 350 282 rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) ) 351 283 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain … … 413 345 ENDIF 414 346 #endif 415 ! 7.Add the trends multiplied by z2dt to the after velocity416 ! ------------------------------------------------------- ----347 ! Add the trends multiplied by z2dt to the after velocity 348 ! ------------------------------------------------------- 417 349 ! ( c a u t i o n : (ua,va) here are the after velocity not the 418 350 ! trend, the leap-frog time stepping will not … … 421 353 DO jj = 2, jpjm1 422 354 DO ji = fs_2, fs_jpim1 ! vector opt. 423 ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)424 va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)355 ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk) 356 va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk) 425 357 END DO 426 358 END DO 427 359 END DO 428 429 ! Sea surface elevation time stepping430 ! -----------------------------------431 IF( lk_vvl ) THEN ! after free surface elevation432 zssha(:,:) = ssha(:,:)433 ELSE434 zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1)435 ENDIF436 #if defined key_obc437 # if defined key_agrif438 IF ( Agrif_Root() ) THEN439 # endif440 zssha(:,:)=zssha(:,:)*obctmsk(:,:)441 CALL lbc_lnk(zssha,'T',1.) ! absolutly compulsory !! (jmm)442 # if defined key_agrif443 ENDIF444 # endif445 #endif446 447 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler (forward) time stepping, no time filter448 ! swap of arrays449 sshb(:,:) = sshn (:,:)450 sshn(:,:) = zssha(:,:)451 ELSE ! Leap-frog time stepping and time filter452 ! time filter and array swap453 sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:)454 sshn(:,:) = zssha(:,:)455 ENDIF456 360 457 361 ! write filtered free surface arrays in restart file 458 362 ! -------------------------------------------------- 459 363 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 460 461 ! print sum trends (used for debugging)462 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg - ssh: ', mask1=tmask )463 364 ! 464 365 END SUBROUTINE dyn_spg_flt … … 481 382 CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) ) 482 383 CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) ) 483 CALL iom_get( numror, jpdom_autoglo, 'sshb', sshb(:,:) ) 484 CALL iom_get( numror, jpdom_autoglo, 'sshn', sshn(:,:) ) 485 IF( neuler == 0 ) THEN 486 sshb(:,:) = sshn(:,:) 487 gcxb(:,:) = gcx (:,:) 488 ENDIF 384 IF( neuler == 0 ) gcxb(:,:) = gcx (:,:) 489 385 ELSE 490 386 gcx (:,:) = 0.e0 491 387 gcxb(:,:) = 0.e0 492 IF( nn_rstssh == 1 ) THEN493 sshb(:,:) = 0.e0494 sshn(:,:) = 0.e0495 ENDIF496 388 ENDIF 497 389 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 498 390 ! Caution : extra-hallow 499 391 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 500 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )392 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) ) 501 393 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 502 CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:) )503 CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:) )504 394 ENDIF 505 395 ! -
trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
r1152 r1438 43 43 !! ------------------------ 44 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop 45 sshn_b, sshb_b, & ! sea surface heigth (now, before)46 u n_b , vn_b ! vertically integrated horizontal velocities (now)45 un_e , vn_e, & ! vertically integrated horizontal velocities (now) 46 ua_e , va_e ! vertically integrated horizontal velocities (after) 47 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop 48 sshb_e, sshn_e, ssha_e, & ! sea surface heigth (before,now,after) 49 ua_e , va_e, & ! vertically integrated horizontal velocities (after) 50 hu_e , hv_e ! depth arrays for the barotropic solution 48 sshn_e, ssha_e, & ! sea surface heigth (now, after) 49 hu_e , hv_e ! depth arrays for the barotropic solution 51 50 #endif 52 51 -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r1241 r1438 1 1 MODULE dynspg_ts 2 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend 5 !!====================================================================== 6 !! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 7 !! " " ! 05-11 (V. Garnier, G. Madec) optimization 8 !! 9.0 ! 06-08 (S. Masson) distributed restart using iom 9 !! " ! 08-01 (R. Benshila) change averaging method 10 !! " ! 07-07 (D. Storkey) calls to BDY routines 11 !!--------------------------------------------------------------------- 3 !! History : 1.0 ! 04-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 05-11 (V. Garnier, G. Madec) optimization 5 !! - ! 06-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 07-07 (D. Storkey) calls to BDY routines 7 !! - ! 08-01 (R. Benshila) change averaging method 8 !!--------------------------------------------------------------------- 12 9 #if defined key_dynspg_ts || defined key_esopa 13 10 !!---------------------------------------------------------------------- … … 19 16 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file 20 17 !!---------------------------------------------------------------------- 21 !! * Modules used22 18 USE oce ! ocean dynamics and tracers 23 19 USE dom_oce ! ocean space and time domain … … 49 45 PUBLIC ts_rst ! routine called by istate.F90 50 46 51 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter 52 & ftsw, ftse ! (only used with een vorticity scheme) 53 47 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter 48 REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme) 54 49 55 50 !! * Substitutions 56 51 # include "domzgr_substitute.h90" 57 52 # include "vectopt_loop_substitute.h90" 58 !!---------------------------------------------------------------------- 59 !! OPA 9.0 , LOCEAN-IPSL (2005)53 !!------------------------------------------------------------------------- 54 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 60 55 !! $Id$ 61 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt62 !!---------------------------------------------------------------------- 56 !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 57 !!------------------------------------------------------------------------- 63 58 64 59 CONTAINS … … 85 80 !! (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b 86 81 !! and zvX_b (X specifying after, now or before). 87 !! -3- Update of sea surface height from time averaged barotropic 88 !! variables. 89 !! - apply lateral boundary conditions on sshn. 90 !! -4- The new general trend becomes : 91 !! ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H 82 !! -3- The new general trend becomes : 83 !! ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H 92 84 !! 93 85 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend … … 96 88 !!--------------------------------------------------------------------- 97 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 98 99 !! * Local declarations 90 !! 100 91 INTEGER :: ji, jj, jk, jit ! dummy loop indices 101 92 INTEGER :: icycle ! temporary scalar … … 107 98 zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays 108 99 zua, zva, zub, zvb, & ! " " 109 z ssha_b, zua_b, zva_b, & ! " "110 zub_e, zvb_e, 111 zu n_e, zvn_e ! " "100 zua_e, zva_e, zssha_e, & ! " " 101 zub_e, zvb_e, zsshb_e, & ! " " 102 zu_sum, zv_sum, zssh_sum 112 103 !! Variable volume 113 104 REAL(wp), DIMENSION(jpi,jpj) :: & 114 105 zspgu_1, zspgv_1, zsshun_e, zsshvn_e ! 2D workspace 115 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace116 106 !!---------------------------------------------------------------------- 117 107 118 108 ! Arrays initialization 119 109 ! --------------------- 120 zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0121 zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0122 110 zhdiv(:,:) = 0.e0 123 111 … … 133 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 134 122 135 ssha_e(:,:) = sshn(:,:) 136 ua_e(:,:) = un_b(:,:) 137 va_e(:,:) = vn_b(:,:) 138 hu_e(:,:) = hu(:,:) 139 hv_e(:,:) = hv(:,:) 140 123 zssha_e(:,:) = sshn(:,:) 124 zua_e (:,:) = un_e(:,:) 125 zva_e (:,:) = vn_e(:,:) 126 hu_e (:,:) = hu (:,:) 127 hv_e (:,:) = hv (:,:) 141 128 IF( ln_dynvor_een ) THEN 142 129 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 … … 170 157 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 171 158 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 172 END DO 159 END DO 173 160 END DO 174 161 … … 193 180 zfact2 = 0.5 * 0.5 194 181 zraur = 1. / rauw ! 1 / volumic mass of pure water 195 182 196 183 ! ----------------------------------------------------------------------------- 197 184 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 215 202 zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 216 203 ! ! Vertically integrated transports (before) 217 zub(ji,1) = zub(ji,1) + fse3u (ji,1,jk) * ub(ji,1,jk)218 zvb(ji,1) = zvb(ji,1) + fse3v (ji,1,jk) * vb(ji,1,jk)204 zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 205 zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 219 206 ! ! Planetary vorticity transport fluxes (now) 220 207 zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) … … 228 215 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 229 216 ! ! Vertically integrated transports (before) 230 zub(:,:) = zub(:,:) + fse3u (:,:,jk) * ub(:,:,jk)231 zvb(:,:) = zvb(:,:) + fse3v (:,:,jk) * vb(:,:,jk)217 zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 218 zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 232 219 ! ! Planetary vorticity (now) 233 220 zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) … … 304 291 !---------------- 305 292 ! Number of iteration of the barotropic loop 306 icycle = 3 * nn_baro / 2293 icycle = 2 * nn_baro + 1 307 294 308 295 ! variables for the barotropic equations 309 sshb_e (:,:) = sshn_b(:,:) ! (barotropic) sea surface height (before and now) 310 sshn_e (:,:) = sshn_b(:,:) 311 zub_e (:,:) = un_b (:,:) ! barotropic transports issued from the barotropic equations (before and now) 312 zvb_e (:,:) = vn_b (:,:) 313 zun_e (:,:) = un_b (:,:) 314 zvn_e (:,:) = vn_b (:,:) 315 zssha_b(:,:) = 0.e0 316 zua_b (:,:) = 0.e0 317 zva_b (:,:) = 0.e0 318 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 319 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 320 IF( lk_vvl ) THEN 321 zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point 322 zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point 323 zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point 324 zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point 296 zu_sum (:,:) = 0.e0 297 zv_sum (:,:) = 0.e0 298 zssh_sum(:,:) = 0.e0 299 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 300 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 301 zsshb_e (:,:) = sshn_e(:,:) ! (barotropic) sea surface height (before and now) 302 ! vertical sum 303 un_e (:,:) = 0.e0 304 vn_e (:,:) = 0.e0 305 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 306 DO jk = 1, jpkm1 307 DO ji = 1, jpij 308 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 309 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 310 END DO 311 END DO 312 ELSE ! No vector opt. 313 DO jk = 1, jpkm1 314 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 315 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 316 END DO 325 317 ENDIF 318 zub_e (:,:) = un_e(:,:) 319 zvb_e (:,:) = vn_e(:,:) 320 326 321 327 322 ! set ssh corrections to 0 … … 352 347 DO jj = 2, jpjm1 353 348 DO ji = fs_2, fs_jpim1 ! vector opt. 354 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun_e(ji ,jj) &355 & -e2u(ji-1,jj ) * zun_e(ji-1,jj) &356 & +e1v(ji ,jj ) * zvn_e(ji ,jj) &357 & -e1v(ji ,jj-1) * zvn_e(ji ,jj-1) ) &349 zhdiv(ji,jj) = ( e2u(ji ,jj ) * un_e(ji ,jj) & 350 & -e2u(ji-1,jj ) * un_e(ji-1,jj) & 351 & +e1v(ji ,jj ) * vn_e(ji ,jj) & 352 & -e1v(ji ,jj-1) * vn_e(ji ,jj-1) ) & 358 353 & / (e1t(ji,jj)*e2t(ji,jj)) 359 354 END DO … … 370 365 371 366 #if defined key_bdy 372 373 374 375 376 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 370 END DO 371 END DO 377 372 #endif 378 373 … … 381 376 DO jj = 2, jpjm1 382 377 DO ji = fs_2, fs_jpim1 ! vector opt. 383 ssha_e(ji,jj) = (sshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) &384 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)378 zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & 379 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 385 380 END DO 386 381 END DO … … 388 383 ! evolution of the barotropic transport ( following the vorticity scheme used) 389 384 ! ---------------------------------------------------------------------------- 390 zwx(:,:) = e2u(:,:) * zun_e(:,:)391 zwy(:,:) = e1v(:,:) * zvn_e(:,:)385 zwx(:,:) = e2u(:,:) * un_e(:,:) 386 zwy(:,:) = e1v(:,:) * vn_e(:,:) 392 387 393 388 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 396 391 ! surface pressure gradient 397 392 IF( lk_vvl) THEN 398 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &399 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)400 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &401 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)393 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 394 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 395 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 402 397 ELSE 403 398 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 412 407 zcvbt =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 413 408 ! after transports 414 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)415 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)409 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 410 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 416 411 END DO 417 412 END DO … … 422 417 ! surface pressure gradient 423 418 IF( lk_vvl) THEN 424 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &425 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)426 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &427 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)419 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 420 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 421 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 422 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 428 423 ELSE 429 424 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 432 427 ! enstrophy conserving formulation for planetary vorticity term 433 428 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 434 429 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 435 430 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 436 431 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 437 432 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 438 433 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 439 434 ! after transports 440 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)441 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)435 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 436 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 442 437 END DO 443 438 END DO … … 449 444 ! surface pressure gradient 450 445 IF( lk_vvl) THEN 451 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &452 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)453 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &454 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)446 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 447 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 448 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 449 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 455 450 ELSE 456 451 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 463 458 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 464 459 ! after transports 465 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)466 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)460 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 461 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 467 462 END DO 468 463 END DO … … 477 472 478 473 ! ... Boundary conditions on ua_e, va_e, ssha_e 479 CALL lbc_lnk( ua_e , 'U', -1. )480 CALL lbc_lnk( va_e , 'V', -1. )481 CALL lbc_lnk( ssha_e, 'T', 1. )474 CALL lbc_lnk( zua_e , 'U', -1. ) 475 CALL lbc_lnk( zva_e , 'V', -1. ) 476 CALL lbc_lnk( zssha_e, 'T', 1. ) 482 477 483 478 ! temporal sum 484 479 !------------- 485 IF( jit >= nn_baro / 2 ) THEN 486 zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:) 487 zua_b (:,:) = zua_b (:,:) + ua_e (:,:) 488 zva_b (:,:) = zva_b (:,:) + va_e (:,:) 489 ENDIF 480 zu_sum (:,:) = zu_sum (:,:) + zua_e (:,:) 481 zv_sum (:,:) = zv_sum (:,:) + zva_e (:,:) 482 zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 490 483 491 484 ! Time filter and swap of dynamics arrays 492 485 ! --------------------------------------- 493 486 IF( jit == 1 ) THEN ! Euler (forward) time stepping 494 sshb_e (:,:) = sshn_e(:,:)495 zub_e (:,:) = zun_e(:,:)496 zvb_e (:,:) = zvn_e(:,:)497 sshn_e (:,:) = ssha_e(:,:)498 zun_e (:,:) =ua_e (:,:)499 zvn_e (:,:) =va_e (:,:)487 zsshb_e(:,:) = sshn_e (:,:) 488 zub_e (:,:) = un_e (:,:) 489 zvb_e (:,:) = vn_e (:,:) 490 sshn_e (:,:) = zssha_e(:,:) 491 un_e (:,:) = zua_e (:,:) 492 vn_e (:,:) = zva_e (:,:) 500 493 ELSE ! Asselin filtering 501 sshb_e (:,:) = atfp * ( sshb_e(:,:) +ssha_e(:,:) ) + atfp1 * sshn_e(:,:)502 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)503 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)504 sshn_e (:,:) = ssha_e(:,:)505 zun_e (:,:) =ua_e (:,:)506 zvn_e (:,:) =va_e (:,:)494 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 495 zub_e (:,:) = atfp * ( zub_e (:,:) + zua_e (:,:) ) + atfp1 * un_e (:,:) 496 zvb_e (:,:) = atfp * ( zvb_e (:,:) + zva_e (:,:) ) + atfp1 * vn_e (:,:) 497 sshn_e (:,:) = zssha_e(:,:) 498 un_e (:,:) = zua_e (:,:) 499 vn_e (:,:) = zva_e (:,:) 507 500 ENDIF 508 501 509 IF( lk_vvl ) THEN ! Variable volume 502 IF( lk_vvl ) THEN ! Variable volume ! needed for BDY ??? 510 503 511 504 ! Sea surface elevation … … 528 521 529 522 ! Boundaries conditions 530 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; 531 532 ! Scale factors at before and after time step523 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 524 525 ! Ocean depth at U- and V-points 533 526 ! ------------------------------------------- 534 CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ; CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e ) 535 536 ! Ocean depth at U- and V-points 537 hu_e(:,:) = 0.e0 538 hv_e(:,:) = 0.e0 539 540 DO jk = 1, jpk 541 hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) 542 hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) 543 END DO 544 545 ENDIF ! End variable volume case 546 527 hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 528 hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 529 530 ! 531 ENDIF 547 532 ! ! ==================== ! 548 533 END DO ! end loop ! 549 534 ! ! ==================== ! 550 535 551 552 536 ! Time average of after barotropic variables 553 zcoef = 1.e0 / ( nn_baro + 1 ) 554 zssha_b(:,:) = zcoef * zssha_b(:,:) 555 zua_b (:,:) = zcoef * zua_b (:,:) 556 zva_b (:,:) = zcoef * zva_b (:,:) 537 zcoef = 1.e0 / ( 2 * nn_baro + 1 ) 538 un_e (:,:) = zcoef * zu_sum (:,:) 539 vn_e (:,:) = zcoef * zv_sum (:,:) 540 sshn_e(:,:) = zcoef * zssh_sum(:,:) 541 557 542 #if defined key_obc 558 543 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) … … 561 546 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:) 562 547 #endif 563 564 565 ! ---------------------------------------------------------------------------566 ! Phase 3 : Update sea surface height from time averaged barotropic variables567 ! ---------------------------------------------------------------------------568 569 570 ! Horizontal divergence of time averaged barotropic transports571 !-------------------------------------------------------------572 IF( .NOT. lk_vvl ) THEN573 DO jj = 2, jpjm1574 DO ji = fs_2, fs_jpim1 ! vector opt.575 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) &576 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) &577 & / ( e1t(ji,jj) * e2t(ji,jj) )578 END DO579 END DO580 ENDIF581 582 #if defined key_obc && ! defined key_vvl583 ! open boundaries (div must be zero behind the open boundary)584 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column585 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east586 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west587 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north588 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south589 #endif590 591 #if defined key_bdy592 DO jj = 1, jpj593 DO ji = 1, jpi594 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)595 END DO596 END DO597 #endif598 599 ! sea surface height600 !-------------------601 IF( lk_vvl ) THEN602 sshbb(:,:) = sshb(:,:)603 sshb (:,:) = sshn(:,:)604 sshn (:,:) = ssha(:,:)605 ELSE606 sshb (:,:) = sshn(:,:)607 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)608 ENDIF609 610 ! ... Boundary conditions on sshn611 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )612 613 548 614 549 ! ----------------------------------------------------------------------------- 615 ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)550 ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step) 616 551 ! ----------------------------------------------------------------------------- 617 552 618 ! Swap on time averaged barotropic variables 619 ! ------------------------------------------ 620 sshb_b(:,:) = sshn_b (:,:) 621 IF ( neuler == 0 .AND. kt == nit000 ) zssha_b(:,:) = sshn(:,:) 622 sshn_b(:,:) = zssha_b(:,:) 623 un_b (:,:) = zua_b (:,:) 624 vn_b (:,:) = zva_b (:,:) 625 553 554 626 555 ! add time averaged barotropic coriolis and surface pressure gradient 627 556 ! terms to the general momentum trend 628 557 ! -------------------------------------------------------------------- 629 558 DO jk=1,jpkm1 630 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b631 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b559 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b 560 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b 632 561 END DO 633 562 … … 636 565 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 637 566 638 ! print sum trends (used for debugging)639 IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )640 567 ! 641 568 END SUBROUTINE dyn_spg_ts … … 655 582 ! 656 583 IF( TRIM(cdrw) == 'READ' ) THEN 657 IF( iom_varid( numror, 'sshn ', ldstop = .FALSE. ) > 0 ) THEN658 CALL iom_get( numror, jpdom_autoglo, 'ssh b' , sshb(:,:) )659 CALL iom_get( numror, jpdom_autoglo, ' sshn' , sshn(:,:) )660 IF( neuler == 0 ) sshb(:,:) = sshn(:,:)584 IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 585 CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) ) ! free surface and 586 CALL iom_get( numror, jpdom_autoglo, 'un_e' , un_e (:,:) ) ! horizontal transports issued 587 CALL iom_get( numror, jpdom_autoglo, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 661 588 ELSE 662 IF( nn_rstssh == 1 ) THEN 663 sshb(:,:) = 0.e0 664 sshn(:,:) = 0.e0 665 ENDIF 666 ENDIF 667 IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 668 CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) ) ! free surface issued 669 CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop 670 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! horizontal transports issued 671 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 672 IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) 673 ELSE 674 sshb_b(:,:) = sshb(:,:) 675 sshn_b(:,:) = sshn(:,:) 676 un_b (:,:) = 0.e0 677 vn_b (:,:) = 0.e0 589 sshn_e(:,:) = sshn(:,:) 590 un_e (:,:) = 0.e0 591 vn_e (:,:) = 0.e0 678 592 ! vertical sum 679 593 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 680 594 DO jk = 1, jpkm1 681 595 DO ji = 1, jpij 682 un_ b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)683 vn_ b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)596 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 597 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 684 598 END DO 685 599 END DO 686 600 ELSE ! No vector opt. 687 601 DO jk = 1, jpkm1 688 un_ b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)689 vn_ b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)602 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 603 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 690 604 END DO 691 605 ENDIF 692 606 ENDIF 693 607 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 694 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) ) 695 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) ) 696 CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued 697 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop 698 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued 699 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 608 CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) ) ! free surface and 609 CALL iom_rstput( kt, nitrst, numrow, 'un_e' , un_e (:,:) ) ! horizontal transports issued 610 CALL iom_rstput( kt, nitrst, numrow, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 700 611 ENDIF 701 612 ! -
trunk/NEMO/OPA_SRC/DYN/dynvor.F90
r1152 r1438 5 5 !! planetary vorticity trends 6 6 !!====================================================================== 7 !! History : 1.0 ! 89-12 (P. Andrich) vor_ens: Original code 8 !! 5.0 ! 91-11 (G. Madec) vor_ene, vor_mix: Original code 9 !! 6.0 ! 96-01 (G. Madec) s-coord, suppress work arrays 10 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 11 !! 8.5 ! 04-02 (G. Madec) vor_een: Original code 12 !! 9.0 ! 03-08 (G. Madec) vor_ctl: Original code 13 !! 9.0 ! 05-11 (G. Madec) dyn_vor: Original code (new step architecture) 14 !! 9.0 ! 06-11 (G. Madec) flux form advection: add metric term 7 !! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code 8 !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code 9 !! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays 10 !! 8.5 ! 2002-08 (G. Madec) F90: Free form and module 11 !! NEMO 1.0 ! 2004-02 (G. Madec) vor_een: Original code 12 !! - ! 2003-08 (G. Madec) add vor_ctl 13 !! - ! 2005-11 (G. Madec) add dyn_vor (new step architecture) 14 !! 2.0 ! 2006-11 (G. Madec) flux form advection: add metric term 15 !! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme 15 16 !!---------------------------------------------------------------------- 16 17 … … 37 38 PUBLIC dyn_vor ! routine called by step.F90 38 39 39 ! !* Namelist nam_dynvor: vorticity term40 ! !!* Namelist nam_dynvor: vorticity term 40 41 LOGICAL, PUBLIC :: ln_dynvor_ene = .FALSE. !: energy conserving scheme 41 42 LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme … … 52 53 # include "vectopt_loop_substitute.h90" 53 54 !!---------------------------------------------------------------------- 54 !! OPA 9.0 , LOCEAN-IPSL (2006)55 !! NEMO/OPA 3,2 , LOCEAN-IPSL (2009) 55 56 !! $Id$ 56 57 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 174 175 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 175 176 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 176 177 ! 177 178 END SUBROUTINE dyn_vor 178 179 … … 206 207 INTEGER , INTENT(in ) :: kt ! ocean time-step index 207 208 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 208 !! =nrvm (relative vorticity or metric)209 ! ! =nrvm (relative vorticity or metric) 209 210 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 210 211 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend … … 222 223 ENDIF 223 224 224 ! Local constant initialization 225 zfact2 = 0.5 * 0.5 225 zfact2 = 0.5 * 0.5 ! Local constant initialization 226 226 227 227 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) … … 229 229 DO jk = 1, jpkm1 ! Horizontal slab 230 230 ! ! =============== 231 ! 231 232 ! Potential vorticity and horizontal fluxes 232 233 ! ----------------------------------------- … … 315 316 INTEGER, INTENT(in) :: kt ! ocean timestep index 316 317 !! 317 INTEGER :: ji, jj, jk ! dummy loop indices318 INTEGER :: ji, jj, jk ! dummy loop indices 318 319 REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! temporary scalars 319 320 REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! " " … … 327 328 ENDIF 328 329 329 ! Local constant initialization 330 zfact1 = 0.5 * 0.25 330 zfact1 = 0.5 * 0.25 ! Local constant initialization 331 331 zfact2 = 0.5 * 0.5 332 332 … … 335 335 DO jk = 1, jpkm1 ! Horizontal slab 336 336 ! ! =============== 337 337 ! 338 338 ! Relative and planetary potential vorticity and horizontal fluxes 339 339 ! ---------------------------------------------------------------- … … 438 438 ENDIF 439 439 440 ! Local constant initialization 441 zfact1 = 0.5 * 0.25 440 zfact1 = 0.5 * 0.25 ! Local constant initialization 442 441 443 442 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) … … 445 444 DO jk = 1, jpkm1 ! Horizontal slab 446 445 ! ! =============== 446 ! 447 447 ! Potential vorticity and horizontal fluxes 448 448 ! ----------------------------------------- … … 465 465 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 466 466 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 467 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &467 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 468 468 & ) 469 469 END DO 470 470 END DO 471 471 END SELECT 472 472 ! 473 473 IF( ln_sco ) THEN 474 474 DO jj = 1, jpj ! caution: don't use (:,:) for this loop … … 487 487 END DO 488 488 ENDIF 489 489 ! 490 490 ! Compute and add the vorticity term trend 491 491 ! ---------------------------------------- … … 514 514 !! 515 515 !! ** Method : Trend evaluated using now fields (centered in time) 516 !! and the Arakawa and Lamb (19 XX) flux form formulation : conserves516 !! and the Arakawa and Lamb (1980) flux form formulation : conserves 517 517 !! both the horizontal kinetic energy and the potential enstrophy 518 !! when horizontal divergence is zero. 519 !! The trend of the vorticity term is given by: 520 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 521 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 522 !! Add this trend to the general momentum trend (ua,va): 523 !! (ua,va) = (ua,va) + ( voru , vorv ) 518 !! when horizontal divergence is zero (see the NEMO documentation) 519 !! Add this trend to the general momentum trend (ua,va). 524 520 !! 525 521 !! ** Action : - Update (ua,va) with the now vorticity term trend … … 531 527 INTEGER , INTENT(in ) :: kt ! ocean time-step index 532 528 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 533 !! =nrvm (relative vorticity or metric)529 ! ! =nrvm (relative vorticity or metric) 534 530 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 535 531 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 536 532 !! 537 INTEGER :: ji, jj, jk! dummy loop indices533 INTEGER :: ji, jj, jk ! dummy loop indices 538 534 REAL(wp) :: zfac12, zua, zva ! temporary scalars 539 535 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace 540 536 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace 537 #if defined key_vvl 538 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3f 539 #else 541 540 REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: ze3f 541 #endif 542 542 !!---------------------------------------------------------------------- 543 543 … … 546 546 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 547 547 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 548 548 ENDIF 549 550 IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t) 549 551 DO jk = 1, jpk 550 552 DO jj = 1, jpjm1 … … 559 561 ENDIF 560 562 561 ! Local constant initialization 562 zfac12 = 1.e0 / 12.e0 563 zfac12 = 1.e0 / 12.e0 ! Local constant initialization 563 564 564 565 … … 571 572 ! ----------------------------------------- 572 573 SELECT CASE( kvor ) ! vorticity considered 573 CASE ( 1 ) ; zwz(:,:) = ff(:,:) * ze3f(:,:,jk) ! planetary vorticity (Coriolis) 574 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) ! relative vorticity 574 CASE ( 1 ) ! planetary vorticity (Coriolis) 575 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 576 CASE ( 2 ) ! relative vorticity 577 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 575 578 CASE ( 3 ) ! metric term 576 579 DO jj = 1, jpjm1 … … 581 584 END DO 582 585 END DO 583 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity) 586 CASE ( 4 ) ! total (relative + planetary vorticity) 587 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 584 588 CASE ( 5 ) ! total (coriolis + metric) 585 589 DO jj = 1, jpjm1 … … 588 592 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 589 593 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 590 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &594 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 591 595 & ) * ze3f(ji,jj,jk) 592 596 END DO … … 599 603 ! Compute and add the vorticity term trend 600 604 ! ---------------------------------------- 601 jj =2602 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ;ztsw(1,:) = 0605 jj = 2 606 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 603 607 DO ji = 2, jpi 604 608 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) … … 636 640 !! 637 641 !! ** Purpose : Control the consistency between cpp options for 638 !! tracer advection schemes642 !! tracer advection schemes 639 643 !!---------------------------------------------------------------------- 640 644 INTEGER :: ioptio ! temporary integer -
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
r1241 r1438 1 1 MODULE wzvmod 2 !! MODULE sshwzv 2 3 !!============================================================================== 3 !! *** MODULE wzvmod***4 !! Ocean d iagnostic variable :vertical velocity4 !! *** MODULE sshwzv *** 5 !! Ocean dynamics : sea surface height and vertical velocity 5 6 !!============================================================================== 6 !! History : 5.0 ! 90-10 (C. Levy, G. Madec) Original code7 !! 7.0 ! 96-01 (G. Madec) Statement function for e38 !! 8.5 ! 02-07 (G. Madec) Free form, F90 9 !! " ! 07-07 (D. Storkey) Zero zhdiv at open boundary (BDY)10 !! ----------------------------------------------------------------------11 !! wzv : Compute the vertical velocity12 !! ----------------------------------------------------------------------13 !! * Modules used7 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! ssh_wzv : after ssh & now vertical velocity 12 !! ssh_nxt : filter ans swap the ssh arrays 13 !! ssh_rst : read/write ssh restart fields in the ocean restart file 14 !!---------------------------------------------------------------------- 14 15 USE oce ! ocean dynamics and tracers variables 15 16 USE dom_oce ! ocean space and time domain variables 16 17 USE sbc_oce ! surface boundary condition: ocean 17 18 USE domvvl ! Variable volume 19 USE iom ! I/O library 20 USE restart ! only for lrst_oce 18 21 USE in_out_manager ! I/O manager 19 22 USE prtctl ! Print control 20 23 USE phycst 21 USE bdy_oce ! unstructured open boundaries22 24 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 25 USE obc_par ! open boundary cond. parameter … … 27 29 PRIVATE 28 30 29 !! * Routine accessibility30 PUBLIC wzv ! routine called by step.F90 and inidtr.F9031 PUBLIC ssh_wzv ! called by step.F90 32 PUBLIC ssh_nxt ! called by step.F90 31 33 32 34 !! * Substitutions 33 35 # include "domzgr_substitute.h90" 34 !!---------------------------------------------------------------------- 35 !! OPA 9.0 , LOCEAN-IPSL (2005) 36 # include "vectopt_loop_substitute.h90" 37 38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 36 40 !! $Id$ 37 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 40 44 CONTAINS 41 45 42 SUBROUTINE wzv( kt )43 !!---------------------------------------------------------------------- 44 !! *** ROUTINEwzv ***45 !! 46 !! ** Purpose : Compute the now vertical velocity after the array swap47 !! 48 !! ** Method : Using the incompressibility hypothesis, the vertical49 !! velocity is computed by integrating the horizontal divergence50 !! from the bottom to the surface.51 !! The boundary conditions are w=0 at the bottom (no flux) and,52 !! in regid-lid case, w=0 at the sea surface.53 !! 54 !! ** action : wn array : the now vertical velocity55 !! ----------------------------------------------------------------------56 !! * Arguments57 INTEGER, INTENT( in ) :: kt ! ocean time-step index58 59 !! * Local declarations60 INTEGER :: jk ! dummy loop indices61 !! Variable volume62 INTEGER :: ji, jj ! dummy loop indices63 REAL(wp) :: z2dt, zraur ! temporary scalar64 REAL(wp), DIMENSION (jpi,jpj) :: zssha, zun, zvn, zhdiv65 #if defined key_bdy 66 INTEGER :: jgrd, jb! temporary scalars67 #endif 46 SUBROUTINE ssh_wzv( kt ) 47 !!---------------------------------------------------------------------- 48 !! *** ROUTINE ssh_wzv *** 49 !! 50 !! ** Purpose : compute the after ssh (ssha), the now vertical velocity 51 !! and update the now vertical coordinate (lk_vvl=T). 52 !! 53 !! ** Method : - 54 !! 55 !! - Using the incompressibility hypothesis, the vertical 56 !! velocity is computed by integrating the horizontal divergence 57 !! from the bottom to the surface minus the scale factor evolution. 58 !! The boundary conditions are w=0 at the bottom (no flux) and. 59 !! 60 !! ** action : ssha : after sea surface height 61 !! wn : now vertical velocity 62 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 63 !! at u-, v-, f-point s 64 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 65 !!---------------------------------------------------------------------- 66 INTEGER, INTENT(in) :: kt ! time step 67 !! 68 INTEGER :: ji, jj, jk ! dummy loop indices 69 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 70 REAL(wp) :: z2dt, zraur ! temporary scalars 71 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 68 72 !!---------------------------------------------------------------------- 69 73 70 74 IF( kt == nit000 ) THEN 71 75 IF(lwp) WRITE(numout,*) 72 IF(lwp) WRITE(numout,*) 'wzv : vertical velocity from continuity eq.' 73 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 74 75 ! bottom boundary condition: w=0 (set once for all) 76 wn(:,:,jpk) = 0.e0 77 ENDIF 78 79 IF( lk_vvl ) THEN ! Variable volume 80 ! 81 z2dt = 2. * rdt ! time step: leap-frog 82 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 83 zraur = 1. / rauw 84 85 ! Vertically integrated quantities 86 ! -------------------------------- 87 zun(:,:) = 0.e0 88 zvn(:,:) = 0.e0 89 ! 90 DO jk = 1, jpkm1 ! Vertically integrated transports (now) 91 zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 92 zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 93 END DO 94 95 ! Horizontal divergence of barotropic transports 96 !-------------------------------------------------- 97 zhdiv(:,:) = 0.e0 98 DO jj = 2, jpjm1 99 DO ji = 2, jpim1 ! vector opt. 100 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) & 101 & - e2u(ji-1,jj ) * zun(ji-1,jj ) & 102 & + e1v(ji ,jj ) * zvn(ji ,jj ) & 103 & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) & 104 & / ( e1t(ji,jj) * e2t(ji,jj) ) 76 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 77 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 78 ! 79 CALL ssh_rst( nit000, 'READ' ) ! read or initialize sshb and sshn 80 ! 81 wn(:,:,jpk) = 0.e0 ! bottom boundary condition: w=0 (set once for all) 82 ! 83 IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only) 84 DO jj = 1, jpjm1 85 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible 86 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 87 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 88 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 89 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 90 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 91 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 92 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 93 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) & 94 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 95 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 96 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 97 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 98 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 99 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 100 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 101 END DO 102 END DO 103 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 104 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 105 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 106 ENDIF 107 ! 108 ENDIF 109 110 ! set time step size (Euler/Leapfrog) 111 z2dt = 2. * rdt 112 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 113 114 zraur = 1. / rauw 115 116 ! !------------------------------! 117 ! ! After Sea Surface Height ! 118 ! !------------------------------! 119 zhdiv(:,:) = 0.e0 120 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 121 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 122 END DO 123 124 ! ! Sea surface elevation time stepping 125 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 126 127 #if defined key_obc 128 # if defined key_agrif 129 IF ( Agrif_Root() ) THEN 130 # endif 131 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 132 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 133 # if defined key_agrif 134 ENDIF 135 # endif 136 #endif 137 138 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 139 IF( lk_vvl ) THEN ! (required only in key_vvl case) 140 DO jj = 1, jpjm1 141 DO ji = 1, fs_jpim1 ! Vector Opt. 142 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 143 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 144 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 145 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 146 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 147 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 148 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) & ! Caution : fmask not used 149 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) & 150 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) ) 105 151 END DO 106 152 END DO 107 108 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 109 ! open boundaries (div must be zero behind the open boundary) 110 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 111 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east 112 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west 113 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 114 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south 115 #endif 116 117 #if defined key_bdy 118 jgrd=1 !: tracer grid. 119 DO jb = 1, nblenrim(jgrd) 120 ji = nbi(jb,jgrd) 121 jj = nbj(jb,jgrd) 122 zhdiv(ji,jj) = 0.e0 153 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 154 CALL lbc_lnk( sshv_a, 'V', 1. ) 155 CALL lbc_lnk( sshf_a, 'F', 1. ) 156 ENDIF 157 158 ! !------------------------------! 159 ! ! Now Vertical Velocity ! 160 ! !------------------------------! 161 ! ! integrate from the bottom the hor. divergence 162 DO jk = jpkm1, 1, -1 163 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 164 & - ( fse3t_a(:,:,jk) & 165 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 166 END DO 167 168 ! !------------------------------! 169 ! ! Update Now Vertical coord. ! 170 ! !------------------------------! 171 IF( lk_vvl ) THEN ! only in vvl case) 172 ! ! now local depth and scale factors (stored in fse3. arrays) 173 DO jk = 1, jpkm1 174 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! depths 175 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 176 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 177 ! 178 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors 179 fse3u (:,:,jk) = fse3u_n (:,:,jk) 180 fse3v (:,:,jk) = fse3v_n (:,:,jk) 181 fse3f (:,:,jk) = fse3f_n (:,:,jk) 182 fse3w (:,:,jk) = fse3w_n (:,:,jk) 183 fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 184 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 123 185 END DO 124 #endif 125 126 CALL lbc_lnk( zhdiv, 'T', 1. ) 127 128 ! Sea surface elevation time stepping 129 ! ----------------------------------- 130 zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 131 132 ! Vertical velocity computed from bottom 133 ! -------------------------------------- 134 DO jk = jpkm1, 1, -1 135 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 136 & - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt 137 END DO 138 139 ELSE ! Fixed volume 140 141 ! Vertical velocity computed from bottom 142 ! -------------------------------------- 143 DO jk = jpkm1, 1, -1 144 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 145 END DO 146 147 ENDIF 148 149 IF(ln_ctl) CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 - : ', mask1=wn) 150 151 END SUBROUTINE wzv 186 ! ! ocean depth (at u- and v-points) 187 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 188 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 189 ! ! masked inverse of the ocean depth (at u- and v-points) 190 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 191 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 192 193 ENDIF 194 ! 195 END SUBROUTINE ssh_wzv 196 197 198 SUBROUTINE ssh_nxt( kt ) 199 !!---------------------------------------------------------------------- 200 !! *** ROUTINE ssh_nxt *** 201 !! 202 !! ** Purpose : achieve the sea surface height time stepping by 203 !! applying Asselin time filter and swapping the arrays 204 !! ssha already computed in ssh_wzv 205 !! 206 !! ** Method : - apply Asselin time fiter to now ssh and swap : 207 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 208 !! sshn = ssha 209 !! 210 !! ** action : - sshb, sshn : before & now sea surface height 211 !! ready for the next time step 212 !!---------------------------------------------------------------------- 213 INTEGER, INTENT( in ) :: kt ! ocean time-step index 214 !! 215 INTEGER :: ji, jj ! dummy loop indices 216 !!---------------------------------------------------------------------- 217 218 IF( kt == nit000 ) THEN 219 IF(lwp) WRITE(numout,*) 220 IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 221 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 222 ENDIF 223 224 ! Time filter and swap of the ssh 225 ! ------------------------------- 226 227 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 228 ! ! ---------------------- ! 229 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 230 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 231 sshu_n(:,:) = sshu_a(:,:) 232 sshv_n(:,:) = sshv_a(:,:) 233 sshf_n(:,:) = sshf_a(:,:) 234 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 235 DO jj = 1, jpj 236 DO ji = 1, jpi ! before <-- now filtered 237 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 238 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 239 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 240 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 241 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 242 sshu_n(ji,jj) = sshu_a(ji,jj) 243 sshv_n(ji,jj) = sshv_a(ji,jj) 244 sshf_n(ji,jj) = sshf_a(ji,jj) 245 END DO 246 END DO 247 ENDIF 248 ! 249 ELSE ! fixed levels : ssh at t-point only 250 ! ! ------------ ! 251 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 252 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 253 ! 254 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 255 DO jj = 1, jpj 256 DO ji = 1, jpi ! before <-- now filtered 257 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 258 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 259 END DO 260 END DO 261 ENDIF 262 ! 263 ENDIF 264 265 ! ! write filtered free surface arrays in restart file 266 IF( lrst_oce ) CALL ssh_rst( kt, 'WRITE' ) 267 ! 268 IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 269 ! 270 END SUBROUTINE ssh_nxt 271 272 273 SUBROUTINE ssh_rst( kt, cdrw ) 274 !!--------------------------------------------------------------------- 275 !! *** ROUTINE ssh_rst *** 276 !! 277 !! ** Purpose : Read or write Sea Surface Height arrays in restart file 278 !! 279 !! ** action : - cdrw = READ : sshb, sshn read in ocean restart file 280 !! - cdrw = WRITE : sshb, sshn written in ocean restart file 281 !!---------------------------------------------------------------------- 282 INTEGER , INTENT(in) :: kt ! ocean time-step 283 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 284 !!---------------------------------------------------------------------- 285 ! 286 IF( TRIM(cdrw) == 'READ' ) THEN 287 IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 288 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) ) 289 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) ) 290 IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 291 ELSE 292 IF( nn_rstssh == 1 ) THEN 293 sshb(:,:) = 0.e0 294 sshn(:,:) = 0.e0 295 ENDIF 296 ENDIF 297 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 298 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb(:,:) ) 299 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn(:,:) ) 300 ENDIF 301 ! 302 END SUBROUTINE ssh_rst 152 303 153 304 !!====================================================================== -
trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90
r1146 r1438 94 94 ! of the square root of the resulting N^2 ( required to compute 95 95 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 96 zn2 = MAX( rn2 (ji,1,jk), 0.e0 )96 zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 97 97 zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 98 98 ! Compute elements required for the inverse time scale of baroclinic … … 113 113 ! of the square root of the resulting N^2 ( required to compute 114 114 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 115 zn2 = MAX( rn2 (ji,jj,jk), 0.e0 )115 zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 116 116 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 117 117 ! Compute elements required for the inverse time scale of baroclinic -
trunk/NEMO/OPA_SRC/TRA/tradmp.F90
r1415 r1438 331 331 INTEGER :: ji, jj, jk ! dummy loop indices 332 332 INTEGER :: ii0, ii1, ij0, ij1 ! " " 333 INTEGER :: i dmp! logical unit for file restoring damping term333 INTEGER :: inum0 ! logical unit for file restoring damping term 334 334 INTEGER :: icot ! logical unit for file distance to the coast 335 335 REAL(wp) :: zinfl, zlon ! temporary scalars -
trunk/NEMO/OPA_SRC/TRA/tranxt.F90
r1146 r1438 14 14 !! 2.0 ! 2006-02 (L. Debreu, C. Mazauric) Agrif implementation 15 15 !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf 16 !! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option 16 17 !!---------------------------------------------------------------------- 17 18 18 19 !!---------------------------------------------------------------------- 19 20 !! tra_nxt : time stepping on temperature and salinity 21 !! tra_nxt_fix : time stepping on temperature and salinity : fixed volume case 22 !! tra_nxt_vvl : time stepping on temperature and salinity : variable volume case 20 23 !!---------------------------------------------------------------------- 21 24 USE oce ! ocean dynamics and tracers variables 22 25 USE dom_oce ! ocean space and time domain variables 23 26 USE zdf_oce ! ??? 27 USE domvvl ! variable volume 24 28 USE dynspg_oce ! surface pressure gradient variables 25 29 USE trdmod_oce ! ocean variables trends … … 39 43 PUBLIC tra_nxt ! routine called by step.F90 40 44 45 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 46 41 47 !! * Substitutions 42 48 # include "domzgr_substitute.h90" … … 67 73 !! at the AGRIF zoom boundaries (lk_agrif=T) 68 74 !! 69 !! - Apply the Asselin time filter on now fields, 70 !! save in (ta,sa) an average over the three time levels 71 !! which will be used to compute rdn and thus the semi-implicit 72 !! hydrostatic pressure gradient (ln_dynhpg_imp = T), and 73 !! swap tracer fields to prepare the next time_step. 74 !! This can be summurized for tempearture as: 75 !! zt = (ta+2tn+tb)/4 ln_dynhpg_imp = T 76 !! zt = 0 otherwise 77 !! tb = tn + atfp*[ tb - 2 tn + ta ] 78 !! tn = ta 79 !! ta = zt (NB: reset to 0 after eos_bn2 call) 75 !! - Update lateral boundary conditions on AGRIF children 76 !! domains (lk_agrif=T) 80 77 !! 81 78 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step … … 87 84 INTEGER, INTENT(in) :: kt ! ocean time-step index 88 85 !! 89 INTEGER :: j i, jj, jk ! dummy loop indices90 REAL(wp) :: z t, zs, zfact ! temporary scalars86 INTEGER :: jk ! dummy loop indices 87 REAL(wp) :: zfact ! temporary scalars 91 88 !!---------------------------------------------------------------------- 92 89 … … 111 108 CALL Agrif_tra ! AGRIF zoom boundaries 112 109 #endif 110 111 ! set time step size (Euler/Leapfrog) 112 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt_t(:) = rdttra(:) ! at nit000 (Euler) 113 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt_t(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog) 114 ENDIF 113 115 114 116 ! trends computation initialisation … … 118 120 ENDIF 119 121 120 ! Asselin time filter and swap of arrays 121 ! 122 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler 1st time step : swap only 123 DO jk = 1, jpkm1 124 tb(:,:,jk) = tn(:,:,jk) ! ta, sa remain at their values which 125 sb(:,:,jk) = sn(:,:,jk) ! correspond to tn, sn after the sawp 126 tn(:,:,jk) = ta(:,:,jk) 127 sn(:,:,jk) = sa(:,:,jk) 128 END DO 129 ! 130 ELSE ! Leap-Frog : filter + swap 131 ! 132 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case 133 DO jk = 1, jpkm1 ! (save the averaged of the 3 time steps 134 DO jj = 1, jpj ! in the after fields) 135 DO ji = 1, jpi 136 zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 137 zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 138 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 139 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 140 tn(ji,jj,jk) = ta(ji,jj,jk) 141 sn(ji,jj,jk) = sa(ji,jj,jk) 142 ta(ji,jj,jk) = zt 143 sa(ji,jj,jk) = zs 144 END DO 145 END DO 146 END DO 147 ELSE ! explicit hpg case 148 DO jk = 1, jpkm1 149 DO jj = 1, jpj 150 DO ji = 1, jpi 151 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 152 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 153 tn(ji,jj,jk) = ta(ji,jj,jk) 154 sn(ji,jj,jk) = sa(ji,jj,jk) 155 END DO 156 END DO 157 END DO 158 ENDIF 159 ! 122 ! Leap-Frog + Asselin filter time stepping 123 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt ) ! variable volume level (vvl) 124 ELSE ; CALL tra_nxt_fix( kt ) ! fixed volume level 160 125 ENDIF 161 126 … … 165 130 #endif 166 131 167 ! trends diagnostics : Asselin filter trend : (tb filtered - tb)/2dt168 IF( l_trdtra ) THEN 132 ! trends computation 133 IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 169 134 DO jk = 1, jpkm1 170 zfact = 1.e0 / ( 2.*rdttra(jk) ) ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK135 zfact = 1.e0 / r2dt_t(jk) 171 136 ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 172 137 ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact … … 174 139 CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) 175 140 END IF 141 176 142 ! ! control print 177 143 IF(ln_ctl) CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt - Tn: ', mask1=tmask, & … … 180 146 END SUBROUTINE tra_nxt 181 147 148 149 SUBROUTINE tra_nxt_fix( kt ) 150 !!---------------------------------------------------------------------- 151 !! *** ROUTINE tra_nxt_fix *** 152 !! 153 !! ** Purpose : fixed volume: apply the Asselin time filter and 154 !! swap the tracer fields. 155 !! 156 !! ** Method : - Apply a Asselin time filter on now fields. 157 !! - save in (ta,sa) an average over the three time levels 158 !! which will be used to compute rdn and thus the semi-implicit 159 !! hydrostatic pressure gradient (ln_dynhpg_imp = T) 160 !! - swap tracer fields to prepare the next time_step. 161 !! This can be summurized for tempearture as: 162 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 163 !! ztm = 0 otherwise 164 !! tb = tn + atfp*[ tb - 2 tn + ta ] 165 !! tn = ta 166 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 167 !! 168 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 169 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 170 !!---------------------------------------------------------------------- 171 INTEGER, INTENT(in) :: kt ! ocean time-step index 172 !! 173 INTEGER :: ji, jj, jk ! dummy loop indices 174 REAL(wp) :: ztm, ztf ! temporary scalars 175 REAL(wp) :: zsm, zsf ! - - 176 !!---------------------------------------------------------------------- 177 178 IF( kt == nit000 ) THEN 179 IF(lwp) WRITE(numout,*) 180 IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' 181 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 182 ENDIF 183 ! 184 ! ! ----------------------- ! 185 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 186 ! ! ----------------------- ! 187 ! 188 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 189 DO jk = 1, jpkm1 ! (only swap) 190 DO jj = 1, jpj 191 DO ji = 1, jpi 192 tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn 193 sn(ji,jj,jk) = sa(ji,jj,jk) 194 END DO 195 END DO 196 END DO 197 ELSE ! general case (Leapfrog + Asselin filter 198 DO jk = 1, jpkm1 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t 202 zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 203 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 204 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 205 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 206 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 207 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 208 sn(ji,jj,jk) = sa(ji,jj,jk) 209 ta(ji,jj,jk) = ztm ! ta <-- mean t 210 sa(ji,jj,jk) = zsm 211 END DO 212 END DO 213 END DO 214 ENDIF 215 ! ! ----------------------- ! 216 ELSE ! explicit hpg case ! 217 ! ! ----------------------- ! 218 ! 219 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 220 DO jk = 1, jpkm1 221 DO jj = 1, jpj 222 DO ji = 1, jpi 223 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 224 sn(ji,jj,jk) = sa(ji,jj,jk) 225 END DO 226 END DO 227 END DO 228 ELSE ! general case (Leapfrog + Asselin filter 229 DO jk = 1, jpkm1 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 233 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 234 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 235 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 236 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 237 sn(ji,jj,jk) = sa(ji,jj,jk) 238 END DO 239 END DO 240 END DO 241 ENDIF 242 ENDIF 243 ! 244 END SUBROUTINE tra_nxt_fix 245 246 247 SUBROUTINE tra_nxt_vvl( kt ) 248 !!---------------------------------------------------------------------- 249 !! *** ROUTINE tra_nxt_vvl *** 250 !! 251 !! ** Purpose : Time varying volume: apply the Asselin time filter 252 !! and swap the tracer fields. 253 !! 254 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 255 !! - save in (ta,sa) a thickness weighted average over the three 256 !! time levels which will be used to compute rdn and thus the semi- 257 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 258 !! - swap tracer fields to prepare the next time_step. 259 !! This can be summurized for tempearture as: 260 !! ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T 261 !! /(e3t_a +2*e3t_n +e3t_b ) 262 !! ztm = 0 otherwise 263 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 264 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 265 !! tn = ta 266 !! ta = zt (NB: reset to 0 after eos_bn2 call) 267 !! 268 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 269 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 270 !!---------------------------------------------------------------------- 271 INTEGER, INTENT(in) :: kt ! ocean time-step index 272 !! 273 INTEGER :: ji, jj, jk ! dummy loop indices 274 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 275 REAL(wp) :: zsm , zsc_f , zsf , zsca, zscn, zscb ! - - 276 REAL(wp) :: ze3mr, ze3fr ! - - 277 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 278 !!---------------------------------------------------------------------- 279 280 IF( kt == nit000 ) THEN 281 IF(lwp) WRITE(numout,*) 282 IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping' 283 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 284 ENDIF 285 286 ! ! ----------------------- ! 287 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 288 ! ! ----------------------- ! 289 ! 290 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 291 DO jk = 1, jpkm1 ! (only swap) 292 DO jj = 1, jpj 293 DO ji = 1, jpi 294 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 295 sn(ji,jj,jk) = sa(ji,jj,jk) 296 END DO 297 END DO 298 END DO 299 ELSE 300 DO jk = 1, jpkm1 301 DO jj = 1, jpj 302 DO ji = 1, jpi 303 ze3t_b = fse3t_b(ji,jj,jk) 304 ze3t_n = fse3t_n(ji,jj,jk) 305 ze3t_a = fse3t_a(ji,jj,jk) 306 ! ! tracer content at Before, now and after 307 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 308 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 309 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 310 ! 311 ! ! Asselin filter on thickness and tracer content 312 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 313 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 314 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 315 ! 316 ! ! filtered tracer including the correction 317 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 318 ztf = ze3fr * ( ztcn + ztc_f ) 319 zsf = ze3fr * ( zscn + zsc_f ) 320 ! ! mean thickness and tracer 321 ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 322 ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) 323 zsm = ze3mr * ( zsca + 2.* zscn + zscb ) 324 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 325 !!gm e3t_m(ji,jj,jk) = 0.25 / ze3mr 326 ! ! swap of arrays 327 tb(ji,jj,jk) = ztf ! tb <-- tn + filter 328 sb(ji,jj,jk) = zsf 329 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 330 sn(ji,jj,jk) = sa(ji,jj,jk) 331 ta(ji,jj,jk) = ztm ! ta <-- mean t 332 sa(ji,jj,jk) = zsm 333 END DO 334 END DO 335 END DO 336 ENDIF 337 ! ! ----------------------- ! 338 ELSE ! explicit hpg case ! 339 ! ! ----------------------- ! 340 ! 341 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 342 DO jk = 1, jpkm1 ! No filter nor thickness weighting computation required 343 DO jj = 1, jpj ! ONLY swap 344 DO ji = 1, jpi 345 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 346 sn(ji,jj,jk) = sa(ji,jj,jk) 347 END DO 348 END DO 349 END DO 350 ! ! general case (Leapfrog + Asselin filter) 351 ELSE ! apply filter on thickness weighted tracer and swap 352 DO jk = 1, jpkm1 353 DO jj = 1, jpj 354 DO ji = 1, jpi 355 ze3t_b = fse3t_b(ji,jj,jk) 356 ze3t_n = fse3t_n(ji,jj,jk) 357 ze3t_a = fse3t_a(ji,jj,jk) 358 ! ! tracer content at Before, now and after 359 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 360 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 361 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 362 ! 363 ! ! Asselin filter on thickness and tracer content 364 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 365 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 366 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 367 ! 368 ! ! filtered tracer including the correction 369 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 370 ztf = ( ztcn + ztc_f ) * ze3fr 371 zsf = ( zscn + zsc_f ) * ze3fr 372 ! ! swap of arrays 373 tb(ji,jj,jk) = ztf ! tb <-- tn filtered 374 sb(ji,jj,jk) = zsf 375 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 376 sn(ji,jj,jk) = sa(ji,jj,jk) 377 END DO 378 END DO 379 END DO 380 ENDIF 381 ENDIF 382 ! 383 END SUBROUTINE tra_nxt_vvl 384 182 385 !!====================================================================== 183 386 END MODULE tranxt -
trunk/NEMO/OPA_SRC/TRA/trazdf.F90
r1239 r1438 78 78 ztrds(:,:,:) = sa(:,:,:) 79 79 ENDIF 80 81 IF( lk_vvl ) CALL dom_vvl_ssh( kt ) ! compute ssha field82 80 83 81 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend -
trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90
r1146 r1438 42 42 # include "vectopt_loop_substitute.h90" 43 43 !!---------------------------------------------------------------------- 44 !! NEMO/OPA 3. 0 , LOCEAN-IPSL (2008)44 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 45 45 !! $Id$ 46 46 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 78 78 INTEGER :: ji, jj, jk, jl ! dummy loop indices 79 79 REAL(wp) :: zlavmr, zave3r, ze3tr ! temporary scalars 80 REAL(wp) :: zta, zsa, ze3tb , zcoef! temporary scalars80 REAL(wp) :: zta, zsa, ze3tb ! temporary scalars 81 81 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwz, zww ! 3D workspace 82 82 !!--------------------------------------------------------------------- … … 105 105 DO jj = 2, jpjm1 106 106 DO ji = fs_2, fs_jpim1 ! vector opt. 107 zave3r = 1.e0 / fse3w (ji,jj,jk)107 zave3r = 1.e0 / fse3w_n(ji,jj,jk) 108 108 zwy(ji,jj,jk) = avt(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) * zave3r 109 109 zww(ji,jj,jk) = fsavs(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) * zave3r … … 115 115 DO jj = 2, jpjm1 116 116 DO ji = fs_2, fs_jpim1 ! vector opt. 117 ze3tr = zlavmr / fse3t (ji,jj,jk)117 ze3tr = zlavmr / fse3t_n(ji,jj,jk) 118 118 zwx(ji,jj,jk) = zwx(ji,jj,jk) + p2dt(jk) * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr 119 119 zwz(ji,jj,jk) = zwz(ji,jj,jk) + p2dt(jk) * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) * ze3tr … … 130 130 DO jj = 2, jpjm1 131 131 DO ji = fs_2, fs_jpim1 ! vector opt. 132 ze3tb = fs ve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )! before e3t132 ze3tb = fse3t_b(ji,jj,jk) / fse3t(ji,jj,jk) ! before e3t 133 133 zta = zwx(ji,jj,jk) - tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ! total trends * 2*rdt 134 134 zsa = zwz(ji,jj,jk) - sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) 135 zcoef = 1.e0 / fse3t(ji,jj,jk) * tmask(ji,jj,jk) 136 ta(ji,jj,jk) = ( ze3tb * tb(ji,jj,jk) + fse3t(ji,jj,jk) * zta ) * zcoef 137 sa(ji,jj,jk) = ( ze3tb * sb(ji,jj,jk) + fse3t(ji,jj,jk) * zsa ) * zcoef 135 ta(ji,jj,jk) = ( ze3tb * tb(ji,jj,jk) + zta ) * tmask(ji,jj,jk) 136 sa(ji,jj,jk) = ( ze3tb * sb(ji,jj,jk) + zsa ) * tmask(ji,jj,jk) 138 137 END DO 139 138 END DO … … 143 142 DO jj = 2, jpjm1 144 143 DO ji = fs_2, fs_jpim1 ! vector opt. 145 ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) * tmask(ji,jj,jk)146 sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) * tmask(ji,jj,jk)144 ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) * tmask(ji,jj,jk) 145 sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) * tmask(ji,jj,jk) 147 146 END DO 148 147 END DO -
trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r1156 r1438 1 1 MODULE trazdf_imp 2 !!====================================================================== ========2 !!====================================================================== 3 3 !! *** MODULE trazdf_imp *** 4 4 !! Ocean active tracers: vertical component of the tracer mixing trend 5 !!============================================================================== 6 !! History : 7 !! 6.0 ! 90-10 (B. Blanke) Original code 8 !! 7.0 ! 91-11 (G. Madec) 9 !! ! 92-06 (M. Imbard) correction on tracer trend loops 10 !! ! 96-01 (G. Madec) statement function for e3 11 !! ! 97-05 (G. Madec) vertical component of isopycnal 12 !! ! 97-07 (G. Madec) geopotential diffusion in s-coord 13 !! ! 00-08 (G. Madec) double diffusive mixing 14 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 15 !! 9.0 ! 06-11 (G. Madec) New step reorganisation 5 !!====================================================================== 6 !! History : OPA ! 1990-10 (B. Blanke) Original code 7 !! 7.0 ! 1991-11 (G. Madec) 8 !! ! 1992-06 (M. Imbard) correction on tracer trend loops 9 !! ! 1996-01 (G. Madec) statement function for e3 10 !! ! 1997-05 (G. Madec) vertical component of isopycnal 11 !! ! 1997-07 (G. Madec) geopotential diffusion in s-coord 12 !! ! 2000-08 (G. Madec) double diffusive mixing 13 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 14 !! 2.0 ! 2006-11 (G. Madec) New step reorganisation 15 !! 3.2 ! 2009-03 (G. Madec) heat and salt content trends 16 !!---------------------------------------------------------------------- 17 16 18 !!---------------------------------------------------------------------- 17 19 !! tra_zdf_imp : Update the tracer trend with the diagonal vertical 18 20 !! part of the mixing tensor. 19 21 !!---------------------------------------------------------------------- 20 !! * Modules used21 22 USE oce ! ocean dynamics and tracers variables 22 23 USE dom_oce ! ocean space and time domain variables … … 36 37 PRIVATE 37 38 38 !! * Routine accessibility 39 PUBLIC tra_zdf_imp ! routine called by step.F90 39 PUBLIC tra_zdf_imp ! routine called by step.F90 40 40 41 41 !! * Substitutions … … 45 45 # include "vectopt_loop_substitute.h90" 46 46 !!---------------------------------------------------------------------- 47 !!---------------------------------------------------------------------- 48 !! OPA 9.0 , LOCEAN-IPSL (2005) 47 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 49 48 !! $Id$ 50 49 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 90 89 !! 91 90 !!--------------------------------------------------------------------- 92 !! * Modules used 93 USE oce , ONLY : zwd => ua, & ! ua used as workspace 94 zws => va ! va " " 95 !! * Arguments 96 INTEGER, INTENT( in ) :: kt ! ocean time-step index 97 REAL(wp), DIMENSION(jpk), INTENT( in ) :: & 98 p2dt ! vertical profile of tracer time-step 99 100 !! * Local declarations 101 INTEGER :: ji, jj, jk ! dummy loop indices 102 REAL(wp) :: zavi, zrhs, znvvl, & ! temporary scalars 103 ze3tb, ze3tn, ze3ta, zvsfvvl ! variable vertical scale factors 104 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 105 zwi, zwt, zavsi ! workspace arrays 91 USE oce , ONLY : zwd => ua ! ua used as workspace 92 USE oce , ONLY : zws => va ! va - - 93 !! 94 INTEGER , INTENT(in) :: kt ! ocean time-step index 95 REAL(wp), DIMENSION(jpk), INTENT(in) :: p2dt ! vertical profile of tracer time-step 96 !! 97 INTEGER :: ji, jj, jk ! dummy loop indices 98 REAL(wp) :: zavi, zrhs, znvvl ! temporary scalars 99 REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors 100 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt, zavsi ! workspace arrays 106 101 !!--------------------------------------------------------------------- 107 102 … … 169 164 DO jj = 2, jpjm1 170 165 DO ji = fs_2, fs_jpim1 ! vector opt. 171 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) ) 172 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 173 ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl ! now scale factor at T-point 166 ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point 167 & + znvvl * fse3t_a(ji,jj,jk) 168 ze3tn = znvvl & ! now scale factor at T-point 169 & + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 174 170 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 175 171 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) … … 182 178 DO jj = 2, jpjm1 183 179 DO ji = fs_2, fs_jpim1 ! vector opt. 184 z vsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )185 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point180 ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point 181 & + znvvl * fse3t_a(ji,jj,1) 186 182 zwi(ji,jj,1) = 0.e0 187 183 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) … … 227 223 DO jj = 2, jpjm1 228 224 DO ji = fs_2, fs_jpim1 229 zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) ) 230 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl 231 ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1) 225 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 226 ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 232 227 ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1) 233 228 END DO … … 236 231 DO jj = 2, jpjm1 237 232 DO ji = fs_2, fs_jpim1 238 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) ) 239 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl 240 ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk) 233 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) 234 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk) 241 235 zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk) ! zrhs=right hand side 242 ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1)236 ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1) 243 237 END DO 244 238 END DO … … 271 265 DO jj = 2, jpjm1 272 266 DO ji = fs_2, fs_jpim1 ! vector opt. 273 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) ) 274 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 275 ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl ! now scale factor at T-point 267 ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point 268 & + znvvl * fse3t_a(ji,jj,jk) 269 ze3tn = znvvl & ! now scale factor at T-point 270 & + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 276 271 zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 277 272 zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) … … 284 279 DO jj = 2, jpjm1 285 280 DO ji = fs_2, fs_jpim1 ! vector opt. 286 zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) ) 287 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 281 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) 288 282 zwi(ji,jj,1) = 0.e0 289 283 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) … … 328 322 DO jj = 2, jpjm1 329 323 DO ji = fs_2, fs_jpim1 330 z vsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )331 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl ! before scale factor at T-point332 ze3tn = ( 1. - znvvl ) + znvvl *fse3t(ji,jj,1)! now scale factor at T-point324 ze3tb = ( 1. - znvvl ) & ! before scale factor at T-point 325 & + znvvl * fse3t_b(ji,jj,1) 326 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,1) ! now scale factor at T-point 333 327 sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1) 334 328 END DO … … 337 331 DO jj = 2, jpjm1 338 332 DO ji = fs_2, fs_jpim1 339 z vsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )340 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl ! before scale factor at T-point341 ze3tn = ( 1. - znvvl ) + znvvl *fse3t(ji,jj,jk)! now scale factor at T-point333 ze3tb = ( 1. - znvvl ) & ! before scale factor at T-point 334 & + znvvl * fse3t_b(ji,jj,jk) 335 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk) ! now scale factor at T-point 342 336 zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk) ! zrhs=right hand side 343 337 sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) … … 361 355 END DO 362 356 END DO 363 357 ! 364 358 END SUBROUTINE tra_zdf_imp 365 359 -
trunk/NEMO/OPA_SRC/ZDF/zdfevd.F90
r1152 r1438 5 5 !! of vertical eddy mixing coefficient 6 6 !!====================================================================== 7 !! History : OPA ! 1997-06 (G. Madec, A. Lazar) Original code 8 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 9 !! - ! 2005-06 (C. Ethe) KPP parameterization 10 !! 3.2 ! 2009-03 (M. Leclair, G. Madec, R. Benshila) test on both before & after 11 !!---------------------------------------------------------------------- 7 12 8 13 !!---------------------------------------------------------------------- 9 !! zdf_evd : update momentum and tracer Kz at the location of 10 !! statically unstable portion of the water column 11 !! (called if ln_zdfevd=T) 14 !! zdf_evd : increase the momentum and tracer Kz at the location of 15 !! statically unstable portion of the water column (ln_zdfevd=T) 12 16 !!---------------------------------------------------------------------- 13 !! * Modules used14 17 USE oce ! ocean dynamics and tracers variables 15 18 USE dom_oce ! ocean space and time domain variables … … 22 25 PRIVATE 23 26 24 !! * Routine accessibility 25 PUBLIC zdf_evd ! called by step.F90 27 PUBLIC zdf_evd ! called by step.F90 26 28 27 29 !! * Substitutions 28 30 # include "domzgr_substitute.h90" 29 31 !!---------------------------------------------------------------------- 30 !! OPA 9.0 , LOCEAN-IPSL (2005)31 !! $Id$ 32 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt32 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 33 !! $Id$ 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 33 35 !!---------------------------------------------------------------------- 34 36 … … 48 50 !! ** Action : Update avt, avmu, avmv in statically instable cases 49 51 !! and avt_evd which is avt due to convection 50 !! References : 51 !! Lazar, A., these de l'universite Paris VI, France, 1997 52 !! History : 53 !! 7.0 ! 97-06 (G. Madec, A. Lazar) Original code 54 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 55 !! 9.0 ! 05-06 (C. Ethe) KPP parameterization 52 !! References : Lazar, A., these de l'universite Paris VI, France, 1997 56 53 !!---------------------------------------------------------------------- 57 !! * Arguments58 54 INTEGER, INTENT( in ) :: kt ! ocean time-step indexocean time step 59 60 !! * Local declarations 55 !! 61 56 INTEGER :: ji, jj, jk ! dummy loop indices 62 57 !!---------------------------------------------------------------------- … … 79 74 DO jk = 1, jpkm1 ! Horizontal slab 80 75 ! ! =============== 81 # if defined key_vectopt_loop 82 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 83 jj = 1 ! big loop forced 84 DO ji = jpi+2, jpij 85 # if defined key_zdfkpp 86 !! no implicit mixing in the boundary layer with KPP 87 IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN 88 # else 89 IF( rn2(ji,jj,jk) <= -1.e-12 ) THEN 90 # endif 91 avt (ji ,jj ,jk) = avevd * tmask(ji ,jj ,jk) 92 avmu(ji ,jj ,jk) = avevd * umask(ji ,jj ,jk) 93 avmu(ji-1,jj ,jk) = avevd * umask(ji-1,jj ,jk) 94 avmv(ji ,jj ,jk) = avevd * vmask(ji ,jj ,jk) 95 avmv(ji ,jj-1,jk) = avevd * vmask(ji ,jj-1,jk) 96 ENDIF 97 END DO 98 # else 76 #if defined key_vectopt_loop 77 DO jj = 1, 1 ! big loop forced 78 DO ji = jpi+2, jpij 79 #else 99 80 DO jj = 2, jpj ! no vector opt. 100 81 DO ji = 2, jpi 101 # if defined key_zdfkpp 82 #endif 83 #if defined key_zdfkpp 102 84 !! no implicit mixing in the boundary layer with KPP 103 IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN104 # 105 IF( rn2(ji,jj,jk) <= -1.e-12 ) THEN106 # 85 IF( ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) THEN 86 #else 87 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 88 #endif 107 89 avt (ji ,jj ,jk) = avevd * tmask(ji ,jj ,jk) 108 90 avmu(ji ,jj ,jk) = avevd * umask(ji ,jj ,jk) … … 113 95 END DO 114 96 END DO 115 # endif116 97 ! ! =============== 117 98 END DO ! End of slab … … 129 110 ! ! =============== 130 111 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 131 # if defined key_vectopt_loop 132 jj = 1 ! big loop forced 133 DO ji = 1, jpij 134 # if defined key_zdfkpp 135 !! no implicit mixing in the boundary layer with KPP 136 IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) & 137 avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 138 # else 139 IF( rn2(ji,jj,jk) <= -1.e-12 ) avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 140 # endif 141 END DO 142 # else 112 #if defined key_vectopt_loop 113 DO jj = 1, 1 ! big loop forced 114 DO ji = 1, jpij 115 #else 143 116 DO jj = 1, jpj ! loop over the whole domain (no lbc_lnk call) 144 117 DO ji = 1, jpi 145 # if defined key_zdfkpp 118 #endif 119 #if defined key_zdfkpp 146 120 !! no implicit mixing in the boundary layer with KPP 147 IF( ( rn2(ji,jj,jk) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) &148 avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 149 # else 150 IF( rn2(ji,jj,jk) <= -1.e-12 ) avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 151 # endif 121 IF( ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) .AND. ( fsdepw(ji,jj,jk) > hkpp(ji,jj) ) ) & 122 #else 123 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 124 #endif 125 avt(ji,jj,jk) = avevd * tmask(ji,jj,jk) 152 126 END DO 153 127 END DO 154 # endif155 128 ! ! =============== 156 129 END DO ! End of slab … … 159 132 160 133 ! update of avt_evd and avmu_evd 161 avt_evd (:,:,:) = avt (:,:,:) - avt_evd(:,:,:)162 avmu_evd (:,:,:) = avmu(:,:,:) - avmu_evd(:,:,:)134 avt_evd (:,:,:) = avt (:,:,:) - avt_evd (:,:,:) 135 avmu_evd(:,:,:) = avmu(:,:,:) - avmu_evd(:,:,:) 163 136 164 137 END SUBROUTINE zdf_evd -
trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
r1268 r1438 224 224 !CDIR NOVERRCHK 225 225 DO ji = fs_2, fs_jpim1 ! vector opt. 226 zrn2 = MAX( rn2 (ji,jj,jk), rsmall )226 zrn2 = MAX( rn2b(ji,jj,jk), rsmall ) 227 227 zmxlm(ji,jj,jk) = MAX( rn_lmin, SQRT( 2. * en(ji,jj,jk) / zrn2 ) ) 228 228 END DO … … 321 321 ! 322 322 ! Computation of total energy produce by LC : cumulative sum over jk 323 zpelc(:,:,1) = MAX( rn2 (:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)323 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 324 324 DO jk = 2, jpk 325 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2 (:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)325 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 326 326 END DO 327 327 ! … … 427 427 zdkv = zcoef * ( vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1) & 428 428 & - vb( ji ,jj-1,jk ) - vb(ji,jj,jk ) ) 429 zesh2 = eboost * ( zdku*zdku + zdkv*zdkv ) - rn2 (ji,jj,jk) ! coefficient (zesh2)429 zesh2 = eboost * ( zdku*zdku + zdkv*zdkv ) - rn2b(ji,jj,jk) ! coefficient (zesh2) 430 430 ! 431 431 ! ! Matrix … … 464 464 & - vb(ji ,jj-1,jk ) - vb(ji,jj,jk ) ) 465 465 zsh2 = zdku * zdku + zdkv * zdkv ! square of shear 466 zri = MAX( rn2 (ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 ) ! local Richardson number466 zri = MAX( rn2b(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 ) ! local Richardson number 467 467 # if defined key_c1d 468 468 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri … … 471 471 IF( zri >= 0.2 ) zpdl = 0.2 / zri 472 472 zpdl = MAX( 0.1, zpdl ) 473 zesh2 = eboost * zsh2 - zpdl * rn2 (ji,jj,jk) ! coefficient (esh2)473 zesh2 = eboost * zsh2 - zpdl * rn2b(ji,jj,jk) ! coefficient (esh2) 474 474 ! 475 475 ! ! Matrix -
trunk/NEMO/OPA_SRC/istate.F90
r1200 r1438 132 132 133 133 ENDIF 134 135 IF( lk_ vvl .OR. lk_agrif ) THEN134 ! 135 IF( lk_agrif ) THEN 136 136 ! read free surface arrays in restart file 137 137 IF( ln_rstart ) THEN 138 138 IF( lk_dynspg_flt ) CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields 139 ! ! gcx, gcxb, sshb, sshn 140 IF( lk_dynspg_ts ) CALL ts_rst ( nit000, 'READ' ) ! read or initialize the following fields 141 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 142 IF( lk_dynspg_exp ) CALL exp_rst( nit000, 'READ' ) ! read or initialize the following fields 143 ! ! sshb, sshn 144 ENDIF 139 ! ! gcx, gcxb for agrif_opa_init 140 ENDIF ! explicit case not coded yet with AGRIF 145 141 ENDIF 146 147 IF( lk_vvl ) THEN148 149 !150 IF( .NOT. lk_dynspg_flt ) sshbb(:,:) = sshb(:,:)151 !152 CALL dom_vvl ! ssh init and vertical grid update153 154 CALL eos_init ! usefull to get the equation state type neos parameter155 156 CALL bn2( tb, sb, rn2 ) ! before Brunt Vaissala frequency157 158 IF( .NOT. ln_rstart ) CALL wzv( nit000 )159 160 ENDIF161 162 ! ! Vertical velocity163 ! ! -----------------164 165 IF( .NOT. lk_vvl ) CALL wzv( nit000 ) ! from horizontal divergence166 142 ! 167 143 END SUBROUTINE istate_init -
trunk/NEMO/OPA_SRC/oce.F90
r1152 r1438 4 4 !! Ocean : dynamics and active tracers defined in memory 5 5 !!====================================================================== 6 !! History : 7 !! 8.5 ! 02-11 (G. Madec) F90: Free form and module8 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization6 !! History : 0.1 ! 2002-11 (G. Madec) F90: Free form and module 7 !! 1.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 8 !! 3.1 ! 2009-02 (G. Madec, M. Leclair) pure z* coordinate 9 9 !!---------------------------------------------------------------------- 10 !! OPA 9.0 , LOCEAN-IPSL (2005)11 !! $Id$12 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt13 !!----------------------------------------------------------------------14 !! * Modules used15 10 USE par_oce ! ocean parameters 16 11 … … 20 15 !! Physics and algorithm flags 21 16 !! --------------------------- 22 LOGICAL, PUBLIC :: l_traldf_rot = .FALSE. !: rotated laplacian operator for lateral diffusion 17 LOGICAL, PUBLIC :: l_traldf_rot = .FALSE. !: rotated laplacian operator for lateral diffusion 23 18 LOGICAL, PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 24 19 INTEGER, PUBLIC :: nn_dynhpg_rst = 0 !: add dynhpg implicit variables in restart ot not 25 20 26 !! dynamics and tracer fields 27 !! -------------------------- 28 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: 29 ! before ! now ! after ! ! the after trends becomes the fields 30 ! fields ! fields ! trends ! ! only in dyn(tra)_zdf and dyn(tra)_nxt 31 ub , un , ua , & !: i-horizontal velocity (m/s) 32 vb , vn , va , & !: j-horizontal velocity (m/s) 33 wn , & !: vertical velocity (m/s) 34 rotb , rotn , & !: relative vorticity (1/s) 35 hdivb , hdivn , & !: horizontal divergence (1/s) 36 tb , tn , ta , & !: potential temperature (celcius) 37 sb , sn , sa !: salinity (psu) 38 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: 39 rhd , & !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units) 40 rhop, & !: potential volumic mass (kg/m3) 41 rn2 !: brunt-vaisala frequency (1/s2) 21 !! dynamics and tracer fields ! before ! now ! after ! the after trends becomes the fields 22 !! -------------------------- ! fields ! fields ! trends ! only after tra_zdf and dyn_spg 23 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: ub , un , ua !: i-horizontal velocity [m/s] 24 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: vb , vn , va !: j-horizontal velocity [m/s] 25 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: wn !: vertical velocity [m/s] 26 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: rotb , rotn !: relative vorticity [s-1] 27 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: hdivb , hdivn !: horizontal divergence [s-1] 28 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: tb , tn , ta !: potential temperature [Celcius] 29 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: sb , sn , sa !: salinity [psu] 30 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: rn2b , rn2 !: brunt-vaisala frequency**2 [s-2] 31 ! 32 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: rhd !: in situ density anomalie rhd=(rho-rau0)/rau0 [no units] 33 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: rhop !: potential volumic mass [kg/m3] 42 34 43 !! advection scheme choice 44 !! ----------------------- 45 CHARACTER(len=3), PUBLIC :: l_adv !: 'ce2' centre scheme used 46 ! !: 'tvd' TVD scheme used 47 ! !: 'mus' MUSCL scheme used 48 ! !: 'mu2' MUSCL2 scheme used 35 !! advection scheme choice 36 !! ----------------------- 37 CHARACTER(len=3), PUBLIC :: l_adv !: flag for the advection scheme used (= 'ce2', 'tvd', 'mus' or ...) 49 38 50 39 !! surface pressure gradient 51 40 !! ------------------------- 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 53 spgu, spgv !: horizontal surface pressure gradient 41 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: spgu, spgv !: horizontal surface pressure gradient 54 42 55 43 !! interpolated gradient (only used in zps case) 56 44 !! --------------------- 57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 58 gtu, gsu, gru, & !: t-, s- and rd horizontal gradient at u- and 59 gtv, gsv, grv !: v-points at bottom ocean level 45 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: gtu, gsu, gru !: horizontal gradient of T, S and rd at bottom u-point 46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: gtv, gsv, grv !: horizontal gradient of T, S and rd at bottom v-point 60 47 61 !! free surface 62 !! ------------ 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !:64 sshb , sshn , & !: before, now sea surface height (meters)65 sshu , sshv , & !: sea surface height at u- and v- point66 sshbb, ssha !: before before sea surface height at t-point48 !! free surface ! before ! now ! after ! 49 !! ------------ ! fields ! fields ! trends ! 50 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshb , sshn , ssha !: sea surface height at t-point [m] 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m] 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m] 53 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshf_b , sshf_n , sshf_a !: sea surface height at f-point [m] 67 54 68 55 #if defined key_dynspg_rl || defined key_esopa 69 56 !! rigid-lid formulation 70 57 !! --------------------- 71 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 72 bsfb, bsfn, & !: before, now barotropic streamfunction (m3/s) 73 bsfd !: now trend of barotropic streamfunction (m3/s2) 58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: bsfb, bsfn !: before, now barotropic streamfunction (m3/s) 59 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: bsfd !: now trend of barotropic streamfunction (m3/s2) 74 60 #endif 75 61 !!---------------------------------------------------------------------- 62 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2008) 63 !! $Id$ 64 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 65 !!====================================================================== 76 66 END MODULE oce -
trunk/NEMO/OPA_SRC/restart.F90
r1409 r1438 143 143 CALL iom_rstput( kt, nitrst, numrow, 'rotn' , rotn ) 144 144 CALL iom_rstput( kt, nitrst, numrow, 'hdivn' , hdivn ) 145 CALL iom_rstput( kt, nitrst, numrow, 'rn2' , rn2 ) 145 146 146 147 #if defined key_zdftke2 … … 150 151 CALL iom_rstput( kt, nitrst, numrow, 'rhd' , rhd ) 151 152 #if defined key_zdftke2 152 CALL iom_rstput( kt, nitrst, numrow, 'rn2' , rn2 )153 153 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt ) 154 154 ENDIF … … 239 239 hdivb(:,:,:) = hdivn(:,:,:) 240 240 ENDIF 241 CALL eos_init ! usefull to get the equation state type neos parameter 242 IF( iom_varid( numror, 'rn2', ldstop = .FALSE. ) > 0 ) THEN 243 CALL iom_get( numror, jpdom_autoglo, 'rn2' , rn2 ) 244 ELSE 245 CALL bn2( tb, sb, rn2 ) ! before Brunt-Vaisala frequency 246 ENDIF 241 247 242 248 #if defined key_zdftke2 … … 257 263 #endif 258 264 #if defined key_zdftke2 259 CALL eos_init ! usefull to get the equation state type neos parameter260 IF( iom_varid( numror, 'rn2', ldstop = .FALSE. ) > 0 ) THEN261 CALL iom_get( numror, jpdom_autoglo, 'rn2' , rn2 )262 ELSE263 IF ( ln_dynhpg_imp ) THEN264 CALL bn2( tb, sb, rn2 ) ! before Brunt-Vaisala frequency265 ENDIF266 ENDIF267 265 IF( iom_varid( numror, 'avt', ldstop = .FALSE. ) > 0 ) THEN 268 266 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) -
trunk/NEMO/OPA_SRC/step.F90
r1417 r1438 4 4 !! Time-stepping : manager of the ocean, tracer and ice time stepping 5 5 !!====================================================================== 6 !! History : ! 91-03 (G. Madec) Original code 7 !! ! 92-06 (M. Imbard) add a first output record 8 !! ! 96-04 (G. Madec) introduction of dynspg 9 !! ! 96-04 (M.A. Foujols) introduction of passive tracer 10 !! 8.0 ! 97-06 (G. Madec) new architecture of call 11 !! 8.2 ! 97-06 (G. Madec, M. Imbard, G. Roullet) free surface 12 !! 8.2 ! 99-02 (G. Madec, N. Grima) hpg implicit 13 !! 8.2 ! 00-07 (J-M Molines, M. Imbard) Open Bondary Conditions 14 !! 9.0 ! 02-06 (G. Madec) free form, suppress macro-tasking 15 !! " " ! 04-08 (C. Talandier) New trends organization 16 !! " " ! 05-01 (C. Ethe) Add the KPP closure scheme 17 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 18 !! " " ! 05-11 (G. Madec) Reorganisation of tra and dyn calls 19 !! " " ! 06-01 (L. Debreu, C. Mazauric) Agrif implementation 20 !! " " ! 06-07 (S. Masson) restart using iom 21 !! " " ! 06-08 (G. Madec) surface module 22 !! " " ! 07-07 (J. Chanut, A. Sellar) Unstructured open boundaries (BDY) 6 !! History : OPA ! 1991-03 (G. Madec) Original code 7 !! ! 1991-11 (G. Madec) 8 !! ! 1992-06 (M. Imbard) add a first output record 9 !! ! 1996-04 (G. Madec) introduction of dynspg 10 !! ! 1996-04 (M.A. Foujols) introduction of passive tracer 11 !! 8.0 ! 1997-06 (G. Madec) new architecture of call 12 !! 8.2 ! 1997-06 (G. Madec, M. Imbard, G. Roullet) free surface 13 !! - ! 1999-02 (G. Madec, N. Grima) hpg implicit 14 !! - ! 2000-07 (J-M Molines, M. Imbard) Open Bondary Conditions 15 !! NEMO 1.0 ! 2002-06 (G. Madec) free form, suppress macro-tasking 16 !! - ! 2004-08 (C. Talandier) New trends organization 17 !! ! 2005-01 (C. Ethe) Add the KPP closure scheme 18 !! ! 2005-11 (G. Madec) Reorganisation of tra and dyn calls 19 !! ! 2006-01 (L. Debreu, C. Mazauric) Agrif implementation 20 !! ! 2006-07 (S. Masson) restart using iom 21 !! 3.2 ! 2009-02 (G. Madec, R. Benshila) reintroduicing z*-coordinate 23 22 !!---------------------------------------------------------------------- 24 23 … … 117 116 USE restart ! ocean restart (rst_wri routine) 118 117 USE prtctl ! Print control (prt_ctl routine) 119 USE domvvl ! variable volume (dom_vvl routine)120 118 121 119 #if defined key_agrif … … 141 139 #if defined key_agrif 142 140 SUBROUTINE stp( ) 141 INTEGER :: kstp ! ocean time-step index 143 142 #else 144 143 SUBROUTINE stp( kstp ) 144 INTEGER, INTENT(in) :: kstp ! ocean time-step index 145 145 #endif 146 146 !!---------------------------------------------------------------------- … … 160 160 !! -8- Outputs and diagnostics 161 161 !!---------------------------------------------------------------------- 162 !! * Arguments163 #if defined key_agrif164 INTEGER :: kstp ! ocean time-step index165 #else166 INTEGER, INTENT(in) :: kstp ! ocean time-step index167 #endif168 162 INTEGER :: jk ! dummy loop indice 169 163 INTEGER :: indic ! error indicator if < 0 … … 202 196 ENDIF 203 197 198 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 199 ! Ocean dynamics : ssh, wn, hdiv, rot ! 200 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 201 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity 202 IF( n_cla == 1 ) CALL div_cla( kstp ) ! Cross Land Advection (Update Hor. divergence) 203 CALL ssh_wzv( kstp ) ! after ssh & vertical velocity 204 204 205 205 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 206 206 ! Ocean physics update 207 207 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 208 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 209 !----------------------------------------------------------------------- 210 IF( neuler == 0 .AND. kstp == nit000 ) THEN 211 CALL bn2( tn, sn, rn2 ) ! now Brunt-Vaisala frequency 212 rn2b(:,:,:) = rn2(:,:,:) 213 ELSE 214 rn2b(:,:,:) = rn2(:,:,:) ! before Brunt-Vaisala frequency 215 CALL bn2( tn, sn, rn2 ) ! now Brunt-Vaisala frequency 216 ENDIF 217 208 218 #if defined key_zdftke2 209 219 IF ( ln_dynhpg_imp ) THEN … … 211 221 ! LATERAL PHYSICS 212 222 !----------------------------------------------------------------------- 213 ! N.B. ua, va, ta, sa arrays are used as workspace in this section214 !-----------------------------------------------------------------------215 223 CALL zdf_mxl( kstp ) ! mixed layer depth 216 IF( lk_ldfslp ) CALL ldf_slp( kstp, rhd, rn2 )! before slope of the lateral mixing224 IF( lk_ldfslp ) CALL ldf_slp( kstp, rhd, rn2b ) ! before slope of the lateral mixing 217 225 # if defined key_traldf_c2d 218 226 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient … … 223 231 ! VERTICAL PHYSICS 224 232 !----------------------------------------------------------------------- 225 ! N.B. ua, va, ta, sa arrays are used as workspace in this section226 !-----------------------------------------------------------------------227 #if defined key_zdftke2228 CALL bn2( tn, sn, rn2 ) ! now Brunt-Vaisala frequency229 #else230 CALL bn2( tb, sb, rn2 ) ! before Brunt-Vaisala frequency231 #endif232 233 ! ! Vertical eddy viscosity and diffusivity coefficients 233 234 IF( lk_zdfric ) CALL zdf_ric( kstp ) ! Richardson number dependent Kz … … 261 262 #if defined key_zdftke2 262 263 IF( .NOT. ln_dynhpg_imp ) THEN 263 CALL bn2( tb, sb, rn2 ) ! before Brunt-Vaisala frequency264 264 CALL eos( tb, sb, rhd, rhop ) ! now (swap=before) in situ density for dynhpg module 265 265 #endif … … 267 267 ! LATERAL PHYSICS 268 268 !----------------------------------------------------------------------- 269 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 270 !----------------------------------------------------------------------- 271 IF( lk_ldfslp ) CALL ldf_slp( kstp, rhd, rn2 ) ! before slope of the lateral mixing 269 IF( lk_ldfslp ) CALL ldf_slp( kstp, rhd, rn2b ) ! before slope of the lateral mixing 272 270 #if defined key_traldf_c2d 273 271 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient … … 363 361 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 364 362 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 365 IF( lk_vvl ) CALL dom_vvl ! vertical mesh at next time step 366 367 368 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 369 ! Computation of diagnostic variables 370 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 371 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 372 !----------------------------------------------------------------------- 373 CALL div_cur( kstp ) ! Horizontal divergence & Relative vorticity 374 IF( n_cla == 1 ) CALL div_cla( kstp ) ! Cross Land Advection (Update Hor. divergence) 375 CALL wzv( kstp ) ! Vertical velocity 363 364 CALL ssh_nxt( kstp ) ! sea surface height at next time step 376 365 377 366 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Note: See TracChangeset
for help on using the changeset viewer.