Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r6140 r7646 14 14 !! 3.3 ! 2010-10 (C. Ethe) merge TRC-TRA 15 15 !! 3.4 ! 2011-04 (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 16 !! 3.7 ! 2016-04 (S. Flavoni) introduce user defined initial state 16 17 !!---------------------------------------------------------------------- 17 18 18 19 !!---------------------------------------------------------------------- 19 20 !! istate_init : initial state setting 20 !! istate_tem : analytical profile for initial Temperature21 !! istate_sal : analytical profile for initial Salinity22 !! istate_eel : initial state setting of EEL R5 configuration23 !! istate_gyre : initial state setting of GYRE configuration24 21 !! istate_uvg : initial velocity in geostropic balance 25 22 !!---------------------------------------------------------------------- 26 USE oce ! ocean dynamics and active tracers 27 USE dom_oce ! ocean space and time domain 28 USE c1d ! 1D vertical configuration 29 USE daymod ! calendar 30 USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine) 31 USE ldftra ! lateral physics: ocean active tracers 32 USE zdf_oce ! ocean vertical physics 33 USE phycst ! physical constants 34 USE dtatsd ! data temperature and salinity (dta_tsd routine) 35 USE dtauvd ! data: U & V current (dta_uvd routine) 23 USE oce ! ocean dynamics and active tracers 24 USE dom_oce ! ocean space and time domain 25 USE daymod ! calendar 26 USE divhor ! horizontal divergence (div_hor routine) 27 USE dtatsd ! data temperature and salinity (dta_tsd routine) 28 USE dtauvd ! data: U & V current (dta_uvd routine) 36 29 USE domvvl ! varying vertical mesh 37 30 USE iscplrst ! ice sheet coupling 31 USE wet_dry ! wetting and drying (needed for wad_istate) 32 USE usrdef_istate ! User defined initial state 38 33 ! 39 34 USE in_out_manager ! I/O manager … … 70 65 IF( nn_timing == 1 ) CALL timing_start('istate_init') 71 66 ! 67 IF(lwp) WRITE(numout,*) 68 IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers' 69 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 72 70 73 IF(lwp) WRITE(numout,*) 74 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 75 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 76 71 !!gm Why not include in the first call of dta_tsd ? 72 !!gm probably associated with the use of internal damping... 77 73 CALL dta_tsd_init ! Initialisation of T & S input data 78 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 74 !!gm to be moved in usrdef of C1D case 75 ! IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 76 !!gm 79 77 80 78 rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk … … 86 84 ! ! ------------------- 87 85 CALL rst_read ! Read the restart file 88 IF (ln_iscpl) CALL iscpl_stp ! extra loate restart to wet and dry86 IF (ln_iscpl) CALL iscpl_stp ! extrapolate restart to wet and dry 89 87 CALL day_init ! model calendar (using both namelist and restart infos) 90 ELSE91 !! Start from rest88 ! 89 ELSE ! Start from rest 92 90 ! ! --------------- 93 numror = 0 ! define numror = 0 -> no restart file to read 94 neuler = 0 ! Set time-step indicator at nit000 (euler forward) 95 CALL day_init ! model calendar (using both namelist and restart infos) 96 ! ! Initialization of ocean to zero 97 ! before fields ! now fields 98 sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp 99 ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp 100 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 101 hdivn(:,:,:) = 0._wp 91 numror = 0 ! define numror = 0 -> no restart file to read 92 neuler = 0 ! Set time-step indicator at nit000 (euler forward) 93 CALL day_init ! model calendar (using both namelist and restart infos) 94 ! ! Initialization of ocean to zero 102 95 ! 103 IF( cp_cfg == 'eel' ) THEN 104 CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields 105 ELSEIF( cp_cfg == 'gyre' ) THEN 106 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 107 ELSE ! Initial T-S, U-V fields read in files 108 IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000 109 CALL dta_tsd( nit000, tsb ) 110 tsn(:,:,:,:) = tsb(:,:,:,:) 111 ! 112 ELSE ! Initial T-S fields defined analytically 113 CALL istate_t_s 114 ENDIF 115 IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 116 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 117 CALL dta_uvd( nit000, zuvd ) 118 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 119 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 120 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 121 ENDIF 96 IF( ln_tsd_init ) THEN 97 CALL dta_tsd( nit000, tsb ) ! read 3D T and S data at nit000 98 ! 99 sshb(:,:) = 0._wp ! set the ocean at rest 100 ub (:,:,:) = 0._wp 101 vb (:,:,:) = 0._wp 102 ! 103 ELSE ! user defined initial T and S 104 CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb ) 122 105 ENDIF 106 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones 107 sshn (:,:) = sshb(:,:) 108 un (:,:,:) = ub (:,:,:) 109 vn (:,:,:) = vb (:,:,:) 110 hdivn(:,:,jpk) = 0._wp ! bottom divergence set one for 0 to zero at jpk level 111 CALL div_hor( 0 ) ! compute interior hdivn value 112 !!gm hdivn(:,:,:) = 0._wp 113 114 !!gm POTENTIAL BUG : 115 !!gm ISSUE : if sshb /= 0 then, in non linear free surface, the e3._n, e3._b should be recomputed 116 !! as well as gdept and gdepw.... !!!!! 117 !! ===>>>> probably a call to domvvl initialisation here.... 118 119 120 ! 121 !!gm to be moved in usrdef of C1D case 122 ! IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 123 ! CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 124 ! CALL dta_uvd( nit000, zuvd ) 125 ! ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 126 ! vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 127 ! CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 128 ! ENDIF 123 129 ! 124 130 !!gm This is to be changed !!!! 125 ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here126 IF( .NOT.ln_linssh ) THEN127 DO jk = 1, jpk128 e3t_b(:,:,jk) = e3t_n(:,:,jk)129 END DO130 ENDIF131 ! ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 132 ! IF( .NOT.ln_linssh ) THEN 133 ! DO jk = 1, jpk 134 ! e3t_b(:,:,jk) = e3t_n(:,:,jk) 135 ! END DO 136 ! ENDIF 131 137 !!gm 132 138 ! 133 ENDIF 139 ENDIF 134 140 ! 135 141 ! Initialize "now" and "before" barotropic velocities: … … 139 145 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 140 146 ! 141 !!gm the use of umsak & vmask is not necessary belo xas un, vn, ub, vb are always masked147 !!gm the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked 142 148 DO jk = 1, jpkm1 143 149 DO jj = 1, jpj … … 162 168 END SUBROUTINE istate_init 163 169 164 165 SUBROUTINE istate_t_s 166 !!--------------------------------------------------------------------- 167 !! *** ROUTINE istate_t_s *** 168 !! 169 !! ** Purpose : Intialization of the temperature field with an 170 !! analytical profile or a file (i.e. in EEL configuration) 171 !! 172 !! ** Method : - temperature: use Philander analytic profile 173 !! - salinity : use to a constant value 35.5 174 !! 175 !! References : Philander ??? 176 !!---------------------------------------------------------------------- 177 INTEGER :: ji, jj, jk 178 REAL(wp) :: zsal = 35.50_wp 179 !!---------------------------------------------------------------------- 180 ! 181 IF(lwp) WRITE(numout,*) 182 IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile' 183 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)' 184 ! 185 DO jk = 1, jpk 186 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) ) & 187 & + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.) ) * tmask(:,:,jk) 188 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 189 END DO 190 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 191 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 192 ! 193 END SUBROUTINE istate_t_s 194 195 196 SUBROUTINE istate_eel 197 !!---------------------------------------------------------------------- 198 !! *** ROUTINE istate_eel *** 199 !! 200 !! ** Purpose : Initialization of the dynamics and tracers for EEL R5 201 !! configuration (channel with or without a topographic bump) 202 !! 203 !! ** Method : - set temprature field 204 !! - set salinity field 205 !! - set velocity field including horizontal divergence 206 !! and relative vorticity fields 207 !!---------------------------------------------------------------------- 208 USE divhor ! hor. divergence (div_hor routine) 209 USE iom 210 ! 211 INTEGER :: inum ! temporary logical unit 212 INTEGER :: ji, jj, jk ! dummy loop indices 213 INTEGER :: ijloc 214 REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars 215 REAL(wp) :: zt1 = 15._wp ! surface temperature value (EEL R5) 216 REAL(wp) :: zt2 = 5._wp ! bottom temperature value (EEL R5) 217 REAL(wp) :: zsal = 35.0_wp ! constant salinity (EEL R2, R5 and R6) 218 REAL(wp) :: zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5) 219 REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain 220 !!---------------------------------------------------------------------- 221 ! 222 SELECT CASE ( jp_cfg ) 223 ! ! ==================== 224 CASE ( 5 ) ! EEL R5 configuration 225 ! ! ==================== 226 ! 227 ! set temperature field with a linear profile 228 ! ------------------------------------------- 229 IF(lwp) WRITE(numout,*) 230 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile' 231 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 232 ! 233 zh1 = gdept_1d( 1 ) 234 zh2 = gdept_1d(jpkm1) 235 ! 236 zslope = ( zt1 - zt2 ) / ( zh1 - zh2 ) 237 zcst = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 ) 238 ! 239 DO jk = 1, jpk 240 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk) 241 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 242 END DO 243 ! 244 ! set salinity field to a constant value 245 ! -------------------------------------- 246 IF(lwp) WRITE(numout,*) 247 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal 248 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 249 ! 250 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 251 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 252 ! 253 ! set the dynamics: U,V, hdiv (and ssh if necessary) 254 ! ---------------- 255 ! Start EEL5 configuration with barotropic geostrophic velocities 256 ! according the sshb and sshn SSH imposed. 257 ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y) 258 ! we use the Coriolis frequency at mid-channel. 259 ub(:,:,:) = zueel * umask(:,:,:) 260 un(:,:,:) = ub(:,:,:) 261 ijloc = mj0(INT(jpjglo-1)/2) 262 zfcor = ff(1,ijloc) 263 ! 264 DO jj = 1, jpjglo 265 zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 266 END DO 267 ! 268 IF(lwp) THEN 269 WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel 270 WRITE(numout,*) ' Geostrophic SSH profile as a function of y:' 271 WRITE(numout,'(12(1x,f6.2))') zssh(1,:) 272 ENDIF 273 ! 274 DO jj = 1, nlcj 275 DO ji = 1, nlci 276 sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1) 277 END DO 278 END DO 279 sshb(nlci+1:jpi, : ) = 0.e0 ! set to zero extra mpp columns 280 sshb( : ,nlcj+1:jpj) = 0.e0 ! set to zero extra mpp rows 281 ! 282 sshn(:,:) = sshb(:,:) ! set now ssh to the before value 283 ! 284 IF( nn_rstssh /= 0 ) THEN 285 nn_rstssh = 0 ! hand-made initilization of ssh 286 CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' ) 287 ENDIF 288 ! 289 !!gm Check here call to div_hor should not be necessary 290 !!gm div_hor call runoffs not sure they are defined at that level 291 CALL div_hor( nit000 ) ! horizontal divergence and relative vorticity (curl) 292 ! N.B. the vertical velocity will be computed from the horizontal divergence field 293 ! in istate by a call to wzv routine 294 295 296 ! ! ========================== 297 CASE ( 2 , 6 ) ! EEL R2 or R6 configuration 298 ! ! ========================== 299 ! 300 ! set temperature field with a NetCDF file 301 ! ---------------------------------------- 302 IF(lwp) WRITE(numout,*) 303 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file' 304 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 305 ! 306 CALL iom_open ( 'eel.initemp', inum ) 307 CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb) 308 CALL iom_close( inum ) 309 ! 310 tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb 311 ! 312 ! set salinity field to a constant value 313 ! -------------------------------------- 314 IF(lwp) WRITE(numout,*) 315 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal 316 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 317 ! 318 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 319 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 320 ! 321 ! ! =========================== 322 CASE DEFAULT ! NONE existing configuration 323 ! ! =========================== 324 WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded' 325 CALL ctl_stop( ctmp1 ) 326 ! 327 END SELECT 328 ! 329 END SUBROUTINE istate_eel 330 331 332 SUBROUTINE istate_gyre 333 !!---------------------------------------------------------------------- 334 !! *** ROUTINE istate_gyre *** 335 !! 336 !! ** Purpose : Initialization of the dynamics and tracers for GYRE 337 !! configuration (double gyre with rotated domain) 338 !! 339 !! ** Method : - set temprature field 340 !! - set salinity field 341 !!---------------------------------------------------------------------- 342 INTEGER :: ji, jj, jk ! dummy loop indices 343 INTEGER :: inum ! temporary logical unit 344 INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization 345 !!---------------------------------------------------------------------- 346 ! 347 SELECT CASE ( ntsinit) 348 ! 349 CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS 350 IF(lwp) WRITE(numout,*) 351 IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS ' 352 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 353 ! 354 DO jk = 1, jpk 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) ) & 358 & * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2 & 359 & + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) ) & 360 & - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.) & 361 & + 7. * (1500. - gdept_n(ji,jj,jk)) / 1500. ) & 362 & * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 363 tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 364 tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) 365 366 tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) ) & 367 & * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2 & 368 & + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000. & 369 & - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60. ) / 650. ) & 370 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 35. ) / 100. ) & 371 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) ) & 372 & * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 373 tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 374 tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) 375 END DO 376 END DO 377 END DO 378 ! 379 CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files 380 IF(lwp) WRITE(numout,*) 381 IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files' 382 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 383 IF(lwp) WRITE(numout,*) ' NetCDF FORMAT' 384 385 ! Read temperature field 386 ! ---------------------- 387 CALL iom_open ( 'data_tem', inum ) 388 CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 389 CALL iom_close( inum ) 390 391 tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 392 tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 393 394 ! Read salinity field 395 ! ------------------- 396 CALL iom_open ( 'data_sal', inum ) 397 CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 398 CALL iom_close( inum ) 399 400 tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 401 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 402 ! 403 END SELECT 404 ! 405 IF(lwp) THEN 406 WRITE(numout,*) 407 WRITE(numout,*) ' Initial temperature and salinity profiles:' 408 WRITE(numout, "(9x,' level gdept_1d temperature salinity ')" ) 409 WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk ) 410 ENDIF 411 ! 412 END SUBROUTINE istate_gyre 413 414 415 SUBROUTINE istate_uvg 416 !!---------------------------------------------------------------------- 417 !! *** ROUTINE istate_uvg *** 418 !! 419 !! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields 420 !! 421 !! ** Method : Using the hydrostatic hypothesis the now hydrostatic 422 !! pressure is computed by integrating the in-situ density from the 423 !! surface to the bottom. 424 !! p=integral [ rau*g dz ] 425 !!---------------------------------------------------------------------- 426 USE divhor ! hor. divergence (div_hor routine) 427 USE lbclnk ! ocean lateral boundary condition (or mpp link) 428 ! 429 INTEGER :: ji, jj, jk ! dummy loop indices 430 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars 431 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn 432 !!---------------------------------------------------------------------- 433 ! 434 CALL wrk_alloc( jpi,jpj,jpk, zprn) 435 ! 436 IF(lwp) WRITE(numout,*) 437 IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy' 438 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 439 440 ! Compute the now hydrostatic pressure 441 ! ------------------------------------ 442 443 zalfg = 0.5 * grav * rau0 444 445 zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value 446 447 DO jk = 2, jpkm1 ! Vertical integration from the surface 448 zprn(:,:,jk) = zprn(:,:,jk-1) & 449 & + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) ) 450 END DO 451 452 ! Compute geostrophic balance 453 ! --------------------------- 454 DO jk = 1, jpkm1 455 DO jj = 2, jpjm1 456 DO ji = fs_2, fs_jpim1 ! vertor opt. 457 zmsv = 1. / MAX( umask(ji-1,jj+1,jk) + umask(ji ,jj+1,jk) & 458 + umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) , 1. ) 459 zphv = ( zprn(ji ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1) & 460 + ( zprn(ji+1,jj+1,jk) - zprn(ji ,jj+1,jk) ) * umask(ji ,jj+1,jk) / e1u(ji ,jj+1) & 461 + ( zprn(ji ,jj ,jk) - zprn(ji-1,jj ,jk) ) * umask(ji-1,jj ,jk) / e1u(ji-1,jj ) & 462 + ( zprn(ji+1,jj ,jk) - zprn(ji ,jj ,jk) ) * umask(ji ,jj ,jk) / e1u(ji ,jj ) 463 zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk) 464 465 zmsu = 1. / MAX( vmask(ji+1,jj ,jk) + vmask(ji ,jj ,jk) & 466 + vmask(ji+1,jj-1,jk) + vmask(ji ,jj-1,jk) , 1. ) 467 zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj ,jk) ) * vmask(ji+1,jj ,jk) / e2v(ji+1,jj ) & 468 + ( zprn(ji ,jj+1,jk) - zprn(ji ,jj ,jk) ) * vmask(ji ,jj ,jk) / e2v(ji ,jj ) & 469 + ( zprn(ji+1,jj ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1) & 470 + ( zprn(ji ,jj ,jk) - zprn(ji ,jj-1,jk) ) * vmask(ji ,jj-1,jk) / e2v(ji ,jj-1) 471 zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk) 472 473 ! Compute the geostrophic velocities 474 un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji ,jj-1) ) 475 vn(ji,jj,jk) = 2. * zphv / ( ff(ji,jj) + ff(ji-1,jj ) ) 476 END DO 477 END DO 478 END DO 479 480 IF(lwp) WRITE(numout,*) ' we force to zero bottom velocity' 481 482 ! Susbtract the bottom velocity (level jpk-1 for flat bottom case) 483 ! to have a zero bottom velocity 484 485 DO jk = 1, jpkm1 486 un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk) 487 vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk) 488 END DO 489 490 CALL lbc_lnk( un, 'U', -1. ) 491 CALL lbc_lnk( vn, 'V', -1. ) 492 493 ub(:,:,:) = un(:,:,:) 494 vb(:,:,:) = vn(:,:,:) 495 496 ! 497 !!gm Check here call to div_hor should not be necessary 498 !!gm div_hor call runoffs not sure they are defined at that level 499 CALL div_hor( nit000 ) ! now horizontal divergence 500 ! 501 CALL wrk_dealloc( jpi,jpj,jpk, zprn) 502 ! 503 END SUBROUTINE istate_uvg 504 505 !!===================================================================== 170 !!====================================================================== 506 171 END MODULE istate
Note: See TracChangeset
for help on using the changeset viewer.
