Changeset 1566
- Timestamp:
- 2009-07-31T16:34:08+02:00 (15 years ago)
- Location:
- trunk/NEMO/OPA_SRC
- Files:
-
- 1 deleted
- 13 edited
-
DOM/domcfg.F90 (modified) (9 diffs)
-
DOM/dommsk.F90 (modified) (21 diffs)
-
DOM/domvvl.F90 (modified) (7 diffs)
-
DOM/domzgr.F90 (modified) (6 diffs)
-
DYN/dynadv_cen2.F90 (modified) (9 diffs)
-
DYN/dynadv_ubs.F90 (modified) (14 diffs)
-
DYN/dynnxt.F90 (modified) (1 diff)
-
DYN/dynspg.F90 (modified) (9 diffs)
-
DYN/dynspg_oce.F90 (modified) (1 diff)
-
DYN/wzvmod.F90 (deleted)
-
LDF/ldftra_c3d.h90 (modified) (1 diff)
-
SOL/solmat.F90 (modified) (1 diff)
-
istate.F90 (modified) (1 diff)
-
stpctl.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/domcfg.F90
r1528 r1566 4 4 !! Ocean initialization : domain configuration initialization 5 5 !!============================================================================== 6 !! History : 1.0 ! 2003-09 (G. Madec) Original code 7 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 !!---------------------------------------------------------------------- 6 9 7 10 !!---------------------------------------------------------------------- 8 11 !! dom_cfg : initialize the domain configuration 9 12 !!---------------------------------------------------------------------- 10 !! * Modules used11 13 USE dom_oce ! ocean space and time domain 12 14 USE phycst ! physical constants … … 17 19 PRIVATE 18 20 19 !! * Routine accessibility20 PUBLIC dom_cfg ! called by opa.F90 21 PUBLIC dom_cfg ! called by opa.F90 22 21 23 !!---------------------------------------------------------------------- 22 !! OPA 9.0 , LOCEAN-IPSL (2005)24 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 23 25 !! $Id$ 24 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt26 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 25 27 !!---------------------------------------------------------------------- 26 28 … … 33 35 !! ** Purpose : set the domain configuration 34 36 !! 35 !! ** Method :36 !!37 !! History :38 !! 9.0 ! 03-09 (G. Madec) Original code39 !!----------------------------------------------------------------------40 !! * Local declarations41 37 !!---------------------------------------------------------------------- 42 38 43 IF(lwp) THEN 39 IF(lwp) THEN ! Control print 44 40 WRITE(numout,*) 45 41 WRITE(numout,*) 'dom_cfg : set the ocean configuration' 46 WRITE(numout,*) '~~~~~~~ ocean model configuration used :', & 47 & ' cp_cfg = ', cp_cfg, ' jp_cfg = ', jp_cfg 42 WRITE(numout,*) '~~~~~~~ ' 43 WRITE(numout,*) ' ocean model configuration used : cp_cfg = ', cp_cfg, ' jp_cfg = ', jp_cfg 44 ! 45 WRITE(numout,*) ' global domain lateral boundaries' 46 ! 47 IF( jperio == 0 ) WRITE(numout,*) ' jperio= 0, closed' 48 IF( jperio == 1 ) WRITE(numout,*) ' jperio= 1, cyclic east-west' 49 IF( jperio == 2 ) WRITE(numout,*) ' jperio= 2, equatorial symmetric' 50 IF( jperio == 3 ) WRITE(numout,*) ' jperio= 3, north fold with T-point pivot' 51 IF( jperio == 4 ) WRITE(numout,*) ' jperio= 4, cyclic east-west and north fold with T-point pivot' 52 IF( jperio == 5 ) WRITE(numout,*) ' jperio= 5, north fold with F-point pivot' 53 IF( jperio == 6 ) WRITE(numout,*) ' jperio= 6, cyclic east-west and north fold with F-point pivot' 48 54 ENDIF 49 50 ! Global domain boundary conditions 51 ! --------------------------------- 52 IF(lwp) THEN 53 WRITE(numout,*) ' global domain lateral boundaries' 54 55 IF( jperio == 0 ) WRITE(numout,*) ' jperio= 0, closed' 56 IF( jperio == 1 ) WRITE(numout,*) ' jperio= 1, cyclic east-west' 57 IF( jperio == 2 ) WRITE(numout,*) ' jperio= 2, equatorial symmetric' 58 IF( jperio == 3 ) WRITE(numout,*) ' jperio= 3, north fold with T-point pivot' 59 IF( jperio == 4 ) WRITE(numout,*) ' jperio= 4, cyclic east-west and', & 60 ' north fold with T-point pivot' 61 IF( jperio == 5 ) WRITE(numout,*) ' jperio= 5, north fold with F-point pivot' 62 IF( jperio == 6 ) WRITE(numout,*) ' jperio= 6, cyclic east-west and', & 63 ' north fold with F-point pivot' 64 ENDIF 65 IF( jperio < 0 .OR. jperio > 6 ) CALL ctl_stop( 'jperio is out of range' ) 66 67 ! global domain versus zoom and/or local domain 68 ! --------------------------------------------- 69 70 CALL dom_glo 71 55 ! 56 IF( jperio < 0 .OR. jperio > 6 ) CALL ctl_stop( 'jperio is out of range' ) 57 ! 58 CALL dom_glo ! global domain versus zoom and/or local domain 59 ! 72 60 END SUBROUTINE dom_cfg 73 61 … … 84 72 !! - mi0 , mi1 : 85 73 !! - mj0, , mj1 : 86 !!87 !! History :88 !! 8.5 ! 02-08 (G. Madec) Original code89 74 !!---------------------------------------------------------------------- 90 !! * Local declarations 91 INTEGER :: ji, jj ! dummy loop argument 75 INTEGER :: ji, jj ! dummy loop argument 92 76 !!---------------------------------------------------------------------- 93 77 94 ! Local domain 95 ! ============ 96 97 ! local domain indices ==> data domain indices 98 DO ji = 1, jpi 78 ! ! ============== ! 79 ! ! Local domain ! 80 ! ! ============== ! 81 DO ji = 1, jpi ! local domain indices ==> data domain indices 99 82 mig(ji) = ji + jpizoom - 1 + nimpp - 1 100 83 END DO … … 102 85 mjg(jj) = jj + jpjzoom - 1 + njmpp - 1 103 86 END DO 104 105 ! data domain indices ==> local domain indices106 ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the107 ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.87 ! 88 ! ! data domain indices ==> local domain indices 89 ! ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the 90 ! !local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain. 108 91 DO ji = 1, jpidta 109 92 mi0(ji) = MAX( 1, MIN( ji - jpizoom + 1 - nimpp + 1, jpi+1 ) ) … … 115 98 END DO 116 99 117 IF(lwp) THEN 100 IF(lwp) THEN ! control print 118 101 WRITE(numout,*) 119 102 WRITE(numout,*) 'dom_glo : domain: data / local ' … … 149 132 25 FORMAT( 100(10x,19i4,/) ) 150 133 151 ! Zoom domain152 ! ===========153 154 ! zoom control134 ! ! ============== ! 135 ! ! Zoom domain ! 136 ! ! ============== ! 137 ! ! zoom control 155 138 IF( jpiglo + jpizoom - 1 > jpidta .OR. & 156 139 jpjglo + jpjzoom - 1 > jpjdta ) & 157 140 & CALL ctl_stop( ' global or zoom domain exceed the data domain ! ' ) 158 141 159 ! set zoom flag160 IF 142 ! ! set zoom flag 143 IF( jpiglo < jpidta .OR. jpjglo < jpjdta ) lzoom = .TRUE. 161 144 162 ! set zoom type flags145 ! ! set zoom type flags 163 146 IF( lzoom .AND. jpizoom /= 1 ) lzoom_w = .TRUE. ! 164 147 IF( lzoom .AND. jpjzoom /= 1 ) lzoom_s = .TRUE. … … 180 163 & CALL ctl_stop( ' Your zoom choice is inconsistent with North fold boundary condition' ) 181 164 182 ! Pre-defined arctic/antarctic zoom of ORCA configuration flag165 ! ! Pre-defined arctic/antarctic zoom of ORCA configuration flag 183 166 IF( cp_cfg == "orca" ) THEN 184 167 SELECT CASE ( jp_cfg ) 185 ! ! =======================186 168 CASE ( 2 ) ! ORCA_R2 configuration 187 ! ! =======================188 169 IF( jpiglo == 142 .AND. jpjglo == 53 .AND. & 189 170 & jpizoom == 21 .AND. jpjzoom == 97 ) lzoom_arct = .TRUE. 190 171 IF( jpiglo == jpidta .AND. jpjglo == 50 .AND. & 191 172 & jpizoom == 1 .AND. jpjzoom == 1 ) lzoom_anta = .TRUE. 192 ! ! =======================173 ! 193 174 CASE ( 05 ) ! ORCA_R05 configuration 194 ! ! =======================195 175 IF( jpiglo == 562 .AND. jpjglo == 202 .AND. & 196 176 & jpizoom == 81 .AND. jpjzoom == 301 ) lzoom_arct = .TRUE. … … 204 184 ! 205 185 ENDIF 206 186 ! 207 187 END SUBROUTINE dom_glo 208 188 -
trunk/NEMO/OPA_SRC/DOM/dommsk.F90
r1528 r1566 1 1 MODULE dommsk 2 !!====================================================================== ========2 !!====================================================================== 3 3 !! *** MODULE dommsk *** 4 4 !! Ocean initialization : domain land/sea mask 5 !!============================================================================== 5 !!====================================================================== 6 !! History : OPA ! 1987-07 (G. Madec) Original code 7 !! - ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon) 8 !! - ! 1996-01 (G. Madec) suppression of common work arrays 9 !! - ! 1996-05 (G. Madec) mask computed from tmask and sup- 10 !! ! pression of the double computation of bmask 11 !! - ! 1997-02 (G. Madec) mesh information put in domhgr.F 12 !! - ! 1997-07 (G. Madec) modification of mbathy and fmask 13 !! - ! 1998-05 (G. Roullet) free surface 14 !! - ! 2000-03 (G. Madec) no slip accurate 15 !! - ! 2001-09 (J.-M. Molines) Open boundaries 16 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 17 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 18 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 19 !!---------------------------------------------------------------------- 6 20 7 21 !!---------------------------------------------------------------------- 8 22 !! dom_msk : compute land/ocean mask 9 !! dom_msk_nsa : update land/ocean mask when no-slip accurate 10 !! option is used. 11 !!---------------------------------------------------------------------- 12 !! * Modules used 23 !! dom_msk_nsa : update land/ocean mask when no-slip accurate option is used. 24 !!---------------------------------------------------------------------- 13 25 USE oce ! ocean dynamics and tracers 14 26 USE dom_oce ! ocean space and time domain … … 22 34 PRIVATE 23 35 24 !! * Routine accessibility 25 PUBLIC dom_msk ! routine called by inidom.F90 26 27 !! * Module variables 28 REAL(wp) :: & 29 shlat = 2. ! type of lateral boundary condition on velocity (namelist namlbc) 36 PUBLIC dom_msk ! routine called by inidom.F90 37 38 REAL(wp) :: shlat = 2. ! type of lateral boundary condition on velocity (namelist namlbc) 30 39 31 40 !! * Substitutions 32 41 # include "vectopt_loop_substitute.h90" 33 !!---------------------------------------------------------------------- -----------34 !! OPA 9.0 , LOCEAN-IPSL (2005)42 !!---------------------------------------------------------------------- 43 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 35 44 !! $Id$ 36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt37 !!---------------------------------------------------------------------- -----------45 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 46 !!---------------------------------------------------------------------- 38 47 39 48 CONTAINS … … 93 102 !! (note that the minimum value of mbathy is 2). 94 103 !! 95 !! ** Action : 96 !! tmask : land/ocean mask at t-point (=0. or 1.) 97 !! umask : land/ocean mask at u-point (=0. or 1.) 98 !! vmask : land/ocean mask at v-point (=0. or 1.) 99 !! fmask : land/ocean mask at f-point (=0. or 1.) 104 !! ** Action : tmask : land/ocean mask at t-point (=0. or 1.) 105 !! umask : land/ocean mask at u-point (=0. or 1.) 106 !! vmask : land/ocean mask at v-point (=0. or 1.) 107 !! fmask : land/ocean mask at f-point (=0. or 1.) 100 108 !! =shlat along lateral boundaries 101 !! bmask : land/ocean mask at barotropic stream 102 !! function point (=0. or 1.) and set to 103 !! 0 along lateral boundaries 104 !! mbathy : number of non-zero w-levels 105 !! 106 !! History : 107 !! ! 87-07 (G. Madec) Original code 108 !! ! 91-12 (G. Madec) 109 !! ! 92-06 (M. Imbard) 110 !! ! 93-03 (M. Guyon) symetrical conditions (M. Guyon) 111 !! ! 96-01 (G. Madec) suppression of common work arrays 112 !! ! 96-05 (G. Madec) mask computed from tmask and sup- 113 !! pression of the double computation of bmask 114 !! ! 97-02 (G. Madec) mesh information put in domhgr.F 115 !! ! 97-07 (G. Madec) modification of mbathy and fmask 116 !! ! 98-05 (G. Roullet) free surface 117 !! ! 00-03 (G. Madec) no slip accurate 118 !! ! 01-09 (J.-M. Molines) Open boundaries 119 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 120 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 109 !! bmask : land/ocean mask at barotropic stream 110 !! function point (=0. or 1.) and set to 0 along lateral boundaries 111 !! mbathy : number of non-zero w-levels 121 112 !!---------------------------------------------------------------------- 122 113 INTEGER :: ji, jj, jk ! dummy loop indices … … 129 120 !!--------------------------------------------------------------------- 130 121 131 ! Namelist namlbc : lateral momentum boundary condition 132 REWIND( numnam ) 122 REWIND( numnam ) ! Namelist namlbc : lateral momentum boundary condition 133 123 READ ( numnam, namlbc ) 134 IF(lwp) THEN 124 125 IF(lwp) THEN ! control print 135 126 WRITE(numout,*) 136 127 WRITE(numout,*) 'dommsk : ocean mask ' 137 128 WRITE(numout,*) '~~~~~~' 138 WRITE(numout,*) ' Namelist namlbc' 139 WRITE(numout,*) ' lateral momentum boundary cond. shlat = ',shlat 140 ENDIF 141 142 IF ( shlat == 0. ) THEN 143 IF(lwp) WRITE(numout,*) ' ocean lateral free-slip ' 144 ELSEIF ( shlat == 2. ) THEN 145 IF(lwp) WRITE(numout,*) ' ocean lateral no-slip ' 146 ELSEIF ( 0. < shlat .AND. shlat < 2. ) THEN 147 IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip ' 148 ELSEIF ( 2. < shlat ) THEN 149 IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip ' 129 WRITE(numout,*) ' Namelist namlbc' 130 WRITE(numout,*) ' lateral momentum boundary cond. shlat = ',shlat 131 ENDIF 132 133 IF( shlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral free-slip ' 134 ELSEIF ( shlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral no-slip ' 135 ELSEIF ( 0. < shlat .AND. shlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral partial-slip ' 136 ELSEIF ( 2. < shlat ) THEN ; IF(lwp) WRITE(numout,*) ' ocean lateral strong-slip ' 150 137 ELSE 151 138 WRITE(ctmp1,*) ' shlat is negative = ', shlat … … 155 142 ! 1. Ocean/land mask at t-point (computed from mbathy) 156 143 ! ----------------------------- 157 ! Tmask has already the right boundary conditions since mbathy is ok158 144 ! N.B. tmask has already the right boundary conditions since mbathy is ok 145 ! 159 146 tmask(:,:,:) = 0.e0 160 147 DO jk = 1, jpk 161 148 DO jj = 1, jpj 162 149 DO ji = 1, jpi 163 IF( FLOAT( mbathy(ji,jj)-jk )+.1 >= 0.e0 ) tmask(ji,jj,jk) = 1.e0150 IF( REAL( mbathy(ji,jj) - jk ) +.1 >= 0.e0 ) tmask(ji,jj,jk) = 1.e0 164 151 END DO 165 152 END DO 166 153 END DO 167 154 155 !!gm ???? 168 156 #if defined key_zdfkpp 169 157 IF( cp_cfg == 'orca' ) THEN … … 184 172 ENDIF 185 173 #endif 174 !!gm end 186 175 187 176 ! Interior domain mask (used for global sum) 188 177 ! -------------------- 189 190 178 tmask_i(:,:) = tmask(:,:,1) 191 179 iif = jpreci ! ??? … … 200 188 201 189 ! north fold mask 190 ! --------------- 202 191 tpol(1:jpiglo) = 1.e0 203 192 fpol(1:jpiglo) = 1.e0 … … 205 194 tpol(jpiglo/2+1:jpiglo) = 0.e0 206 195 fpol( 1 :jpiglo) = 0.e0 207 ! T-point pivot: only half of the nlcj-1 row 208 IF( mjg(nlej) == jpjglo ) THEN 196 IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row 209 197 DO ji = iif+1, iil-1 210 198 tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji)) … … 219 207 ! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask) 220 208 ! ------------------------------------------- 221 222 ! Computation223 209 DO jk = 1, jpk 224 210 DO jj = 1, jpjm1 … … 233 219 END DO 234 220 END DO 235 236 ! Lateral boundary conditions 237 CALL lbc_lnk( umask, 'U', 1. ) 221 CALL lbc_lnk( umask, 'U', 1. ) ! Lateral boundary conditions 238 222 CALL lbc_lnk( vmask, 'V', 1. ) 239 223 CALL lbc_lnk( fmask, 'F', 1. ) … … 242 226 ! 4. ocean/land mask for the elliptic equation 243 227 ! -------------------------------------------- 244 245 ! Computation246 228 bmask(:,:) = tmask(:,:,1) ! elliptic equation is written at t-point 247 248 ! Boundary conditions249 ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi229 ! 230 ! ! Boundary conditions 231 ! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi 250 232 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 251 233 bmask( 1 ,:) = 0.e0 252 234 bmask(jpi,:) = 0.e0 253 235 ENDIF 254 255 ! south symmetric : bmask must be set to 0. on row 1 256 IF( nperio == 2 ) THEN 236 IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1 257 237 bmask(:, 1 ) = 0.e0 258 238 ENDIF 259 260 ! north fold : 261 IF( nperio == 3 .OR. nperio == 4 ) THEN 262 ! T-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj and on half jpjglo-1 row 263 DO ji = 1, jpi 239 ! ! north fold : 240 IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row 241 DO ji = 1, jpi 264 242 ii = ji + nimpp - 1 265 243 bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii) … … 267 245 END DO 268 246 ENDIF 269 IF( nperio == 5 .OR. nperio == 6 ) THEN 270 ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj 247 IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj 271 248 bmask(:,jpj) = 0.e0 272 249 ENDIF 273 274 ! Mpp boundary conditions: bmask is set to zero on the overlap 275 ! region for all elliptic solvers 276 277 IF( lk_mpp ) THEN 250 ! 251 IF( lk_mpp ) THEN ! mpp specificities 252 ! ! bmask is set to zero on the overlap region 278 253 IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0.e0 279 254 IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0.e0 280 255 IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0.e0 281 256 IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0.e0 282 283 ! north fold : bmask must be set to 0. on rows jpj-1 and jpj 284 IF( npolj == 3 .OR. npolj == 4 ) THEN 257 ! 258 IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj 285 259 DO ji = 1, nlci 286 260 ii = ji + nimpp - 1 … … 289 263 END DO 290 264 ENDIF 291 IF( npolj == 5 .OR. npolj == 6 ) THEN 265 IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj 292 266 DO ji = 1, nlci 293 267 bmask(ji,nlcj ) = 0.e0 … … 299 273 ! mask for second order calculation of vorticity 300 274 ! ---------------------------------------------- 301 302 275 CALL dom_msk_nsa 303 276 304 277 305 278 ! Lateral boundary conditions on velocity (modify fmask) 306 ! --------------------------------------- 307 279 ! --------------------------------------- 308 280 DO jk = 1, jpk 309 310 zwf(:,:) = fmask(:,:,jk) 311 281 zwf(:,:) = fmask(:,:,jk) 312 282 DO jj = 2, jpjm1 313 283 DO ji = fs_2, fs_jpim1 ! vector opt. … … 318 288 END DO 319 289 END DO 320 321 290 DO jj = 2, jpjm1 322 291 IF( fmask(1,jj,jk) == 0. ) THEN … … 326 295 fmask(jpi,jj,jk) = shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 327 296 ENDIF 328 END DO 329 297 END DO 330 298 DO ji = 2, jpim1 331 299 IF( fmask(ji,1,jk) == 0. ) THEN … … 337 305 END DO 338 306 END DO 339 340 341 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN 342 ! ! ======================= 343 ! Increased lateral friction in ! ORCA_R2 configuration 344 ! the vicinity of some straits ! ======================= 345 ! 307 ! 308 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration 309 ! ! Increased lateral friction near of some straits 346 310 IF( n_cla == 0 ) THEN 347 311 ! ! Gibraltar strait : partial slip (fmask=0.5) … … 365 329 ! 366 330 ENDIF 367 368 ! Lateral boundary conditions on fmask369 CALL lbc_lnk( fmask, 'F', 1. ) 331 ! 332 CALL lbc_lnk( fmask, 'F', 1. ) ! Lateral boundary conditions on fmask 333 370 334 371 335 ! Mbathy set to the number of w-level (minimum value 2) … … 377 341 END DO 378 342 379 ! Control print 380 ! ------------- 381 IF( nprint == 1 .AND. lwp ) THEN 343 IF( nprint == 1 .AND. lwp ) THEN ! Control print 382 344 imsk(:,:) = INT( tmask_i(:,:) ) 383 345 WRITE(numout,*) ' tmask_i : ' … … 423 385 & 1, jpj, 1, 1, numout ) 424 386 ENDIF 425 387 ! 426 388 END SUBROUTINE dom_msk 427 389 … … 441 403 !! ** Action : 442 404 !! 443 !! History :444 !! ! 00-03 (G. Madec) no slip accurate445 405 !!---------------------------------------------------------------------- 446 406 INTEGER :: ji, jj, jk, jl ! dummy loop indices 447 INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd 407 INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd 408 REAL(wp) :: zaa 448 409 INTEGER, DIMENSION(jpi*jpj*jpk,3) :: icoord 449 REAL(wp) :: zaa450 410 !!--------------------------------------------------------------------- 451 411 -
trunk/NEMO/OPA_SRC/DOM/domvvl.F90
r1528 r1566 24 24 PRIVATE 25 25 26 PUBLIC dom_vvl! called by domain.F9026 PUBLIC dom_vvl ! called by domain.F90 27 27 28 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ee_t, ee_u, ee_v, ee_f !: ??? 28 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ee_t, ee_u, ee_v, ee_f !: ??? 29 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: mut, muu, muv, muf !: ??? 29 30 30 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: mut, muu, muv, muf !: ??? 31 32 REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time-step, = 2 rdttra 33 ! ! except at nit000 (=rdttra) if neuler=0 31 REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time-step, = 2 rdttra 32 ! ! except at nit000 (=rdttra) if neuler=0 34 33 35 34 !! * Substitutions … … 50 49 !! ** Purpose : compute coefficients muX at T-U-V-F points to spread 51 50 !! ssh over the whole water column (scale factors) 52 !!53 51 !!---------------------------------------------------------------------- 54 52 INTEGER :: ji, jj, jk … … 62 60 ENDIF 63 61 64 #if defined key_zco 65 CALL ctl_stop( 'dom_vvl_ini : options key_zco is incompatible with variable volume option key_vvl') 66 #endif 62 IF( lk_zco ) CALL ctl_stop( 'dom_vvl : key_zco is incompatible with variable volume option key_vvl') 67 63 68 64 fsdept(:,:,:) = gdept (:,:,:) … … 77 73 fse3vw(:,:,:) = e3vw (:,:,:) 78 74 79 ! mu computation 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 75 ! !== mu computation ==! 76 ee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 83 77 ee_u(:,:) = fse3u_0(:,:,1) 84 78 ee_v(:,:) = fse3v_0(:,:,1) 85 79 ee_f(:,:) = fse3f_0(:,:,1) 86 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors80 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 87 81 ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 88 82 ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) … … 92 86 END DO 93 87 END DO 94 ! ! Compute and mask the inverse of the local depth at T, U, V and F points88 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 95 89 ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1) 96 90 ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1) 97 91 ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1) 98 DO jj = 1, jpjm1 ! f-point case fmask cannot be used92 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 99 93 ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 100 94 END DO 101 CALL lbc_lnk( ee_f, 'F', 1. ) ! lateral boundary condition on ee_f95 CALL lbc_lnk( ee_f, 'F', 1. ) ! lateral boundary condition on ee_f 102 96 ! 103 DO jk = 1, jpk 104 mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk) !at T levels105 muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk) !at T levels106 muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk) !at T levels97 DO jk = 1, jpk ! mu coefficients 98 mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk) ! T-point at T levels 99 muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk) ! U-point at T levels 100 muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk) ! V-point at T levels 107 101 END DO 108 DO jk = 1, jpk 109 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask102 DO jk = 1, jpk ! F-point : fmask=shlat at coasts, use the product of umask 103 DO jj = 1, jpjm1 110 104 muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 111 105 END DO 112 106 muf(:,jpj,jk) = 0.e0 113 107 END DO 114 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition on ee_f108 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition 115 109 116 110 117 ! Reference ocean depth at U- and V-points 118 hu_0(:,:) = 0.e0 111 hu_0(:,:) = 0.e0 ! Reference ocean depth at U- and V-points 119 112 hv_0(:,:) = 0.e0 120 113 DO jk = 1, jpk … … 123 116 END DO 124 117 125 ! before and now Sea Surface Height at u-, v-, f-points 126 DO jj = 1, jpjm1 118 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 127 119 DO ji = 1, jpim1 128 120 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) … … 143 135 END DO 144 136 END DO 145 ! Boundaries conditions 146 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 137 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) ! lateral boundary conditions 147 138 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 148 139 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) -
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r1528 r1566 4 4 !! Ocean initialization : domain initialization 5 5 !!============================================================================== 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate 7 !! ! 1997-07 (G. Madec) lbc_lnk call 8 !! ! 1997-04 (J.-O. Beismann) 9 !! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module 10 !! - ! 2002-09 (A. de Miranda) rigid-lid + islands 11 !! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module 12 !! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function 13 !! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco 14 !! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate 7 !! ! 1997-07 (G. Madec) lbc_lnk call 8 !! ! 1997-04 (J.-O. Beismann) 9 !! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module 10 !! - ! 2002-09 (A. de Miranda) rigid-lid + islands 11 !! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module 12 !! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function 13 !! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco 14 !! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style 15 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 15 16 !!---------------------------------------------------------------------- 16 17 … … 623 624 ENDIF 624 625 625 ! Set to zero mbathy over islands if necessary626 IF(lwp) WRITE(numout,*)627 IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands'628 IF(lwp) WRITE(numout,*) ' ----------------------------'629 !630 mbathy(:,:) = MAX( 0, mbathy(:,:) )631 !632 626 ! Boundary condition on mbathy 633 627 IF( .NOT.lk_mpp ) THEN … … 655 649 ENDIF 656 650 657 ! control print 658 IF( lwp .AND. nprint == 1 ) THEN 651 IF( lwp .AND. nprint == 1 ) THEN ! control print 659 652 WRITE(numout,*) 660 653 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels ' … … 1027 1020 REAL(wp), INTENT(in ) :: bb ! Stretching coefficient 1028 1021 REAL(wp) :: pf1 ! sigma value 1029 1030 !!---------------------------------------------------------------------- 1031 ! 1032 IF ( theta == 0 ) then !uniform sigma 1033 pf1 = -(pk1-0.5)/REAL(jpkm1) 1034 ELSE ! stretched sigma 1022 !!---------------------------------------------------------------------- 1023 ! 1024 IF ( theta == 0 ) then ! uniform sigma 1025 pf1 = -(pk1-0.5) / REAL( jpkm1 ) 1026 ELSE ! stretched sigma 1035 1027 pf1 = (1.0-bb) * (sinh( theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(theta) + & 1036 1037 1038 ENDIF 1039 1028 & bb * ( (tanh( theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*theta) ) / & 1029 & (2*tanh(0.5*theta) ) ) 1030 ENDIF 1031 ! 1040 1032 END FUNCTION fssig1 1041 1033 … … 1078 1070 REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace 1079 1071 REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - - 1080 1081 LOGICAL :: ln_s_sigma = .false. !use hybrid s_sigma coordinates & stretching function fssig1,used with ln_sco = .true.1072 !! 1073 LOGICAL :: ln_s_sigma = .false. !use hybrid s_sigma coordinates & stretching function fssig1,used with ln_sco = .true. 1082 1074 REAL(wp) :: bb = 0.8 ! stretching parameter for song and haidvogel stretching, bb=0; top only, bb =1; top and bottom 1083 1075 REAL(wp) :: hc = 150 ! Critical depth for s-sigma coordinates 1084 1076 !!gm never do that !!!! ==> Pb at compilation phase on several computer 1085 1077 REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3 = 0.0d0 1086 1078 REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigt3 = 0.0d0 … … 1093 1085 REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwu3 = 0.0d0 1094 1086 REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwv3 = 0.0d0 1087 !!gm end 1095 1088 !! 1096 1089 NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max, ln_s_sigma, bb, hc -
trunk/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r1528 r1566 5 5 !! using a 2nd order centred scheme 6 6 !!====================================================================== 7 !! History : 9.0 ! 06-08 (G. Madec, S. Theetten) Original code 7 !! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code 8 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 9 !!---------------------------------------------------------------------- 9 10 … … 14 15 USE oce ! ocean dynamics and tracers 15 16 USE dom_oce ! ocean space and time domain 17 USE trdmod_oce ! ocean variables trends 18 USE trdmod ! ocean dynamics trends 16 19 USE in_out_manager ! I/O manager 17 USE trdmod ! ocean dynamics trends18 USE trdmod_oce ! ocean variables trends19 20 USE prtctl ! Print control 20 21 … … 22 23 PRIVATE 23 24 24 !! * Routine accessibility 25 PUBLIC dyn_adv_cen2 ! routine called by step.F90 25 PUBLIC dyn_adv_cen2 ! routine called by step.F90 26 26 27 27 !! * Substitutions … … 29 29 # include "vectopt_loop_substitute.h90" 30 30 !!---------------------------------------------------------------------- 31 !! OPA 9.0 , LODYC-IPSL (2006)31 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 32 32 !! $Id$ 33 33 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 41 41 !! 42 42 !! ** Purpose : Compute the now momentum advection trend in flux form 43 !! and the general trend of the momentum equation.43 !! and the general trend of the momentum equation. 44 44 !! 45 45 !! ** Method : Trend evaluated using now fields (centered in time) 46 46 !! 47 !! ** Action : - Update (ua,va)with the now vorticity term trend47 !! ** Action : (ua,va) updated with the now vorticity term trend 48 48 !!---------------------------------------------------------------------- 49 USE oce, ONLY: zfu => ta , &! use ta as 3D workspace50 zfv => sa! use sa as 3D workspace51 52 INTEGER, INTENT( in ) :: kt 53 54 INTEGER :: ji, jj, jk 55 REAL(wp) :: z ua, zva, zbu, zbv! temporary scalars49 USE oce, ONLY: zfu => ta ! use ta as 3D workspace 50 USE oce, ONLY: zfv => sa ! use sa as 3D workspace 51 !! 52 INTEGER, INTENT( in ) :: kt ! ocean time-step index 53 !! 54 INTEGER :: ji, jj, jk ! dummy loop indices 55 REAL(wp) :: zbu, zbv ! temporary scalars 56 56 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu_t, zfu_f, zfu_uw ! 3D workspace 57 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfv_t, zfv_f, zfv_vw ! " "58 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfw ! " "57 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfv_t, zfv_f, zfv_vw ! - - 58 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfw ! - - 59 59 !!---------------------------------------------------------------------- 60 60 … … 70 70 ENDIF 71 71 72 ! I. Horizontal advection 73 ! ----------------------- 74 ! ! =============== 75 DO jk = 1, jpkm1 ! Horizontal slab 76 ! ! =============== 77 ! horizontal volume fluxes 72 ! ! ====================== ! 73 ! ! Horizontal advection ! 74 DO jk = 1, jpkm1 ! ====================== ! 75 ! ! horizontal volume fluxes 78 76 zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 79 77 zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 80 81 ! horizontal momentum fluxes at T- and F-point 82 DO jj = 1, jpjm1 78 ! 79 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point 83 80 DO ji = 1, fs_jpim1 ! vector opt. 84 81 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) … … 88 85 END DO 89 86 END DO 90 91 ! divergence of horizontal momentum fluxes 92 DO jj = 2, jpjm1 87 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 93 88 DO ji = fs_2, fs_jpim1 ! vector opt. 94 89 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 95 90 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 96 ! horizontal advective trends 97 zua = - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 98 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 99 zva = - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 100 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 101 ! add it to the general tracer trends 102 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 103 va(ji,jj,jk) = va(ji,jj,jk) + zva 91 ! 92 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 93 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 94 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 95 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 104 96 END DO 105 97 END DO 106 ! ! =============== 107 END DO ! End of slab 108 ! ! =============== 109 110 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 98 END DO 99 ! 100 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 111 101 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 112 102 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 113 103 CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt ) 104 zfu_t(:,:,:) = ua(:,:,:) 105 zfv_t(:,:,:) = va(:,:,:) 114 106 ENDIF 115 107 ! 116 108 117 ! II. Vertical advection 118 ! ---------------------- 119 120 IF( l_trddyn ) THEN ! Save ua and va trends 121 zfu_t(:,:,:) = ua(:,:,:) 122 zfv_t(:,:,:) = va(:,:,:) 123 ENDIF 124 125 ! Second order centered tracer flux at w-point 126 DO jk = 1, jpkm1 127 ! Vertical volume fluxes 109 ! ! ==================== ! 110 ! ! Vertical advection ! 111 DO jk = 1, jpkm1 ! ==================== ! 112 ! ! Vertical volume fluxesÊ 128 113 zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk) 129 ! Vertical advective fluxes130 IF( jk == 1 ) THEN 131 zfu_uw(:,:,jpk) = 0.e0 ! Bottomvalue : flux set to zero114 ! 115 IF( jk == 1 ) THEN ! surface/bottom advective fluxes 116 zfu_uw(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 132 117 zfv_vw(:,:,jpk) = 0.e0 133 ! ! Surface value 134 IF( lk_vvl ) THEN 135 ! variable volume : flux set to zero 118 ! ! Surface value : 119 IF( lk_vvl ) THEN ! variable volume : flux set to zero 136 120 zfu_uw(:,:, 1 ) = 0.e0 137 121 zfv_vw(:,:, 1 ) = 0.e0 138 ELSE 139 ! free surface-constant volume 122 ELSE ! constant volume : advection through the surface 140 123 DO jj = 2, jpjm1 141 124 DO ji = fs_2, fs_jpim1 … … 145 128 END DO 146 129 ENDIF 147 ELSE 148 ! ! interior fluxes 130 ELSE ! interior fluxes 149 131 DO jj = 2, jpjm1 150 132 DO ji = fs_2, fs_jpim1 ! vector opt. … … 155 137 ENDIF 156 138 END DO 157 158 ! momentum flux divergence at u-, v-points added to the general trend 159 DO jk = 1, jpkm1 139 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence 160 140 DO jj = 2, jpjm1 161 141 DO ji = fs_2, fs_jpim1 ! vector opt. 162 zua =- ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) &142 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 163 143 & / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 164 zva =- ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) &144 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 165 145 & / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 166 ! add it to the general tracer trends167 ua(ji,jj,jk) = ua(ji,jj,jk) + zua168 va(ji,jj,jk) = va(ji,jj,jk) + zva169 146 END DO 170 147 END DO 171 148 END DO 172 173 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic149 ! 150 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 174 151 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 175 152 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 176 153 CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt ) 177 154 ENDIF 178 179 ! ! Control print 155 ! ! Control print 180 156 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask, & 181 157 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) -
trunk/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r1528 r1566 5 5 !! trend using a 3rd order upstream biased scheme 6 6 !!====================================================================== 7 !! History : 9.0 ! 06-08 (R. Benshila, L. Debreu) Original code 7 !! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code 8 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 9 !!---------------------------------------------------------------------- 9 10 … … 27 28 REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp ! =0 2nd order ; =1/8 4th order centred 28 29 29 !! * Routine accessibility 30 PUBLIC dyn_adv_ubs ! routine called by step.F90 30 PUBLIC dyn_adv_ubs ! routine called by step.F90 31 31 32 32 !! * Substitutions … … 34 34 # include "vectopt_loop_substitute.h90" 35 35 !!---------------------------------------------------------------------- 36 !! OPA 9.0 , LODYC-IPSL (2006)36 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 37 37 !! $Id$ 38 38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 46 46 !! 47 47 !! ** Purpose : Compute the now momentum advection trend in flux form 48 !! and the general trend of the momentum equation.48 !! and the general trend of the momentum equation. 49 49 !! 50 50 !! ** Method : The scheme is the one implemeted in ROMS. It depends … … 64 64 !! gamma1=1/4 and gamma2=1/8. 65 65 !! 66 !! ** Action : - update (ua,va)with the 3D advective momentum trends66 !! ** Action : - (ua,va) updated with the 3D advective momentum trends 67 67 !! 68 68 !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 69 69 !!---------------------------------------------------------------------- 70 USE oce, ONLY: zfu => ta , &! use ta as 3D workspace71 zfv => sa! use sa as 3D workspace72 70 USE oce, ONLY: zfu => ta ! use ta as 3D workspace 71 USE oce, ONLY: zfv => sa ! use sa as 3D workspace 72 !! 73 73 INTEGER, INTENT(in) :: kt ! ocean time-step index 74 75 INTEGER :: ji, jj, jk 76 REAL(wp) :: z ua, zva, zbu, zbv! temporary scalars77 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! temporary scalars74 !! 75 INTEGER :: ji, jj, jk ! dummy loop indices 76 REAL(wp) :: zbu, zbv ! temporary scalars 77 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! temporary scalars 78 78 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu_t, zfu_f ! temporary workspace 79 79 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfv_t, zfv_f ! " " … … 92 92 zfu_f(:,:,:) = 0.e0 93 93 zfv_f(:,:,:) = 0.e0 94 94 ! 95 95 zlu_uu(:,:,:,:) = 0.e0 96 96 zlv_vv(:,:,:,:) = 0.e0 … … 103 103 ENDIF 104 104 105 ! ! =============== 106 DO jk = 1, jpkm1 ! Horizontal slab 107 ! ! =============== 108 ! Laplacian 109 ! --------- 105 ! ! =========================== ! 106 DO jk = 1, jpkm1 ! Laplacian of the velocity ! 107 ! ! =========================== ! 108 ! ! horizontal volume fluxes 110 109 zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 111 110 zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 112 113 DO jj = 2, jpjm1 111 ! 112 DO jj = 2, jpjm1 ! laplacian 114 113 DO ji = fs_2, fs_jpim1 ! vector opt. 115 114 zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj,jk)-2.*ub (ji,jj,jk)+ub (ji-1,jj,jk) ) * umask(ji,jj,jk) … … 124 123 END DO 125 124 END DO 126 ! ! =============== 127 END DO ! End of slab 128 ! ! =============== 129 130 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) 131 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) 132 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) 133 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) 134 135 CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 136 CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 137 CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 138 CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 139 140 ! I. Horizontal advection 141 ! ----------------------- 142 ! ! =============== 143 DO jk = 1, jpkm1 ! Horizontal slab 144 ! ! =============== 145 ! horizontal volume fluxes 125 END DO 126 !!gm BUG !!! just below this should be +1 in all the communications 127 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 128 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 129 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 130 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 131 132 !!gm corrected: 133 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. ) 134 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. ) 135 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 136 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 137 !!gm end 138 139 ! ! ====================== ! 140 ! ! Horizontal advection ! 141 DO jk = 1, jpkm1 ! ====================== ! 142 ! ! horizontal volume fluxes 146 143 zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 147 144 zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 148 149 ! horizontal momentum fluxes at T- and F-point 150 DO jj = 1, jpjm1 145 ! 146 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point 151 147 DO ji = 1, fs_jpim1 ! vector opt. 152 148 zui = ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) 153 149 zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 154 150 ! 155 151 IF (zui > 0) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) 156 152 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) … … 159 155 ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1) 160 156 ENDIF 161 157 ! 162 158 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) & 163 159 & - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj ,jk,2) ) ) & … … 166 162 & - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji ,jj+1,jk,2) ) ) & 167 163 & * ( zvj - gamma1 * zl_v) 168 164 ! 169 165 zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) ) 170 166 zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) ) … … 175 171 ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1) 176 172 ENDIF 177 173 ! 178 174 zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) & 179 175 & * ( un(ji,jj,jk) + un(ji ,jj+1,jk) - gamma1 * zl_u ) … … 182 178 END DO 183 179 END DO 184 185 ! divergence of horizontal momentum fluxes 186 DO jj = 2, jpjm1 180 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 187 181 DO ji = fs_2, fs_jpim1 ! vector opt. 188 182 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 189 183 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 190 ! horizontal advective trends 191 zua = - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 192 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 193 zva = - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 194 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 195 ! add it to the general tracer trends 196 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 197 va(ji,jj,jk) = va(ji,jj,jk) + zva 198 END DO 199 END DO 200 ! ! =============== 201 END DO ! End of slab 202 ! ! =============== 203 204 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 184 ! 185 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 186 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 187 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 188 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 189 END DO 190 END DO 191 END DO 192 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 205 193 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 206 194 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 207 195 CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt ) 208 ENDIF209 210 ! II. Vertical advection211 ! ----------------------212 213 IF( l_trddyn ) THEN ! Save ua and va trends214 196 zfu_t(:,:,:) = ua(:,:,:) 215 197 zfv_t(:,:,:) = va(:,:,:) 216 198 ENDIF 217 199 218 ! Second order centered tracer flux at w-point 219 DO jk = 1, jpkm1 220 ! Vertical volume fluxes 200 ! ! ==================== ! 201 ! ! Vertical advection ! 202 DO jk = 1, jpkm1 ! ==================== ! 203 ! ! Vertical volume fluxesÊ 221 204 zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk) 222 ! Vertical advective fluxes223 IF( jk == 1 ) THEN 224 zfu_uw(:,:,jpk) = 0.e0 ! Bottomvalue : flux set to zero205 ! 206 IF( jk == 1 ) THEN ! surface/bottom advective fluxes 207 zfu_uw(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 225 208 zfv_vw(:,:,jpk) = 0.e0 226 ! ! Surface value 227 IF( lk_vvl ) THEN 228 ! variable volume : flux set to zero 209 ! ! Surface value : 210 IF( lk_vvl ) THEN ! variable volume : flux set to zero 229 211 zfu_uw(:,:, 1 ) = 0.e0 230 212 zfv_vw(:,:, 1 ) = 0.e0 231 ELSE 232 ! free surface-constant volume 213 ELSE ! constant volume : advection through the surface 233 214 DO jj = 2, jpjm1 234 215 DO ji = fs_2, fs_jpim1 … … 238 219 END DO 239 220 ENDIF 240 ELSE 241 ! ! interior fluxes 221 ELSE ! interior fluxes 242 222 DO jj = 2, jpjm1 243 223 DO ji = fs_2, fs_jpim1 ! vector opt. … … 248 228 ENDIF 249 229 END DO 250 251 ! momentum flux divergence at u-, v-points added to the general trend 252 DO jk = 1, jpkm1 230 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence 253 231 DO jj = 2, jpjm1 254 232 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zua =- ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) &233 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 256 234 & / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 257 zva =- ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) &235 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 258 236 & / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 259 ! add it to the general tracer trends 260 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 261 va(ji,jj,jk) = va(ji,jj,jk) + zva 262 END DO 263 END DO 264 END DO 265 266 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 237 END DO 238 END DO 239 END DO 240 ! 241 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 267 242 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 268 243 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 269 244 CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt ) 270 245 ENDIF 271 272 ! ! Control print 246 ! ! Control print 273 247 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask, & 274 248 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 275 249 ! 276 250 END SUBROUTINE dyn_adv_ubs 277 251 -
trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
r1502 r1566 90 90 !! 91 91 INTEGER :: ji, jj, jk ! dummy loop indices 92 #if ! defined key_dynspg_flt 92 93 REAL(wp) :: z2dt ! temporary scalar 94 #endif 93 95 REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar 94 96 REAL(wp) :: zve3a , zve3n , zve3b ! - - -
trunk/NEMO/OPA_SRC/DYN/dynspg.F90
r1528 r1566 4 4 !! Ocean dynamics: surface pressure gradient control 5 5 !!====================================================================== 6 !! History : 9.0 ! 05-12 (C. Talandier, G. Madec) Original code7 !! 9.0 ! 05-12 (V. Garnier) dyn_spg_ctl: Original code6 !! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code 7 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 8 !!---------------------------------------------------------------------- 9 9 … … 27 27 PRIVATE 28 28 29 PUBLIC dyn_spg! routine called by step module29 PUBLIC dyn_spg ! routine called by step module 30 30 31 !! * module variables32 31 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from lk_dynspg_... 33 32 … … 36 35 # include "vectopt_loop_substitute.h90" 37 36 !!---------------------------------------------------------------------- 38 !! OPA 9.0 , LOCEAN-IPSL (2005)37 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 39 38 !! $Id$ 40 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 39 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 41 40 !!---------------------------------------------------------------------- 42 41 … … 47 46 !! *** ROUTINE dyn_spg *** 48 47 !! 49 !! ** Purpose : compute the lateral ocean dynamics physics. 48 !! ** Purpose : achieve the momentum time stepping by computing the 49 !! last trend, the surface pressure gradient, and performing 50 !! the Leap-Frog integration. 51 !!gm In the current version only the filtered solution provide 52 !!gm the after velocity, in the 2 other (ua,va) are still the trends 53 !! 54 !! ** Method : Three schemes: 55 !! - explicit computation : the spg is evaluated at now 56 !! - filtered computation : the Roulet & madec (2000) technique is used 57 !! - split-explicit computation: a time splitting technique is used 58 !! 59 !! N.B. : When key_esopa is used all the scheme are tested, regardless 60 !! of the physical meaning of the results. 50 61 !!---------------------------------------------------------------------- 51 INTEGER, INTENT( in ) :: kt! ocean time-step index52 INTEGER, INTENT( out ) :: kindic! solver flag62 INTEGER, INTENT(in ) :: kt ! ocean time-step index 63 INTEGER, INTENT( out) :: kindic ! solver flag 53 64 !! 54 REAL(wp) :: z2dt 65 REAL(wp) :: z2dt ! temporary scalar 55 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrdu, ztrdv ! 3D workspace 56 67 !!---------------------------------------------------------------------- 68 69 70 !!gm NOTA BENE : the dynspg_exp and dynspg_ts should be modified so that 71 !!gm they return the after velocity, not the trends (as in trazdf_imp...) 72 !!gm In this case, change/simplify dynnxt 73 74 57 75 58 76 IF( kt == nit000 ) CALL dyn_spg_ctl ! initialisation & control of options … … 65 83 SELECT CASE ( nspg ) ! compute surf. pressure gradient trend and add it to the general trend 66 84 ! 67 CASE ( 0 ) ; CALL dyn_spg_exp 68 CASE ( 1 ) ; CALL dyn_spg_ts 69 CASE ( 2 ) ; CALL dyn_spg_flt 85 CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit 86 CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 87 CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered 70 88 ! 71 CASE ( -1 ) 72 CALL dyn_spg_exp( kt )73 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &74 &tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )75 CALL dyn_spg_ts( kt )76 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &77 &tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )78 CALL dyn_spg_flt( kt, kindic )79 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &80 &tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )89 CASE ( -1 ) ! esopa: test all possibility with control print 90 CALL dyn_spg_exp( kt ) 91 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, & 92 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 93 CALL dyn_spg_ts ( kt ) 94 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, & 95 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 96 CALL dyn_spg_flt( kt, kindic ) 97 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, & 98 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 81 99 END SELECT 82 100 ! 83 IF( l_trddyn ) THEN ! save the horizontal diffusivetrends for further diagnostics101 IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics 84 102 SELECT CASE ( nspg ) 85 103 CASE ( 0, 1 ) … … 106 124 !! 107 125 !! ** Purpose : Control the consistency between cpp options for 108 !! surface pressure gradient schemes126 !! surface pressure gradient schemes 109 127 !!---------------------------------------------------------------------- 110 !! * Local declarations111 128 INTEGER :: ioptio 112 129 !!---------------------------------------------------------------------- 113 130 114 ! Parameter control and print 115 ! --------------------------- 116 ! Control print 117 IF(lwp) THEN 131 IF(lwp) THEN ! Control print 118 132 WRITE(numout,*) 119 133 WRITE(numout,*) 'dyn_spg_ctl : choice of the surface pressure gradient scheme' … … 124 138 ENDIF 125 139 126 ! Control of surface pressure gradient scheme options 127 ! --------------------------------------------------- 140 ! ! Control of surface pressure gradient scheme options 128 141 ioptio = 0 129 142 IF(lk_dynspg_exp) ioptio = ioptio + 1 130 143 IF(lk_dynspg_ts ) ioptio = ioptio + 1 131 144 IF(lk_dynspg_flt) ioptio = ioptio + 1 132 145 ! 133 146 IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ioptio == 0 ) & 134 147 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 135 148 ! 136 149 IF( lk_esopa ) nspg = -1 137 150 IF( lk_dynspg_exp) nspg = 0 138 151 IF( lk_dynspg_ts ) nspg = 1 139 152 IF( lk_dynspg_flt) nspg = 2 140 153 ! 141 154 IF( lk_esopa ) nspg = -1 142 143 IF(lwp) THEN155 ! 156 IF(lwp) THEN 144 157 WRITE(numout,*) 145 158 IF( nspg == -1 ) WRITE(numout,*) ' ESOPA test All scheme used' … … 149 162 ENDIF 150 163 151 ! Control of timestep choice 152 ! -------------------------- 164 ! ! Control of timestep choice 153 165 IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN 154 166 IF( n_cla == 1 ) & … … 157 169 158 170 #if defined key_obc 159 ! Conservation of ocean volume (key_dynspg_flt) 160 ! --------------------------------------------- 161 IF( lk_dynspg_flt ) ln_vol_cst = .true. 171 ! ! Conservation of ocean volume (key_dynspg_flt) 172 IF( lk_dynspg_flt ) ln_vol_cst = .true. 162 173 163 ! Application of Flather's algorithm at open boundaries 164 ! ----------------------------------------------------- 165 IF( lk_dynspg_flt ) ln_obc_fla = .false. 166 IF( lk_dynspg_exp ) ln_obc_fla = .true. 167 IF( lk_dynspg_ts ) ln_obc_fla = .true. 174 ! ! Application of Flather's algorithm at open boundaries 175 IF( lk_dynspg_flt ) ln_obc_fla = .false. 176 IF( lk_dynspg_exp ) ln_obc_fla = .true. 177 IF( lk_dynspg_ts ) ln_obc_fla = .true. 168 178 #endif 169 179 ! 170 180 END SUBROUTINE dyn_spg_ctl 171 181 -
trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
r1528 r1566 1 1 MODULE dynspg_oce 2 !! ----------------------------------------------------------------------2 !!====================================================================== 3 3 !! *** MODULE dynspg_oce *** 4 4 !! 5 !! ** Purpose : Define in memory all the ocean space domain variables 5 !! Ocean dynamics: Define in memory surface pressure gradient variables 6 !!====================================================================== 7 !! History : 1.0 ! 05-12 (C. Talandier, G. Madec) Original code 8 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 6 9 !!---------------------------------------------------------------------- 7 !! Modules used 8 USE par_oce ! ocean parameters 10 USE par_oce ! ocean parameters 9 11 10 12 IMPLICIT NONE 11 13 PUBLIC 12 !!----------------------------------------------------------------------13 !! OPA 9.0 , LOCEAN-IPSL (2005)14 !! $Id$15 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt16 !!----------------------------------------------------------------------17 14 18 !! Surface pressure gradient logicals 19 !! ---------------------------------- 20 #if defined key_dynspg_exp || defined key_esopa 21 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_exp = .TRUE. !: Explicit free surface flag 15 ! !!! Surface pressure gradient logicals 16 #if defined key_dynspg_exp || defined key_esopa 17 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_exp = .TRUE. !: Explicit free surface flag 22 18 #else 23 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_exp = .FALSE. !: Explicit free surface flag19 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_exp = .FALSE. !: Explicit free surface flag 24 20 #endif 25 21 #if defined key_dynspg_ts || defined key_esopa 26 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_ts = .TRUE. !: Free surface with time splitting flag22 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_ts = .TRUE. !: Free surface with time splitting flag 27 23 #else 28 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_ts = .FALSE. !: Free surface with time splitting flag24 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_ts = .FALSE. !: Free surface with time splitting flag 29 25 #endif 30 26 #if defined key_dynspg_flt || defined key_esopa 31 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_flt = .TRUE. !: Filtered free surface cst volume flag27 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_flt = .TRUE. !: Filtered free surface cst volume flag 32 28 #else 33 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_flt = .FALSE. !: Filtered free surface cst volume flag29 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_flt = .FALSE. !: Filtered free surface cst volume flag 34 30 #endif 35 31 36 #if defined key_dynspg_ts || defined key_vvl || defined key_esopa 37 !! Time splitting variables ! variables of the explicit barotropic loop 38 !! ------------------------ 39 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ua_e , va_e ! barotropic velocities (after) 40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshn_e, ssha_e, sshn_b ! sea surface heigth (now, after, average) 41 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hu_e , hv_e ! depth arrays for the barotropic solution 42 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hur_e , hvr_e ! inverse of depth arrays 32 !!gm BUG : always required in _ts, only some of them in vvl 33 ! #if defined key_dynspg_ts || defined key_esopa 34 !!gm end 35 #if defined key_dynspg_ts || defined key_vvl || defined key_esopa 36 ! !!! Time splitting scheme (sub-time step variables) 37 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ua_e , va_e ! barotropic velocities (after) 38 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshn_e, ssha_e, sshn_b ! sea surface heigth (now, after, average) 39 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hu_e , hv_e ! now ocean depth ( = Ho+sshn_e ) 40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hur_e , hvr_e ! inverse of the now depth ( = 1/(Ho+sshn_e) ) 43 41 #endif 44 42 45 43 !!---------------------------------------------------------------------- 46 44 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 45 !! $Id$ 46 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 47 !!====================================================================== 47 48 END MODULE dynspg_oce -
trunk/NEMO/OPA_SRC/LDF/ldftra_c3d.h90
r1152 r1566 37 37 LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout 38 38 39 !! * Local declarations40 INTEGER :: ji, jj, jk ! dummy loop indices41 39 !!---------------------------------------------------------------------- 42 40 -
trunk/NEMO/OPA_SRC/SOL/solmat.F90
r1556 r1566 66 66 !! * Local declarations 67 67 INTEGER :: ji, jj ! dummy loop indices 68 INTEGER :: ii, ij, iiend, ijend ! temporary integers69 68 REAL(wp) :: zcoefs, zcoefw, zcoefe, zcoefn ! temporary scalars 70 69 REAL(wp) :: z2dt, zcoef -
trunk/NEMO/OPA_SRC/istate.F90
r1528 r1566 29 29 USE zdf_oce ! ocean vertical physics 30 30 USE phycst ! physical constants 31 USE wzvmod ! verctical velocity (wzv routine)32 31 USE dtatem ! temperature data (dta_tem routine) 33 32 USE dtasal ! salinity data (dta_sal routine) -
trunk/NEMO/OPA_SRC/stpctl.F90
r1561 r1566 145 145 IF( lk_dynspg_flt ) THEN ! elliptic solver statistics (if required) 146 146 ! 147 IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps ! Solver147 IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps ! Solver 148 148 ! 149 IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN 149 IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN ! create a abort file if problem found 150 150 IF(lwp) THEN 151 151 WRITE(numout,*) ' stpctl: the elliptic solver DO not converge or explode'
Note: See TracChangeset
for help on using the changeset viewer.
