- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5332 r6808 29 29 USE daymod ! calendar 30 30 USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine) 31 USE ldftra _oce ! ocean active tracers: lateral physics31 USE ldftra ! lateral physics: ocean active tracers 32 32 USE zdf_oce ! ocean vertical physics 33 33 USE phycst ! physical constants 34 34 USE dtatsd ! data temperature and salinity (dta_tsd routine) 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 USE zpshde ! partial step: hor. derivative (zps_hde routine)37 USE eosbn2 ! equation of state (eos bn2 routine)38 36 USE domvvl ! varying vertical mesh 39 USE dynspg_oce ! pressure gradient schemes 40 USE dynspg_flt ! filtered free surface 41 USE sol_oce ! ocean solver variables 37 USE iscplrst ! ice sheet coupling 42 38 ! 43 39 USE in_out_manager ! I/O manager … … 54 50 55 51 !! * Substitutions 56 # include "domzgr_substitute.h90"57 52 # include "vectopt_loop_substitute.h90" 58 53 !!---------------------------------------------------------------------- … … 76 71 ! 77 72 78 IF(lwp) WRITE(numout,*) ' '73 IF(lwp) WRITE(numout,*) 79 74 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 80 75 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 81 76 82 CALL dta_tsd_init! Initialisation of T & S input data83 IF( lk_c1d ) CALL dta_uvd_init! Initialization of U & V input data77 CALL dta_tsd_init ! Initialisation of T & S input data 78 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 84 79 85 80 rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk … … 90 85 IF( ln_rstart ) THEN ! Restart from a file 91 86 ! ! ------------------- 92 CALL rst_read ! Read the restart file 93 CALL day_init ! model calendar (using both namelist and restart infos) 87 CALL rst_read ! Read the restart file 88 IF (ln_iscpl) CALL iscpl_stp ! extraloate restart to wet and dry 89 CALL day_init ! model calendar (using both namelist and restart infos) 94 90 ELSE 95 91 ! ! Start from rest … … 103 99 ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp 104 100 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 105 rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp 106 hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp 101 hdivn(:,:,:) = 0._wp 107 102 ! 108 103 IF( cp_cfg == 'eel' ) THEN … … 119 114 ENDIF 120 115 IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 121 CALL wrk_alloc( jpi, jpj, jpk, 2,zuvd )116 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 122 117 CALL dta_uvd( nit000, zuvd ) 123 118 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 124 119 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 125 CALL wrk_dealloc( jpi, jpj, jpk, 2,zuvd )120 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 126 121 ENDIF 127 122 ENDIF 128 123 ! 129 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities 130 #if ! defined key_c1d 131 IF( ln_zps .AND. .NOT. ln_isfcav) & 132 & CALL zps_hde ( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 133 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 134 IF( ln_zps .AND. ln_isfcav) & 135 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 136 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 137 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 138 #endif 139 ! 140 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 141 IF( lk_vvl ) THEN 124 !!gm This is to be changed !!!! 125 ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 126 IF( .NOT.ln_linssh ) THEN 142 127 DO jk = 1, jpk 143 fse3t_b(:,:,jk) = fse3t_n(:,:,jk)144 END DO128 e3t_b(:,:,jk) = e3t_n(:,:,jk) 129 END DO 145 130 ENDIF 131 !!gm 146 132 ! 147 133 ENDIF 148 !149 IF( lk_agrif ) THEN ! read free surface arrays in restart file150 IF( ln_rstart ) THEN151 IF( lk_dynspg_flt ) THEN ! read or initialize the following fields152 ! ! gcx, gcxb for agrif_opa_init153 IF( sol_oce_alloc() > 0 ) CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')154 CALL flt_rst( nit000, 'READ' )155 ENDIF156 ENDIF ! explicit case not coded yet with AGRIF157 ENDIF158 !159 134 ! 160 135 ! Initialize "now" and "before" barotropic velocities: 161 ! Do it whatever the free surface method, these arrays 162 ! being eventually used 163 ! 164 ! 165 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 166 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 167 ! 136 ! Do it whatever the free surface method, these arrays being eventually used 137 ! 138 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 139 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 140 ! 141 !!gm the use of umsak & vmask is not necessary belox as un, vn, ub, vb are always masked 168 142 DO jk = 1, jpkm1 169 143 DO jj = 1, jpj 170 144 DO ji = 1, jpi 171 un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)172 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)145 un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 146 vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 173 147 ! 174 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)175 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)148 ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 149 vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 176 150 END DO 177 151 END DO 178 152 END DO 179 153 ! 180 un_b(:,:) = un_b(:,:) * hur (:,:) 181 vn_b(:,:) = vn_b(:,:) * hvr (:,:) 182 ! 183 ub_b(:,:) = ub_b(:,:) * hur_b(:,:) 184 vb_b(:,:) = vb_b(:,:) * hvr_b(:,:) 185 ! 154 un_b(:,:) = un_b(:,:) * r1_hu_n(:,:) 155 vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:) 156 ! 157 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 158 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 186 159 ! 187 160 IF( nn_timing == 1 ) CALL timing_stop('istate_init') … … 202 175 !! References : Philander ??? 203 176 !!---------------------------------------------------------------------- 204 INTEGER :: ji, jj, jk205 REAL(wp) :: zsal = 35.50 177 INTEGER :: ji, jj, jk 178 REAL(wp) :: zsal = 35.50_wp 206 179 !!---------------------------------------------------------------------- 207 180 ! … … 211 184 ! 212 185 DO jk = 1, jpk 213 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH(( fsdept(:,:,jk)-80.)/30.) ) &214 & + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk)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) 215 188 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 216 189 END DO … … 233 206 !! and relative vorticity fields 234 207 !!---------------------------------------------------------------------- 235 USE div cur ! hor. divergence & rel. vorticity (div_cur routine)208 USE divhor ! hor. divergence (div_hor routine) 236 209 USE iom 237 210 ! … … 265 238 ! 266 239 DO jk = 1, jpk 267 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)240 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk) 268 241 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 269 242 END DO 270 !271 IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &272 & 1 , jpi , 5 , 1 , jpk , &273 & 1 , 1. , numout )274 243 ! 275 244 ! set salinity field to a constant value … … 282 251 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 283 252 ! 284 ! set the dynamics: U,V, hdiv , rot(and ssh if necessary)253 ! set the dynamics: U,V, hdiv (and ssh if necessary) 285 254 ! ---------------- 286 255 ! Start EEL5 configuration with barotropic geostrophic velocities … … 318 287 ENDIF 319 288 ! 320 CALL div_cur( nit000 ) ! horizontal divergence and relative vorticity (curl) 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) 321 292 ! N.B. the vertical velocity will be computed from the horizontal divergence field 322 293 ! in istate by a call to wzv routine … … 338 309 ! 339 310 tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb 340 !341 IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &342 & 1 , jpi , 5 , 1 , jpk , &343 & 1 , 1. , numout )344 311 ! 345 312 ! set salinity field to a constant value … … 371 338 !! 372 339 !! ** Method : - set temprature field 373 !! - set salinity field340 !! - set salinity field 374 341 !!---------------------------------------------------------------------- 375 342 INTEGER :: ji, jj, jk ! dummy loop indices … … 388 355 DO jj = 1, jpj 389 356 DO ji = 1, jpi 390 tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( ( fsdept(ji,jj,jk) - 400) / 700 ) ) &391 & * (-TANH( (500- fsdept(ji,jj,jk)) / 150 ) + 1) / 2 &392 & + ( 15. * ( 1. - TANH( ( fsdept(ji,jj,jk)-50.) / 1500.) ) &393 & - 1.4 * TANH(( fsdept(ji,jj,jk)-100.) / 100.) &394 & + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) &395 & * (-TANH( ( fsdept(ji,jj,jk) - 500) / 150) + 1) / 2357 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 396 363 tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 397 364 tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) 398 365 399 tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( ( fsdept(ji,jj,jk) - 305) / 460 ) ) &400 & * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2 &401 & + ( 35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000. &402 & - 1.62 * TANH( ( fsdept(ji,jj,jk) - 60. ) / 650. ) &403 & + 0.2 * TANH( ( fsdept(ji,jj,jk) - 35. ) / 100. ) &404 & + 0.2 * TANH( ( fsdept(ji,jj,jk) - 1000.) / 5000.) ) &405 & * (-TANH(( fsdept(ji,jj,jk) - 500) / 150) + 1) / 2366 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 406 373 tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 407 374 tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) … … 457 424 !! p=integral [ rau*g dz ] 458 425 !!---------------------------------------------------------------------- 459 USE dynspg ! surface pressure gradient (dyn_spg routine) 460 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 426 USE divhor ! hor. divergence (div_hor routine) 461 427 USE lbclnk ! ocean lateral boundary condition (or mpp link) 462 428 ! 463 429 INTEGER :: ji, jj, jk ! dummy loop indices 464 INTEGER :: indic ! ???465 430 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars 466 431 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn 467 432 !!---------------------------------------------------------------------- 468 433 ! 469 CALL wrk_alloc( jpi, jpj, jpk,zprn)434 CALL wrk_alloc( jpi,jpj,jpk, zprn) 470 435 ! 471 436 IF(lwp) WRITE(numout,*) … … 478 443 zalfg = 0.5 * grav * rau0 479 444 480 zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value445 zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value 481 446 482 447 DO jk = 2, jpkm1 ! Vertical integration from the surface 483 448 zprn(:,:,jk) = zprn(:,:,jk-1) & 484 & + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )449 & + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) ) 485 450 END DO 486 451 … … 529 494 vb(:,:,:) = vn(:,:,:) 530 495 531 ! WARNING !!!!! 532 ! after initializing u and v, we need to calculate the initial streamfunction bsf. 533 ! Otherwise, only the trend will be computed and the model will blow up (inconsistency). 534 ! to do that, we call dyn_spg with a special trick: 535 ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the 536 ! right value assuming the velocities have been set up in one time step. 537 ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.) 538 ! sets up s false trend to calculate the barotropic streamfunction. 539 540 ua(:,:,:) = ub(:,:,:) / rdt 541 va(:,:,:) = vb(:,:,:) / rdt 542 543 ! calls dyn_spg. we assume euler time step, starting from rest. 544 indic = 0 545 CALL dyn_spg( nit000, indic ) ! surface pressure gradient 546 547 ! the new velocity is ua*rdt 548 549 CALL lbc_lnk( ua, 'U', -1. ) 550 CALL lbc_lnk( va, 'V', -1. ) 551 552 ub(:,:,:) = ua(:,:,:) * rdt 553 vb(:,:,:) = va(:,:,:) * rdt 554 ua(:,:,:) = 0.e0 555 va(:,:,:) = 0.e0 556 un(:,:,:) = ub(:,:,:) 557 vn(:,:,:) = vb(:,:,:) 558 559 ! Compute the divergence and curl 560 561 CALL div_cur( nit000 ) ! now horizontal divergence and curl 562 563 hdivb(:,:,:) = hdivn(:,:,:) ! set the before to the now value 564 rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value 565 ! 566 CALL wrk_dealloc( jpi, jpj, jpk, zprn) 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) 567 502 ! 568 503 END SUBROUTINE istate_uvg
Note: See TracChangeset
for help on using the changeset viewer.
