- Timestamp:
- 2016-08-01T15:37:15+02:00 (8 years ago)
- Location:
- branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 21 deleted
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r5836 r6827 21 21 USE phycst ! physical constants 22 22 USE in_out_manager ! I/O manager 23 USE sbc_oce ! ocean surface boundary conditions24 23 USE lib_fortran, ONLY: glob_sum, DDPDD 25 24 USE lbclnk ! lateral boundary condition - MPP exchanges … … 31 30 32 31 PUBLIC dom_clo ! routine called by domain module 33 PUBLIC sbc_clo ! routine called by step module34 PUBLIC clo_rnf ! routine called by sbcrnf module35 PUBLIC clo_ups ! routine called in traadv_cen2(_jki) module36 32 PUBLIC clo_bat ! routine called in domzgr module 37 33 … … 185 181 186 182 187 SUBROUTINE sbc_clo( kt )188 !!---------------------------------------------------------------------189 !! *** ROUTINE sbc_clo ***190 !!191 !! ** Purpose : Special handling of closed seas192 !!193 !! ** Method : Water flux is forced to zero over closed sea194 !! Excess is shared between remaining ocean, or195 !! put as run-off in open ocean.196 !!197 !! ** Action : emp updated surface freshwater fluxes and associated heat content at kt198 !!----------------------------------------------------------------------199 INTEGER, INTENT(in) :: kt ! ocean model time step200 !201 INTEGER :: ji, jj, jc, jn ! dummy loop indices202 REAL(wp), PARAMETER :: rsmall = 1.e-20_wp ! Closed sea correction epsilon203 REAL(wp) :: zze2, ztmp, zcorr !204 REAL(wp) :: zcoef, zcoef1 !205 COMPLEX(wp) :: ctmp206 REAL(wp), DIMENSION(jpncs) :: zfwf ! 1D workspace207 !!----------------------------------------------------------------------208 !209 IF( nn_timing == 1 ) CALL timing_start('sbc_clo')210 ! !------------------!211 IF( kt == nit000 ) THEN ! Initialisation !212 ! !------------------!213 IF(lwp) WRITE(numout,*)214 IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '215 IF(lwp) WRITE(numout,*)'~~~~~~~'216 217 surf(:) = 0.e0_wp218 !219 surf(jpncs+1) = glob_sum( e1e2t(:,:) ) ! surface of the global ocean220 !221 ! ! surface of closed seas222 IF( lk_mpp_rep ) THEN ! MPP reproductible calculation223 DO jc = 1, jpncs224 ctmp = CMPLX( 0.e0, 0.e0, wp )225 DO jj = ncsj1(jc), ncsj2(jc)226 DO ji = ncsi1(jc), ncsi2(jc)227 ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)228 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )229 END DO230 END DO231 IF( lk_mpp ) CALL mpp_sum( ctmp )232 surf(jc) = REAL(ctmp,wp)233 END DO234 ELSE ! Standard calculation235 DO jc = 1, jpncs236 DO jj = ncsj1(jc), ncsj2(jc)237 DO ji = ncsi1(jc), ncsi2(jc)238 surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj) ! surface of closed seas239 END DO240 END DO241 END DO242 IF( lk_mpp ) CALL mpp_sum ( surf, jpncs ) ! mpp: sum over all the global domain243 ENDIF244 245 IF(lwp) WRITE(numout,*)' Closed sea surfaces'246 DO jc = 1, jpncs247 IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)248 END DO249 250 ! jpncs+1 : surface of sea, closed seas excluded251 DO jc = 1, jpncs252 surf(jpncs+1) = surf(jpncs+1) - surf(jc)253 END DO254 !255 ENDIF256 ! !--------------------!257 ! ! update emp !258 zfwf = 0.e0_wp !--------------------!259 IF( lk_mpp_rep ) THEN ! MPP reproductible calculation260 DO jc = 1, jpncs261 ctmp = CMPLX( 0.e0, 0.e0, wp )262 DO jj = ncsj1(jc), ncsj2(jc)263 DO ji = ncsi1(jc), ncsi2(jc)264 ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)265 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )266 END DO267 END DO268 IF( lk_mpp ) CALL mpp_sum( ctmp )269 zfwf(jc) = REAL(ctmp,wp)270 END DO271 ELSE ! Standard calculation272 DO jc = 1, jpncs273 DO jj = ncsj1(jc), ncsj2(jc)274 DO ji = ncsi1(jc), ncsi2(jc)275 zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)276 END DO277 END DO278 END DO279 IF( lk_mpp ) CALL mpp_sum ( zfwf(:) , jpncs ) ! mpp: sum over all the global domain280 ENDIF281 282 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! Black Sea case for ORCA_R2 configuration283 zze2 = ( zfwf(3) + zfwf(4) ) * 0.5_wp284 zfwf(3) = zze2285 zfwf(4) = zze2286 ENDIF287 288 zcorr = 0._wp289 290 DO jc = 1, jpncs291 !292 ! The following if avoids the redistribution of the round off293 IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN294 !295 IF( ncstt(jc) == 0 ) THEN ! water/evap excess is shared by all open ocean296 zcoef = zfwf(jc) / surf(jpncs+1)297 zcoef1 = rcp * zcoef298 emp(:,:) = emp(:,:) + zcoef299 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)300 ! accumulate closed seas correction301 zcorr = zcorr + zcoef302 !303 ELSEIF( ncstt(jc) == 1 ) THEN ! Excess water in open sea, at outflow location, excess evap shared304 IF ( zfwf(jc) <= 0.e0_wp ) THEN305 DO jn = 1, ncsnr(jc)306 ji = mi0(ncsir(jc,jn))307 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean308 IF ( ji > 1 .AND. ji < jpi &309 .AND. jj > 1 .AND. jj < jpj ) THEN310 zcoef = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )311 zcoef1 = rcp * zcoef312 emp(ji,jj) = emp(ji,jj) + zcoef313 qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)314 ENDIF315 END DO316 ELSE317 zcoef = zfwf(jc) / surf(jpncs+1)318 zcoef1 = rcp * zcoef319 emp(:,:) = emp(:,:) + zcoef320 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)321 ! accumulate closed seas correction322 zcorr = zcorr + zcoef323 ENDIF324 ELSEIF( ncstt(jc) == 2 ) THEN ! Excess e-p-r (either sign) goes to open ocean, at outflow location325 DO jn = 1, ncsnr(jc)326 ji = mi0(ncsir(jc,jn))327 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean328 IF( ji > 1 .AND. ji < jpi &329 .AND. jj > 1 .AND. jj < jpj ) THEN330 zcoef = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )331 zcoef1 = rcp * zcoef332 emp(ji,jj) = emp(ji,jj) + zcoef333 qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)334 ENDIF335 END DO336 ENDIF337 !338 DO jj = ncsj1(jc), ncsj2(jc)339 DO ji = ncsi1(jc), ncsi2(jc)340 zcoef = zfwf(jc) / surf(jc)341 zcoef1 = rcp * zcoef342 emp(ji,jj) = emp(ji,jj) - zcoef343 qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)344 END DO345 END DO346 !347 END IF348 END DO349 350 IF ( ABS(zcorr) > rsmall ) THEN ! remove the global correction from the closed seas351 DO jc = 1, jpncs ! only if it is large enough352 DO jj = ncsj1(jc), ncsj2(jc)353 DO ji = ncsi1(jc), ncsi2(jc)354 emp(ji,jj) = emp(ji,jj) - zcorr355 qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)356 END DO357 END DO358 END DO359 ENDIF360 !361 emp (:,:) = emp (:,:) * tmask(:,:,1)362 !363 CALL lbc_lnk( emp , 'T', 1._wp )364 !365 IF( nn_timing == 1 ) CALL timing_stop('sbc_clo')366 !367 END SUBROUTINE sbc_clo368 369 370 SUBROUTINE clo_rnf( p_rnfmsk )371 !!---------------------------------------------------------------------372 !! *** ROUTINE sbc_rnf ***373 !!374 !! ** Purpose : allow the treatment of closed sea outflow grid-points375 !! to be the same as river mouth grid-points376 !!377 !! ** Method : set to 1 the runoff mask (mskrnf, see sbcrnf module)378 !! at the closed sea outflow grid-point.379 !!380 !! ** Action : update (p_)mskrnf (set 1 at closed sea outflow)381 !!----------------------------------------------------------------------382 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: p_rnfmsk ! river runoff mask (rnfmsk array)383 !384 INTEGER :: jc, jn, ji, jj ! dummy loop indices385 !!----------------------------------------------------------------------386 !387 DO jc = 1, jpncs388 IF( ncstt(jc) >= 1 ) THEN ! runoff mask set to 1 at closed sea outflows389 DO jn = 1, 4390 DO jj = mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) )391 DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) )392 p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp )393 END DO394 END DO395 END DO396 ENDIF397 END DO398 !399 END SUBROUTINE clo_rnf400 401 402 SUBROUTINE clo_ups( p_upsmsk )403 !!---------------------------------------------------------------------404 !! *** ROUTINE sbc_rnf ***405 !!406 !! ** Purpose : allow the treatment of closed sea outflow grid-points407 !! to be the same as river mouth grid-points408 !!409 !! ** Method : set to 0.5 the upstream mask (upsmsk, see traadv_cen2410 !! module) over the closed seas.411 !!412 !! ** Action : update (p_)upsmsk (set 0.5 over closed seas)413 !!----------------------------------------------------------------------414 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: p_upsmsk ! upstream mask (upsmsk array)415 !416 INTEGER :: jc, ji, jj ! dummy loop indices417 !!----------------------------------------------------------------------418 !419 DO jc = 1, jpncs420 DO jj = ncsj1(jc), ncsj2(jc)421 DO ji = ncsi1(jc), ncsi2(jc)422 p_upsmsk(ji,jj) = 0.5_wp ! mixed upstream/centered scheme over closed seas423 END DO424 END DO425 END DO426 !427 END SUBROUTINE clo_ups428 429 430 183 SUBROUTINE clo_bat( pbat, kbat ) 431 184 !!--------------------------------------------------------------------- -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r6140 r6827 33 33 USE ioipsl , ONLY : ymds2ju ! for calendar 34 34 USE prtctl ! Print control 35 USE trc_oce , ONLY : lk_offline ! offline flag36 35 USE timing ! Timing 37 USE restart ! restart38 36 39 37 IMPLICIT NONE … … 88 86 ndt05 = NINT(0.5 * rdt ) 89 87 90 IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' )91 88 92 89 ! set the calandar from ndastp (read in restart file and namelist) … … 281 278 ENDIF 282 279 283 IF( .NOT. lk_offline ) CALL rst_opn( kt ) ! Open the restart file if needed and control lrst_oce284 IF( lrst_oce ) CALL day_rst( kt, 'WRITE' ) ! write day restart information285 280 ! 286 281 IF( nn_timing == 1 ) CALL timing_stop('day') … … 288 283 END SUBROUTINE day 289 284 290 291 SUBROUTINE day_rst( kt, cdrw )292 !!---------------------------------------------------------------------293 !! *** ROUTINE ts_rst ***294 !!295 !! ** Purpose : Read or write calendar in restart file:296 !!297 !! WRITE(READ) mode:298 !! kt : number of time step since the begining of the experiment at the299 !! end of the current(previous) run300 !! adatrj(0) : number of elapsed days since the begining of the experiment at the301 !! end of the current(previous) run (REAL -> keep fractions of day)302 !! ndastp : date at the end of the current(previous) run (coded as yyyymmdd integer)303 !!304 !! According to namelist parameter nrstdt,305 !! nrstdt = 0 no control on the date (nit000 is arbitrary).306 !! nrstdt = 1 we verify that nit000 is equal to the last307 !! time step of previous run + 1.308 !! In both those options, the exact duration of the experiment309 !! since the beginning (cumulated duration of all previous restart runs)310 !! is not stored in the restart and is assumed to be (nit000-1)*rdt.311 !! This is valid is the time step has remained constant.312 !!313 !! nrstdt = 2 the duration of the experiment in days (adatrj)314 !! has been stored in the restart file.315 !!----------------------------------------------------------------------316 INTEGER , INTENT(in) :: kt ! ocean time-step317 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag318 !319 REAL(wp) :: zkt, zndastp, zdayfrac, ksecs, ktime320 INTEGER :: ihour, iminute321 !!----------------------------------------------------------------------322 323 IF( TRIM(cdrw) == 'READ' ) THEN324 325 IF( iom_varid( numror, 'kt', ldstop = .FALSE. ) > 0 ) THEN326 ! Get Calendar informations327 CALL iom_get( numror, 'kt', zkt ) ! last time-step of previous run328 IF(lwp) THEN329 WRITE(numout,*) ' *** Info read in restart : '330 WRITE(numout,*) ' previous time-step : ', NINT( zkt )331 WRITE(numout,*) ' *** restart option'332 SELECT CASE ( nrstdt )333 CASE ( 0 ) ; WRITE(numout,*) ' nrstdt = 0 : no control of nit000'334 CASE ( 1 ) ; WRITE(numout,*) ' nrstdt = 1 : no control the date at nit000 (use ndate0 read in the namelist)'335 CASE ( 2 ) ; WRITE(numout,*) ' nrstdt = 2 : calendar parameters read in restart'336 END SELECT337 WRITE(numout,*)338 ENDIF339 ! Control of date340 IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 ) &341 & CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart', &342 & ' verify the restart file or rerun with nrstdt = 0 (namelist)' )343 ! define ndastp and adatrj344 IF ( nrstdt == 2 ) THEN345 ! read the parameters corresponding to nit000 - 1 (last time step of previous run)346 CALL iom_get( numror, 'ndastp', zndastp )347 ndastp = NINT( zndastp )348 CALL iom_get( numror, 'adatrj', adatrj )349 CALL iom_get( numror, 'ntime', ktime )350 nn_time0=INT(ktime)351 ! calculate start time in hours and minutes352 zdayfrac=adatrj-INT(adatrj)353 ksecs = NINT(zdayfrac*86400) ! Nearest second to catch rounding errors in adatrj354 ihour = INT(ksecs/3600)355 iminute = ksecs/60-ihour*60356 357 ! Add to nn_time0358 nhour = nn_time0 / 100359 nminute = ( nn_time0 - nhour * 100 )360 nminute=nminute+iminute361 362 IF( nminute >= 60 ) THEN363 nminute=nminute-60364 nhour=nhour+1365 ENDIF366 nhour=nhour+ihour367 IF( nhour >= 24 ) THEN368 nhour=nhour-24369 adatrj=adatrj+1370 ENDIF371 nn_time0 = nhour * 100 + nminute372 adatrj = INT(adatrj) ! adatrj set to integer as nn_time0 updated373 ELSE374 ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)375 ndastp = ndate0 ! ndate0 read in the namelist in dom_nam376 nhour = nn_time0 / 100377 nminute = ( nn_time0 - nhour * 100 )378 IF( nhour*3600+nminute*60-ndt05 .lt. 0 ) ndastp=ndastp-1 ! Start hour is specified in the namelist (default 0)379 adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday380 ! note this is wrong if time step has changed during run381 ENDIF382 ELSE383 ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day)384 ndastp = ndate0 ! ndate0 read in the namelist in dom_nam385 nhour = nn_time0 / 100386 nminute = ( nn_time0 - nhour * 100 )387 IF( nhour*3600+nminute*60-ndt05 .lt. 0 ) ndastp=ndastp-1 ! Start hour is specified in the namelist (default 0)388 adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday389 ENDIF390 IF( ABS(adatrj - REAL(NINT(adatrj),wp)) < 0.1 / rday ) adatrj = REAL(NINT(adatrj),wp) ! avoid truncation error391 !392 IF(lwp) THEN393 WRITE(numout,*) ' *** Info used values : '394 WRITE(numout,*) ' date ndastp : ', ndastp395 WRITE(numout,*) ' number of elapsed days since the begining of run : ', adatrj396 WRITE(numout,*) ' nn_time0 : ',nn_time0397 WRITE(numout,*)398 ENDIF399 !400 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN401 !402 IF( kt == nitrst ) THEN403 IF(lwp) WRITE(numout,*)404 IF(lwp) WRITE(numout,*) 'rst_write : write oce restart file kt =', kt405 IF(lwp) WRITE(numout,*) '~~~~~~~'406 ENDIF407 ! calendar control408 CALL iom_rstput( kt, nitrst, numrow, 'kt' , REAL( kt , wp) ) ! time-step409 CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) ) ! date410 CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj ) ! number of elapsed days since411 ! ! the begining of the run [s]412 CALL iom_rstput( kt, nitrst, numrow, 'ntime' , REAL( nn_time0, wp) ) ! time413 ENDIF414 !415 END SUBROUTINE day_rst416 285 417 286 !!====================================================================== -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r6140 r6827 166 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point 167 167 ! 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_f, ff_t !: coriolis factor [1/s] 169 169 170 170 !!---------------------------------------------------------------------- … … 332 332 & e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj) , & 333 333 & e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj) , & 334 & ff (jpi,jpj), STAT=ierr(3) )334 & ff_f(jpi,jpj) , ff_t(jpi,jpj) , STAT=ierr(3) ) 335 335 ! 336 336 ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) , & -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r6140 r6827 24 24 USE oce ! ocean variables 25 25 USE dom_oce ! domain: ocean 26 USE sbc_oce ! surface boundary condition: ocean27 26 USE phycst ! physical constants 28 27 USE closea ! closed seas … … 33 32 USE domwri ! domain: write the meshmask file 34 33 USE domvvl ! variable volume 35 USE c1d ! 1D vertical configuration36 USE dyncor_c1d ! Coriolis term (c1d case) (cor_c1d routine)37 34 ! 38 35 USE in_out_manager ! I/O manager 36 USE iom ! 39 37 USE wrk_nemo ! Memory Allocation 40 38 USE lib_mpp ! distributed memory computing library … … 137 135 ENDIF 138 136 ! 139 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 140 ! 141 CALL dom_stp ! time step 142 IF( nmsh /= 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file 143 IF( nmsh /= 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file 144 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 137 CALL cfg_write ! create the configuration file 145 138 ! 146 139 IF( nn_timing == 1 ) CALL timing_stop('dom_init') … … 465 458 END SUBROUTINE dom_stiff 466 459 460 461 SUBROUTINE cfg_write 462 !!---------------------------------------------------------------------- 463 !! *** ROUTINE cfg_write *** 464 !! 465 !! ** Purpose : Create the "domain_cfg" file, a NetCDF file which 466 !! contains all the ocean domain informations required to 467 !! define an ocean configuration. 468 !! 469 !! ** Method : Write in a file all the arrays required to set up an 470 !! ocean configuration. 471 !! 472 !! ** output file : domain_cfg.nc : domain size, characteristics, 473 !horizontal mesh, 474 !! Coriolis parameter, depth and vertical 475 !scale factors 476 !!---------------------------------------------------------------------- 477 INTEGER :: ji, jj, jk ! dummy loop indices 478 INTEGER :: izco, izps, isco, icav 479 INTEGER :: inum ! temprary units for 'domain_cfg.nc' file 480 CHARACTER(len=21) :: clnam ! filename (mesh and mask informations) 481 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! workspace 482 !!---------------------------------------------------------------------- 483 ! 484 IF(lwp) WRITE(numout,*) 485 IF(lwp) WRITE(numout,*) 'cfg_write : create the "domain_cfg.nc" file containing all required configuration information' 486 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 487 ! 488 ! ! ============================= ! 489 ! ! create 'domain_cfg.nc' file ! 490 ! ! ============================= ! 491 ! 492 clnam = 'domain_cfg' ! filename (configuration information) 493 CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib ) 494 495 ! !== global domain size ==! 496 ! 497 CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 ) 498 CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 ) 499 CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk , wp), ktype = jp_i4 ) 500 ! 501 ! !== domain characteristics ==! 502 ! 503 ! ! lateral boundary of the global 504 ! domain 505 CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 ) 506 ! 507 ! ! type of vertical coordinate 508 IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF 509 IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF 510 IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF 511 CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 ) 512 CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 ) 513 CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 ) 514 ! 515 ! ! ocean cavities under iceshelves 516 IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF 517 CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 518 ! 519 ! !== horizontal mesh ! 520 ! 521 CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 ) ! latitude 522 CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 ) 523 CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 ) 524 CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 ) 525 ! 526 CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 ) ! longitude 527 CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 ) 528 CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 ) 529 CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 ) 530 ! 531 CALL iom_rstput( 0, 0, inum, 'e1t' , e1t , ktype = jp_r8 ) ! i-scale factors (e1.) 532 CALL iom_rstput( 0, 0, inum, 'e1u' , e1u , ktype = jp_r8 ) 533 CALL iom_rstput( 0, 0, inum, 'e1v' , e1v , ktype = jp_r8 ) 534 CALL iom_rstput( 0, 0, inum, 'e1f' , e1f , ktype = jp_r8 ) 535 ! 536 CALL iom_rstput( 0, 0, inum, 'e2t' , e2t , ktype = jp_r8 ) ! j-scale factors (e2.) 537 CALL iom_rstput( 0, 0, inum, 'e2u' , e2u , ktype = jp_r8 ) 538 CALL iom_rstput( 0, 0, inum, 'e2v' , e2v , ktype = jp_r8 ) 539 CALL iom_rstput( 0, 0, inum, 'e2f' , e2f , ktype = jp_r8 ) 540 ! 541 CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 ) ! coriolis factor 542 CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 ) 543 ! 544 ! !== vertical mesh - 3D mask ==! 545 ! 546 CALL iom_rstput( 0, 0, inum, 'gdept_1d', gdept_1d, ktype = jp_r8 ) ! reference 1D-coordinate 547 CALL iom_rstput( 0, 0, inum, 'gdepw_1d', gdepw_1d, ktype = jp_r8 ) 548 CALL iom_rstput( 0, 0, inum, 'e3t_1d' , e3t_1d , ktype = jp_r8 ) 549 CALL iom_rstput( 0, 0, inum, 'e3w_1d' , e3w_1d , ktype = jp_r8 ) 550 ! 551 CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 ) ! depth (t- & w-points) 552 CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 ) 553 ! 554 CALL iom_rstput( 0, 0, inum, 'e3t_0' , e3t_0 , ktype = jp_r8 ) ! vertical scale factors (e 555 CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8 ) 556 CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8 ) 557 CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8 ) 558 CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8 ) 559 CALL iom_rstput( 0, 0, inum, 'e3uw_0' , e3uw_0 , ktype = jp_r8 ) 560 CALL iom_rstput( 0, 0, inum, 'e3vw_0' , e3vw_0 , ktype = jp_r8 ) 561 ! 562 ! !== ocean top and bottom level ==! 563 ! 564 CALL iom_rstput( 0, 0, inum, 'bottom level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 ) ! nb of ocean T-points 565 CALL iom_rstput( 0, 0, inum, 'top level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 ) ! nb of ocean T-points (ISF) 566 ! 567 IF( ln_sco ) THEN ! s-coordinate: store grid stiffness ratio (Not required anyway) 568 CALL dom_stiff 569 !SF CALL dom_stiff( z2d ) !commented because at compilation error: The 570 !number of actual arguments cannot be greater than the number of dummy 571 !arguments. [DOM_STIFF] 572 !CALL dom_stiff( z2d ) 573 ! --------------^ compilation aborted for 574 !/workgpfs/rech/gzi/rgzi011/commit-simplif2-TOOLS/NEMOGCM/CONFIG/TEST/BLD/ppsrc/nemo/domain.f90 575 CALL iom_rstput( 0, 0, inum, 'stiffness', z2d ) ! ! Max. grid stiffness ratio 576 ENDIF 577 ! 578 ! ! ============================ 579 ! ! close the files 580 ! ! ============================ 581 CALL iom_close( inum ) 582 ! 583 END SUBROUTINE cfg_write 584 585 586 467 587 !!====================================================================== 468 588 END MODULE domain -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
r6140 r6827 16 16 USE lib_mpp ! distributed memory computing library 17 17 USE timing ! Timing 18 USE c1d ! 1D configuration19 USE domc1d ! 1D configuration: column location20 18 21 19 IMPLICIT NONE … … 82 80 !!---------------------------------------------------------------------- 83 81 ! ! recalculate jpizoom/jpjzoom given lat/lon 84 IF( lk_c1d .AND. ln_c1d_locpt ) CALL dom_c1d( rn_lat1d, rn_lon1d )85 82 ! 86 83 ! ! ============== ! -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r6140 r6827 377 377 CASE ( 0, 1, 4 ) ! mesh on the sphere 378 378 ! 379 ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 379 ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 380 ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) ) ! - - - at t-point 380 381 ! 381 382 CASE ( 2 ) ! f-plane at ppgphi0 382 383 ! 383 ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 384 ! 385 IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff(1,1) 384 ff_f(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 385 ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 386 ! 387 IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff_f(1,1) 386 388 ! 387 389 CASE ( 3 ) ! beta-plane … … 399 401 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 400 402 ! 401 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 403 ff_f(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 404 ff_t(:,:) = ( zf0 + zbeta * gphit(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 402 405 ! 403 406 IF(lwp) THEN 404 407 WRITE(numout,*) 405 WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff (nldi,nldj)406 WRITE(numout,*) ' Coriolis parameter varies from ', ff (nldi,nldj),' to ', ff(nldi,nlej)408 WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff_f(nldi,nldj) 409 WRITE(numout,*) ' Coriolis parameter varies from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej) 407 410 ENDIF 408 411 IF( lk_mpp ) THEN 409 zminff=ff (nldi,nldj)410 zmaxff=ff (nldi,nlej)412 zminff=ff_f(nldi,nldj) 413 zmaxff=ff_f(nldi,nlej) 411 414 CALL mpp_min( zminff ) ! min over the global domain 412 415 CALL mpp_max( zmaxff ) ! max over the global domain … … 420 423 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 421 424 ! 422 ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 425 ff_f(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 426 ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 423 427 ! 424 428 IF(lwp) THEN 425 429 WRITE(numout,*) 426 430 WRITE(numout,*) ' Beta-plane and rotated domain : ' 427 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff (nldi,nldj),' to ', ff(nldi,nlej)431 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff_f(nldi,nldj),' to ', ff_f(nldi,nlej) 428 432 ENDIF 429 433 ! 430 434 IF( lk_mpp ) THEN 431 zminff=ff (nldi,nldj)432 zmaxff=ff (nldi,nlej)435 zminff=ff_f(nldi,nldj) 436 zmaxff=ff_f(nldi,nlej) 433 437 CALL mpp_min( zminff ) ! min over the global domain 434 438 CALL mpp_max( zmaxff ) ! max over the global domain -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6351 r6827 22 22 USE phycst ! physical constant 23 23 USE dom_oce ! ocean space and time domain 24 USE sbc_oce ! ocean surface boundary condition25 USE wet_dry ! wetting and drying26 USE restart ! ocean restart27 24 ! 28 25 USE in_out_manager ! I/O manager … … 37 34 38 35 PUBLIC dom_vvl_init ! called by domain.F90 39 PUBLIC dom_vvl_sf_nxt ! called by step.F9040 PUBLIC dom_vvl_sf_swp ! called by step.F9041 PUBLIC dom_vvl_interpol ! called by dynnxt.F9042 36 43 37 ! !!* Namelist nam_vvl … … 133 127 ! 134 128 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 135 CALL dom_vvl_rst( nit000, 'READ' )136 129 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 137 130 ! … … 246 239 247 240 248 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )249 !!----------------------------------------------------------------------250 !! *** ROUTINE dom_vvl_sf_nxt ***251 !!252 !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,253 !! tranxt and dynspg routines254 !!255 !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.256 !! - z_tilde_case: after scale factor increment =257 !! high frequency part of horizontal divergence258 !! + retsoring towards the background grid259 !! + thickness difusion260 !! Then repartition of ssh INCREMENT proportionnaly261 !! to the "baroclinic" level thickness.262 !!263 !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case264 !! - tilde_e3t_a: after increment of vertical scale factor265 !! in z_tilde case266 !! - e3(t/u/v)_a267 !!268 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.269 !!----------------------------------------------------------------------270 INTEGER, INTENT( in ) :: kt ! time step271 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence272 !273 INTEGER :: ji, jj, jk ! dummy loop indices274 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers275 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars276 LOGICAL :: ll_do_bclinic ! local logical277 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t278 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv279 !!----------------------------------------------------------------------280 !281 IF( ln_linssh ) RETURN ! No calculation in linear free surface282 !283 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt')284 !285 CALL wrk_alloc( jpi,jpj,zht, z_scale, zwu, zwv, zhdiv )286 CALL wrk_alloc( jpi,jpj,jpk, ze3t )287 288 IF( kt == nit000 ) THEN289 IF(lwp) WRITE(numout,*)290 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'291 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'292 ENDIF293 294 ll_do_bclinic = .TRUE.295 IF( PRESENT(kcall) ) THEN296 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.297 ENDIF298 299 ! ******************************* !300 ! After acale factors at t-points !301 ! ******************************* !302 ! ! --------------------------------------------- !303 ! ! z_star coordinate and barotropic z-tilde part !304 ! ! --------------------------------------------- !305 !306 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )307 DO jk = 1, jpkm1308 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)309 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)310 END DO311 !312 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !313 ! ! ------baroclinic part------ !314 ! I - initialization315 ! ==================316 317 ! 1 - barotropic divergence318 ! -------------------------319 zhdiv(:,:) = 0._wp320 zht(:,:) = 0._wp321 DO jk = 1, jpkm1322 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)323 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)324 END DO325 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )326 327 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only)328 ! --------------------------------------------------329 IF( ln_vvl_ztilde ) THEN330 IF( kt > nit000 ) THEN331 DO jk = 1, jpkm1332 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) &333 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )334 END DO335 ENDIF336 ENDIF337 338 ! II - after z_tilde increments of vertical scale factors339 ! =======================================================340 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms341 342 ! 1 - High frequency divergence term343 ! ----------------------------------344 IF( ln_vvl_ztilde ) THEN ! z_tilde case345 DO jk = 1, jpkm1346 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )347 END DO348 ELSE ! layer case349 DO jk = 1, jpkm1350 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)351 END DO352 ENDIF353 354 ! 2 - Restoring term (z-tilde case only)355 ! ------------------356 IF( ln_vvl_ztilde ) THEN357 DO jk = 1, jpk358 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)359 END DO360 ENDIF361 362 ! 3 - Thickness diffusion term363 ! ----------------------------364 zwu(:,:) = 0._wp365 zwv(:,:) = 0._wp366 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes367 DO jj = 1, jpjm1368 DO ji = 1, fs_jpim1 ! vector opt.369 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) &370 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )371 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) &372 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )373 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)374 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)375 END DO376 END DO377 END DO378 DO jj = 1, jpj ! b - correction for last oceanic u-v points379 DO ji = 1, jpi380 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)381 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)382 END DO383 END DO384 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes385 DO jj = 2, jpjm1386 DO ji = fs_2, fs_jpim1 ! vector opt.387 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &388 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) &389 & ) * r1_e1e2t(ji,jj)390 END DO391 END DO392 END DO393 ! ! d - thickness diffusion transport: boundary conditions394 ! (stored for tracer advction and continuity equation)395 CALL lbc_lnk( un_td , 'U' , -1._wp)396 CALL lbc_lnk( vn_td , 'V' , -1._wp)397 398 ! 4 - Time stepping of baroclinic scale factors399 ! ---------------------------------------------400 ! Leapfrog time stepping401 ! ~~~~~~~~~~~~~~~~~~~~~~402 IF( neuler == 0 .AND. kt == nit000 ) THEN403 z2dt = rdt404 ELSE405 z2dt = 2.0_wp * rdt406 ENDIF407 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )408 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)409 410 ! Maximum deformation control411 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~412 ze3t(:,:,jpk) = 0._wp413 DO jk = 1, jpkm1414 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)415 END DO416 z_tmax = MAXVAL( ze3t(:,:,:) )417 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain418 z_tmin = MINVAL( ze3t(:,:,:) )419 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain420 ! - ML - test: for the moment, stop simulation for too large e3_t variations421 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN422 IF( lk_mpp ) THEN423 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )424 CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )425 ELSE426 ijk_max = MAXLOC( ze3t(:,:,:) )427 ijk_max(1) = ijk_max(1) + nimpp - 1428 ijk_max(2) = ijk_max(2) + njmpp - 1429 ijk_min = MINLOC( ze3t(:,:,:) )430 ijk_min(1) = ijk_min(1) + nimpp - 1431 ijk_min(2) = ijk_min(2) + njmpp - 1432 ENDIF433 IF (lwp) THEN434 WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax435 WRITE(numout, *) 'at i, j, k=', ijk_max436 WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin437 WRITE(numout, *) 'at i, j, k=', ijk_min438 CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')439 ENDIF440 ENDIF441 ! - ML - end test442 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below443 tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )444 tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )445 446 !447 ! "tilda" change in the after scale factor448 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~449 DO jk = 1, jpkm1450 dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)451 END DO452 ! III - Barotropic repartition of the sea surface height over the baroclinic profile453 ! ==================================================================================454 ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)455 ! - ML - baroclinicity error should be better treated in the future456 ! i.e. locally and not spread over the water column.457 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible)458 zht(:,:) = 0.459 DO jk = 1, jpkm1460 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)461 END DO462 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )463 DO jk = 1, jpkm1464 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)465 END DO466 467 ENDIF468 469 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !470 ! ! ---baroclinic part--------- !471 DO jk = 1, jpkm1472 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)473 END DO474 ENDIF475 476 IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging477 !478 IF( lwp ) WRITE(numout, *) 'kt =', kt479 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN480 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )481 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain482 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax483 END IF484 !485 zht(:,:) = 0.0_wp486 DO jk = 1, jpkm1487 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)488 END DO489 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )490 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain491 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax492 !493 zht(:,:) = 0.0_wp494 DO jk = 1, jpkm1495 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)496 END DO497 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )498 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain499 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax500 !501 zht(:,:) = 0.0_wp502 DO jk = 1, jpkm1503 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)504 END DO505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )506 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax508 !509 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) )510 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain511 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax512 !513 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) )514 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain515 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax516 !517 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) )518 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain519 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax520 END IF521 522 ! *********************************** !523 ! After scale factors at u- v- points !524 ! *********************************** !525 526 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )527 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )528 529 ! *********************************** !530 ! After depths at u- v points !531 ! *********************************** !532 533 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)534 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)535 DO jk = 2, jpkm1536 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)537 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)538 END DO539 ! ! Inverse of the local depth540 !!gm BUG ? don't understand the use of umask_i here .....541 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )542 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )543 !544 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv )545 CALL wrk_dealloc( jpi,jpj,jpk, ze3t )546 !547 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt')548 !549 END SUBROUTINE dom_vvl_sf_nxt550 551 552 SUBROUTINE dom_vvl_sf_swp( kt )553 !!----------------------------------------------------------------------554 !! *** ROUTINE dom_vvl_sf_swp ***555 !!556 !! ** Purpose : compute time filter and swap of scale factors557 !! compute all depths and related variables for next time step558 !! write outputs and restart file559 !!560 !! ** Method : - swap of e3t with trick for volume/tracer conservation561 !! - reconstruct scale factor at other grid points (interpolate)562 !! - recompute depths and water height fields563 !!564 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step565 !! - Recompute:566 !! e3(u/v)_b567 !! e3w_n568 !! e3(u/v)w_b569 !! e3(u/v)w_n570 !! gdept_n, gdepw_n and gde3w_n571 !! h(u/v) and h(u/v)r572 !!573 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.574 !! Leclair, M., and G. Madec, 2011, Ocean Modelling.575 !!----------------------------------------------------------------------576 INTEGER, INTENT( in ) :: kt ! time step577 !578 INTEGER :: ji, jj, jk ! dummy loop indices579 REAL(wp) :: zcoef ! local scalar580 !!----------------------------------------------------------------------581 !582 IF( ln_linssh ) RETURN ! No calculation in linear free surface583 !584 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp')585 !586 IF( kt == nit000 ) THEN587 IF(lwp) WRITE(numout,*)588 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'589 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'590 ENDIF591 !592 ! Time filter and swap of scale factors593 ! =====================================594 ! - ML - e3(t/u/v)_b are allready computed in dynnxt.595 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN596 IF( neuler == 0 .AND. kt == nit000 ) THEN597 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)598 ELSE599 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &600 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )601 ENDIF602 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)603 ENDIF604 gdept_b(:,:,:) = gdept_n(:,:,:)605 gdepw_b(:,:,:) = gdepw_n(:,:,:)606 607 e3t_n(:,:,:) = e3t_a(:,:,:)608 e3u_n(:,:,:) = e3u_a(:,:,:)609 e3v_n(:,:,:) = e3v_a(:,:,:)610 611 ! Compute all missing vertical scale factor and depths612 ! ====================================================613 ! Horizontal scale factor interpolations614 ! --------------------------------------615 ! - ML - e3u_b and e3v_b are allready computed in dynnxt616 ! - JC - hu_b, hv_b, hur_b, hvr_b also617 618 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )619 620 ! Vertical scale factor interpolations621 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )622 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )623 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )624 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' )625 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )626 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )627 628 ! t- and w- points depth (set the isf depth as it is in the initial step)629 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)630 gdepw_n(:,:,1) = 0.0_wp631 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)632 DO jk = 2, jpk633 DO jj = 1,jpj634 DO ji = 1,jpi635 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt636 ! 1 for jk = mikt637 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))638 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)639 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) &640 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) )641 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)642 END DO643 END DO644 END DO645 646 ! Local depth and Inverse of the local depth of the water647 ! -------------------------------------------------------648 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:)649 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:)650 !651 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)652 DO jk = 2, jpkm1653 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)654 END DO655 656 ! write restart file657 ! ==================658 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )659 !660 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')661 !662 END SUBROUTINE dom_vvl_sf_swp663 664 665 241 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 666 242 !!--------------------------------------------------------------------- … … 684 260 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 685 261 ! 686 IF(ln_wd) THEN 687 zlnwd = 1.0_wp 688 ELSE 689 zlnwd = 0.0_wp 690 END IF 262 zlnwd = 0.0_wp 691 263 ! 692 264 SELECT CASE ( pout ) !== type of interpolation ==! … … 772 344 ! 773 345 END SUBROUTINE dom_vvl_interpol 774 775 776 SUBROUTINE dom_vvl_rst( kt, cdrw )777 !!---------------------------------------------------------------------778 !! *** ROUTINE dom_vvl_rst ***779 !!780 !! ** Purpose : Read or write VVL file in restart file781 !!782 !! ** Method : use of IOM library783 !! if the restart does not contain vertical scale factors,784 !! they are set to the _0 values785 !! if the restart does not contain vertical scale factors increments (z_tilde),786 !! they are set to 0.787 !!----------------------------------------------------------------------788 INTEGER , INTENT(in) :: kt ! ocean time-step789 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag790 !791 INTEGER :: ji, jj, jk792 INTEGER :: id1, id2, id3, id4, id5 ! local integers793 !!----------------------------------------------------------------------794 !795 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst')796 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise797 ! ! ===============798 IF( ln_rstart ) THEN !* Read the restart file799 CALL rst_read_open ! open the restart file if necessary800 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn )801 !802 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )803 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )804 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )805 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )806 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )807 ! ! --------- !808 ! ! all cases !809 ! ! --------- !810 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist811 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )812 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )813 ! needed to restart if land processor not computed814 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'815 WHERE ( tmask(:,:,:) == 0.0_wp )816 e3t_n(:,:,:) = e3t_0(:,:,:)817 e3t_b(:,:,:) = e3t_0(:,:,:)818 END WHERE819 IF( neuler == 0 ) THEN820 e3t_b(:,:,:) = e3t_n(:,:,:)821 ENDIF822 ELSE IF( id1 > 0 ) THEN823 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'824 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'825 IF(lwp) write(numout,*) 'neuler is forced to 0'826 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )827 e3t_n(:,:,:) = e3t_b(:,:,:)828 neuler = 0829 ELSE IF( id2 > 0 ) THEN830 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'831 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'832 IF(lwp) write(numout,*) 'neuler is forced to 0'833 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )834 e3t_b(:,:,:) = e3t_n(:,:,:)835 neuler = 0836 ELSE837 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'838 IF(lwp) write(numout,*) 'Compute scale factor from sshn'839 IF(lwp) write(numout,*) 'neuler is forced to 0'840 DO jk = 1, jpk841 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &842 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &843 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))844 END DO845 e3t_b(:,:,:) = e3t_n(:,:,:)846 neuler = 0847 ENDIF848 ! ! ----------- !849 IF( ln_vvl_zstar ) THEN ! z_star case !850 ! ! ----------- !851 IF( MIN( id3, id4 ) > 0 ) THEN852 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )853 ENDIF854 ! ! ----------------------- !855 ELSE ! z_tilde and layer cases !856 ! ! ----------------------- !857 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist858 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )859 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )860 ELSE ! one at least array is missing861 tilde_e3t_b(:,:,:) = 0.0_wp862 tilde_e3t_n(:,:,:) = 0.0_wp863 ENDIF864 ! ! ------------ !865 IF( ln_vvl_ztilde ) THEN ! z_tilde case !866 ! ! ------------ !867 IF( id5 > 0 ) THEN ! required array exists868 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )869 ELSE ! array is missing870 hdiv_lf(:,:,:) = 0.0_wp871 ENDIF872 ENDIF873 ENDIF874 !875 ELSE !* Initialize at "rest"876 e3t_b(:,:,:) = e3t_0(:,:,:)877 e3t_n(:,:,:) = e3t_0(:,:,:)878 sshn(:,:) = 0.0_wp879 880 IF( ln_wd ) THEN881 DO jj = 1, jpj882 DO ji = 1, jpi883 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN884 e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1885 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1886 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1887 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)888 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)889 ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj)890 ENDIF891 ENDDO892 ENDDO893 END IF894 895 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN896 tilde_e3t_b(:,:,:) = 0.0_wp897 tilde_e3t_n(:,:,:) = 0.0_wp898 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp899 END IF900 ENDIF901 !902 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file903 ! ! ===================904 IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'905 ! ! --------- !906 ! ! all cases !907 ! ! --------- !908 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) )909 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) )910 ! ! ----------------------- !911 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !912 ! ! ----------------------- !913 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )914 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )915 END IF916 ! ! -------------!917 IF( ln_vvl_ztilde ) THEN ! z_tilde case !918 ! ! ------------ !919 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )920 ENDIF921 !922 ENDIF923 !924 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst')925 !926 END SUBROUTINE dom_vvl_rst927 346 928 347 … … 990 409 ! 991 410 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 992 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )993 411 ! 994 412 IF(lwp) THEN ! Print the choice -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r6689 r6827 249 249 CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 ) 250 250 251 CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 ) ! ! coriolis factor 251 CALL iom_rstput( 0, 0, inum3, 'ff_f', ff_f, ktype = jp_r8 ) ! ! coriolis factor 252 CALL iom_rstput( 0, 0, inum3, 'ff_t', ff_t, ktype = jp_r8 ) ! ! coriolis factor 252 253 253 254 ! note that mbkt is set to 1 over land ==> use surface tmask -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6492 r6827 37 37 USE oce ! ocean variables 38 38 USE dom_oce ! ocean domain 39 USE wet_dry ! wetting and drying40 39 USE closea ! closed seas 41 USE c1d ! 1D vertical configuration42 40 ! 43 41 USE in_out_manager ! I/O manager … … 147 145 CALL zgr_z ! Reference z-coordinate system (always called) 148 146 CALL zgr_bat ! Bathymetry fields (levels and meters) 149 IF( lk_c1d ) CALL lbc_lnk( bathy , 'T', 1._wp ) ! 1D config.: same bathy value over the 3x3 domain150 147 IF( ln_zco ) CALL zgr_zco ! z-coordinate 151 148 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate … … 155 152 ! ----------------------------------- 156 153 IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain 157 IF( .NOT.lk_c1d )CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points154 CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 158 155 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 159 156 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 160 !161 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain162 ibat = mbathy(2,2)163 mbathy(:,:) = ibat164 END IF165 157 ! 166 158 IF( nprint == 1 .AND. lwp ) THEN … … 1950 1942 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1951 1943 1952 IF( .NOT.ln_wd ) THEN1953 1944 DO jj = 1, jpj 1954 1945 DO ji = 1, jpi … … 1956 1947 END DO 1957 1948 END DO 1958 END IF1959 1949 ! ! ============================= 1960 1950 ! ! Define the envelop bathymetry (hbatt) … … 1963 1953 zenv(:,:) = bathy(:,:) 1964 1954 ! 1965 IF( .NOT.ln_wd ) THEN1966 1955 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1967 1956 DO jj = 1, jpj … … 1986 1975 END DO 1987 1976 END DO 1988 END IF1989 1977 1990 1978 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero … … 2080 2068 IF(lwp) THEN 2081 2069 WRITE(numout,*) 2082 IF( .NOT.ln_wd ) THEN2083 2070 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2084 ELSE2085 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min2086 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld2087 ENDIF2088 2071 ENDIF 2089 2072 hbatu(:,:) = rn_sbot_min … … 2099 2082 END DO 2100 2083 2101 IF( ln_wd ) THEN !avoid the zero depth on T- (U-,V-,F-) points2102 DO jj = 1, jpj2103 DO ji = 1, jpi2104 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) &2105 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin12106 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) &2107 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin12108 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) &2109 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin12110 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) &2111 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin12112 END DO2113 END DO2114 END IF2115 2084 ! 2116 2085 ! Apply lateral boundary condition … … 2146 2115 2147 2116 !!bug: key_helsinki a verifer 2148 IF( .NOT.ln_wd ) THEN2149 2117 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2150 2118 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2151 2119 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2152 2120 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2153 END IF2154 2121 2155 2122 IF( nprint == 1 .AND. lwp ) THEN … … 2193 2160 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2194 2161 ! 2195 IF( .NOT.ln_wd ) THEN2196 2162 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2197 2163 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp … … 2201 2167 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2202 2168 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2203 END IF2204 2169 2205 2170 #if defined key_agrif … … 2234 2199 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2235 2200 END DO 2236 IF( ln_wd ) THEN2237 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN2238 mbathy(ji,jj) = 02239 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN2240 mbathy(ji,jj) = 22241 ENDIF2242 ELSE2243 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 02244 ENDIF2245 2201 END DO 2246 2202 END DO … … 2425 2381 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2426 2382 DO jk = 1, jpk 2427 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN2428 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )2429 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )2430 ELSE2431 2383 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2432 2384 & / ztmpu 2433 2385 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2434 2386 & / ztmpu 2435 END IF 2436 2437 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2438 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2439 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2440 ELSE 2387 2441 2388 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2442 2389 & / ztmpv 2443 2390 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2444 2391 & / ztmpv 2445 END IF 2446 2447 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2448 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj ,jk) + z_esigt3(ji+1,jj ,jk) & 2449 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2450 ELSE 2392 2451 2393 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2452 2394 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2453 2395 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2454 2396 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2455 END IF2456 2397 2457 2398 ! … … 2582 2523 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2583 2524 DO jk = 1, jpk 2584 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN2585 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )2586 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )2587 ELSE2588 2525 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2589 2526 & / ztmpu 2590 2527 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2591 2528 & / ztmpu 2592 END IF2593 2594 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN2595 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )2596 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )2597 ELSE2598 2529 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2599 2530 & / ztmpv 2600 2531 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2601 2532 & / ztmpv 2602 END IF 2603 2604 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2605 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2606 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2607 ELSE 2533 2608 2534 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2609 2535 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2610 2536 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2611 2537 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2612 END IF2613 2538 2614 2539 ! Code prior to wetting and drying option (for reference) -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
r6519 r6827 20 20 !!-------------------------------------------------------------------- 21 21 USE dom_oce ! ocean space and time domain 22 USE c1d ! 1D vertical configuration23 USE flo_oce ! floats module declarations24 22 USE lbclnk ! lateal boundary condition / mpp exchanges 25 23 USE iom_def ! iom variables definitions 26 24 USE iom_nf90 ! NetCDF format with native NetCDF library 27 25 USE in_out_manager ! I/O manager 28 USE lib_mpp ! MPP library 29 #if defined key_iomput 30 USE sbc_oce, ONLY : nn_fsbc ! ocean space and time domain 31 USE trc_oce, ONLY : nn_dttrc ! !: frequency of step on passive tracers 32 USE icb_oce, ONLY : nclasses, class_num ! !: iceberg classes 33 #if defined key_lim3 34 USE ice , ONLY : jpl 35 #elif defined key_lim2 36 USE par_ice_2 37 #endif 26 USE lib_mpp ! MPP library 38 27 USE domngb ! ocean space and time domain 39 28 USE phycst ! physical constants 40 USE dianam ! build name of file29 !SF USE dianam ! build name of file 41 30 USE xios 42 # endif43 31 USE ioipsl, ONLY : ju2ymds ! for calendar 44 USE crs ! Grid coarsening45 32 46 33 IMPLICIT NONE … … 139 126 ENDIF 140 127 141 IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN142 CALL dom_grid_crs ! Save the parent grid information & Switch to coarse grid domain143 !144 CALL set_grid( "T", glamt_crs, gphit_crs )145 CALL set_grid( "U", glamu_crs, gphiu_crs )146 CALL set_grid( "V", glamv_crs, gphiv_crs )147 CALL set_grid( "W", glamt_crs, gphit_crs )148 CALL set_grid_znl( gphit_crs )149 !150 CALL dom_grid_glo ! Return to parent grid domain151 !152 IF( ln_cfmeta ) THEN ! Add additional grid metadata153 CALL iom_set_domain_attr("grid_T", area = e1e2t_crs(nldi:nlei, nldj:nlej))154 CALL iom_set_domain_attr("grid_U", area = e1u_crs(nldi:nlei, nldj:nlej) * e2u_crs(nldi:nlei, nldj:nlej))155 CALL iom_set_domain_attr("grid_V", area = e1v_crs(nldi:nlei, nldj:nlej) * e2v_crs(nldi:nlei, nldj:nlej))156 CALL iom_set_domain_attr("grid_W", area = e1e2t_crs(nldi:nlei, nldj:nlej))157 CALL set_grid_bounds( "T", glamf_crs, gphif_crs, glamt_crs, gphit_crs )158 CALL set_grid_bounds( "U", glamv_crs, gphiv_crs, glamu_crs, gphiu_crs )159 CALL set_grid_bounds( "V", glamu_crs, gphiu_crs, glamv_crs, gphiv_crs )160 CALL set_grid_bounds( "W", glamf_crs, gphif_crs, glamt_crs, gphit_crs )161 ENDIF162 ENDIF163 164 128 ! vertical grid definition 165 129 CALL iom_set_axis_attr( "deptht", gdept_1d ) … … 180 144 CALL iom_set_axis_attr( "depthw", bounds=z_bnds ) 181 145 182 # if defined key_floats183 CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,nfloat) /) )184 # endif185 #if defined key_lim3 || defined key_lim2186 CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) )187 #endif188 CALL iom_set_axis_attr( "icbcla", class_num )189 146 CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) ) 190 147 CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) ) … … 882 839 ENDIF 883 840 884 ! C1D case : always call lbc_lnk to replicate the central value over the whole 3X3 domain885 IF( lk_c1d .AND. PRESENT(pv_r2d) ) CALL lbc_lnk( pv_r2d,'Z',1. )886 IF( lk_c1d .AND. PRESENT(pv_r3d) ) CALL lbc_lnk( pv_r3d,'Z',1. )887 888 841 !--- Apply scale_factor and offset 889 842 zscf = iom_file(kiomid)%scf(idvar) ! scale factor … … 1459 1412 ! frequency of the call of iom_put (attribut: freq_op) 1460 1413 WRITE(cl1,'(i1)') 1 ; CALL iom_set_field_attr('field_definition', freq_op = cl1//'ts', freq_offset='0ts') 1461 WRITE(cl1,'(i1)') nn_fsbc ; CALL iom_set_field_attr('SBC' , freq_op = cl1//'ts', freq_offset='0ts')1462 WRITE(cl1,'(i1)') nn_fsbc ; CALL iom_set_field_attr('SBC_scalar' , freq_op = cl1//'ts', freq_offset='0ts')1463 WRITE(cl1,'(i1)') nn_dttrc ; CALL iom_set_field_attr('ptrc_T' , freq_op = cl1//'ts', freq_offset='0ts')1464 WRITE(cl1,'(i1)') nn_dttrc ; CALL iom_set_field_attr('diad_T' , freq_op = cl1//'ts', freq_offset='0ts')1465 1414 1466 1415 ! output file names (attribut: name) -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r6505 r6827 39 39 USE dom_oce ! ocean space and time domain 40 40 USE phycst ! physical constants 41 USE stopar ! Stochastic T/S fluctuations42 USE stopts ! Stochastic T/S fluctuations43 41 ! 44 42 USE in_out_manager ! I/O manager … … 335 333 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 336 334 ! 337 ! Stochastic equation of state338 IF ( ln_sto_eos ) THEN339 ALLOCATE(zn0_sto(1:2*nn_sto_eos))340 ALLOCATE(zn_sto(1:2*nn_sto_eos))341 ALLOCATE(zsign(1:2*nn_sto_eos))342 DO jsmp = 1, 2*nn_sto_eos, 2343 zsign(jsmp) = 1._wp344 zsign(jsmp+1) = -1._wp345 END DO346 !347 DO jk = 1, jpkm1348 DO jj = 1, jpj349 DO ji = 1, jpi350 !351 ! compute density (2*nn_sto_eos) times:352 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)353 ! (2) for t-dt, s-ds (with the opposite fluctuation)354 DO jsmp = 1, nn_sto_eos*2355 jdof = (jsmp + 1) / 2356 zh = pdep(ji,jj,jk) * r1_Z0 ! depth357 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature358 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)359 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity360 ztm = tmask(ji,jj,jk) ! tmask361 !362 zn3 = EOS013*zt &363 & + EOS103*zs+EOS003364 !365 zn2 = (EOS022*zt &366 & + EOS112*zs+EOS012)*zt &367 & + (EOS202*zs+EOS102)*zs+EOS002368 !369 zn1 = (((EOS041*zt &370 & + EOS131*zs+EOS031)*zt &371 & + (EOS221*zs+EOS121)*zs+EOS021)*zt &372 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &373 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001374 !375 zn0_sto(jsmp) = (((((EOS060*zt &376 & + EOS150*zs+EOS050)*zt &377 & + (EOS240*zs+EOS140)*zs+EOS040)*zt &378 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &379 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &380 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &381 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000382 !383 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)384 END DO385 !386 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities387 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp388 DO jsmp = 1, nn_sto_eos*2389 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface390 !391 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked)392 END DO393 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos394 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos395 END DO396 END DO397 END DO398 DEALLOCATE(zn0_sto,zn_sto,zsign)399 ! Non-stochastic equation of state400 ELSE401 335 DO jk = 1, jpkm1 402 336 DO jj = 1, jpj … … 437 371 END DO 438 372 END DO 439 ENDIF440 373 441 374 CASE( np_seos ) !== simplified EOS ==! -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r6152 r6827 48 48 USE mppini ! shared/distributed memory setting (mpp_init routine) 49 49 USE domain ! domain initialization (dom_init routine) 50 #if defined key_nemocice_decomp51 USE ice_domain_size, only: nx_global, ny_global52 #endif53 USE tideini ! tidal components initialization (tide_ini routine)54 USE bdyini ! open boundary cond. setting (bdy_init routine)55 USE bdydta ! open boundary cond. setting (bdy_dta_init routine)56 USE bdytides ! open boundary cond. setting (bdytide_init routine)57 USE sbctide, ONLY : lk_tide58 USE istate ! initial state setting (istate_init routine)59 USE ldfdyn ! lateral viscosity setting (ldfdyn_init routine)60 USE ldftra ! lateral diffusivity setting (ldftra_init routine)61 USE zdfini ! vertical physics setting (zdf_init routine)62 50 USE phycst ! physical constant (par_cst routine) 63 USE trdini ! dyn/tra trends initialization (trd_init routine)64 USE asminc ! assimilation increments65 USE asmbkg ! writing out state trajectory66 USE diaptr ! poleward transports (dia_ptr_init routine)67 USE diadct ! sections transports (dia_dct_init routine)68 USE diaobs ! Observation diagnostics (dia_obs_init routine)69 USE diacfl ! CFL diagnostics (dia_cfl_init routine)70 51 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 71 USE step ! NEMO time-stepping (stp routine)72 USE icbini ! handle bergs, initialisation73 USE icbstp ! handle bergs, calving, themodynamics and transport74 USE cpl_oasis3 ! OASIS3 coupling75 USE c1d ! 1D configuration76 USE step_c1d ! Time stepping loop for the 1D configuration77 USE dyndmp ! Momentum damping78 USE stopar ! Stochastic param.: ???79 USE stopts ! Stochastic param.: ???80 #if defined key_top81 USE trcini ! passive tracer initialisation82 #endif83 52 USE lib_mpp ! distributed memory computing 84 USE diurnal_bulk ! diurnal bulk SST85 USE step_diu ! diurnal bulk SST timestepping (called from here if run offline)86 53 #if defined key_iomput 87 54 USE xios ! xIOserver 88 55 #endif 89 USE crsini ! initialise grid coarsening utility90 56 USE lbcnfd , ONLY : isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 91 USE sbc_oce, ONLY : lk_oasis92 USE diatmb ! Top,middle,bottom output93 USE dia25h ! 25h mean output94 USE wet_dry ! Wetting and drying setting (wad_init routine)95 57 96 58 IMPLICIT NONE … … 137 99 CALL Agrif_Declare_Var_dom ! AGRIF: set the meshes for DOM 138 100 CALL Agrif_Declare_Var ! " " " " " DYN/TRA 139 # if defined key_top140 CALL Agrif_Declare_Var_top ! " " " " " TOP141 # endif142 # if defined key_lim2143 CALL Agrif_Declare_Var_lim2 ! " " " " " LIM144 # endif145 101 #endif 146 102 ! check that all process are still there... If some process have an error, … … 151 107 152 108 ! !-----------------------! 153 ! !== time stepping ==!154 ! !-----------------------!155 istp = nit000156 #if defined key_c1d157 DO WHILE ( istp <= nitend .AND. nstop == 0 )158 CALL stp_c1d( istp )159 istp = istp + 1160 END DO161 #else162 IF( lk_asminc ) THEN163 IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields164 IF( ln_asmdin ) THEN ! Direct initialization165 IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 ) ! Tracers166 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 ) ! Dynamics167 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 ) ! SSH168 ENDIF169 ENDIF170 171 #if defined key_agrif172 CALL Agrif_Regrid()173 #endif174 175 DO WHILE ( istp <= nitend .AND. nstop == 0 )176 #if defined key_agrif177 CALL stp ! AGRIF: time stepping178 #else179 IF ( .NOT. ln_diurnal_only ) THEN180 CALL stp( istp ) ! standard time stepping181 ELSE182 CALL stp_diurnal( istp ) ! time step only the diurnal SST183 ENDIF184 #endif185 istp = istp + 1186 IF( lk_mpp ) CALL mpp_max( nstop )187 END DO188 #endif189 190 IF( ln_diaobs ) CALL dia_obs_wri191 !192 IF( ln_icebergs ) CALL icb_end( nitend )193 194 ! !------------------------!195 109 ! !== finalize the run ==! 196 110 ! !------------------------! 197 IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA198 111 ! 199 112 IF( nstop /= 0 .AND. lwp ) THEN ! error print … … 204 117 #if defined key_agrif 205 118 IF( .NOT. Agrif_Root() ) THEN 206 CALL Agrif_ParentGrid_To_ChildGrid() 207 IF( ln_diaobs ) CALL dia_obs_wri 119 CALL Agrif_ParentGrid_To_ChildGrid() 208 120 IF( nn_timing == 1 ) CALL timing_finalize 209 121 CALL Agrif_ChildGrid_To_ParentGrid() … … 216 128 #if defined key_iomput 217 129 CALL xios_finalize ! end mpp communications with xios 218 IF( lk_oasis ) CALL cpl_finalize ! end coupling and mpp communications with OASIS219 130 #else 220 IF( lk_oasis ) THEN 221 CALL cpl_finalize ! end coupling and mpp communications with OASIS 222 ELSE 223 IF( lk_mpp ) CALL mppstop ! end mpp communications 224 ENDIF 131 IF( lk_mpp ) CALL mppstop ! end mpp communications 225 132 #endif 226 133 ! … … 294 201 #if defined key_iomput 295 202 IF( Agrif_Root() ) THEN 296 IF( lk_oasis ) THEN297 CALL cpl_init( "oceanx", ilocal_comm ) ! nemo local communicator given by oasis298 CALL xios_initialize( "not used",local_comm=ilocal_comm ) ! send nemo communicator to xios299 ELSE300 203 CALL xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm ) ! nemo local communicator given by xios 301 ENDIF302 204 ENDIF 303 205 ! Nodes selection (control print return in cltxt) 304 206 narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm ) 305 207 #else 306 IF( lk_oasis ) THEN307 IF( Agrif_Root() ) THEN308 CALL cpl_init( "oceanx", ilocal_comm ) ! nemo local communicator given by oasis309 ENDIF310 ! Nodes selection (control print return in cltxt)311 narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )312 ELSE313 208 ilocal_comm = 0 314 209 ! Nodes selection (control print return in cltxt) 315 210 narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop ) 316 ENDIF317 211 #endif 318 212 narea = narea + 1 ! mynode return the rank of proc (0 --> jpnij -1 ) … … 402 296 CALL phy_cst ! Physical constants 403 297 CALL eos_init ! Equation of state 404 IF( lk_c1d ) CALL c1d_init ! 1D column configuration405 CALL wad_init ! Wetting and drying options406 298 CALL dom_cfg ! Domain configuration 407 299 CALL dom_init ! Domain 408 IF( ln_crs ) CALL crs_init ! coarsened grid: domain initialization409 IF( ln_nnogather ) CALL nemo_northcomms! northfold neighbour lists (must be done after the masks are defined)410 300 IF( ln_ctl ) CALL prt_ctl_init ! Print control 411 412 CALL diurnal_sst_bulk_init ! diurnal sst413 IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init ! cool skin414 415 ! IF ln_diurnal_only, then we only want a subset of the initialisation routines416 IF ( ln_diurnal_only ) THEN417 CALL istate_init ! ocean initial state (Dynamics and tracers)418 CALL sbc_init ! Forcings : surface module419 CALL tra_qsr_init ! penetrative solar radiation qsr420 IF( ln_diaobs ) THEN ! Observation & model comparison421 CALL dia_obs_init ! Initialize observational data422 CALL dia_obs( nit000 - 1 ) ! Observation operator for restart423 ENDIF424 ! ! Assimilation increments425 IF( lk_asminc ) CALL asm_inc_init ! Initialize assimilation increments426 427 IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler428 RETURN429 ENDIF430 431 CALL istate_init ! ocean initial state (Dynamics and tracers)432 433 ! ! external forcing434 !!gm to be added : creation and call of sbc_apr_init435 IF( lk_tide ) CALL tide_init ! tidal harmonics436 CALL sbc_init ! surface boundary conditions (including sea-ice)437 !!gm ==>> bdy_init should call bdy_dta_init and bdytide_init NOT in nemogcm !!!438 IF( lk_bdy ) CALL bdy_init ! Open boundaries initialisation439 IF( lk_bdy ) CALL bdy_dta_init ! Open boundaries initialisation of external data arrays440 IF( lk_bdy .AND. lk_tide ) &441 & CALL bdytide_init ! Open boundaries initialisation of tidal harmonic forcing442 443 ! ! Ocean physics444 ! ! Vertical physics445 CALL zdf_init ! namelist read446 CALL zdf_bfr_init ! bottom friction447 IF( lk_zdfric ) CALL zdf_ric_init ! Richardson number dependent Kz448 IF( lk_zdftke ) CALL zdf_tke_init ! TKE closure scheme449 IF( lk_zdfgls ) CALL zdf_gls_init ! GLS closure scheme450 IF( lk_zdftmx ) CALL zdf_tmx_init ! tidal vertical mixing451 IF( lk_zdfddm ) CALL zdf_ddm_init ! double diffusive mixing452 453 ! ! Lateral physics454 CALL ldf_tra_init ! Lateral ocean tracer physics455 CALL ldf_eiv_init ! eddy induced velocity param.456 CALL ldf_dyn_init ! Lateral ocean momentum physics457 458 ! ! Active tracers459 CALL tra_qsr_init ! penetrative solar radiation qsr460 CALL tra_bbc_init ! bottom heat flux461 IF( lk_trabbl ) CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme462 CALL tra_dmp_init ! internal tracer damping463 CALL tra_adv_init ! horizontal & vertical advection464 CALL tra_ldf_init ! lateral mixing465 CALL tra_zdf_init ! vertical mixing and after tracer fields466 467 ! ! Dynamics468 IF( lk_c1d ) CALL dyn_dmp_init ! internal momentum damping469 CALL dyn_adv_init ! advection (vector or flux form)470 CALL dyn_vor_init ! vorticity term including Coriolis471 CALL dyn_ldf_init ! lateral mixing472 CALL dyn_hpg_init ! horizontal gradient of Hydrostatic pressure473 CALL dyn_zdf_init ! vertical diffusion474 CALL dyn_spg_init ! surface pressure gradient475 476 #if defined key_top477 ! ! Passive tracers478 CALL trc_init479 #endif480 IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing481 482 ! ! Icebergs483 CALL icb_init( rdt, nit000) ! initialise icebergs instance484 485 ! ! Misc. options486 CALL sto_par_init ! Stochastic parametrization487 IF( ln_sto_eos ) CALL sto_pts_init ! RRandom T/S fluctuations488 489 ! ! Diagnostics490 IF( lk_floats ) CALL flo_init ! drifting Floats491 CALL dia_cfl_init ! Initialise CFL diagnostics492 IF( lk_diaar5 ) CALL dia_ar5_init ! ar5 diag493 CALL dia_ptr_init ! Poleward TRansports initialization494 IF( lk_diadct ) CALL dia_dct_init ! Sections tranports495 CALL dia_hsb_init ! heat content, salt content and volume budgets496 CALL trd_init ! Mixed-layer/Vorticity/Integral constraints trends497 CALL dia_obs_init ! Initialize observational data498 IF( ln_diaobs ) CALL dia_obs( nit000 - 1 ) ! Observation operator for restart499 500 ! ! Assimilation increments501 IF( lk_asminc ) CALL asm_inc_init ! Initialize assimilation increments502 IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler503 CALL dia_tmb_init ! TMB outputs504 CALL dia_25h_init ! 25h mean outputs505 506 301 ! 507 302 END SUBROUTINE nemo_init … … 600 395 ENDIF 601 396 ! 602 IF( nbench == 1 ) THEN ! Benchmark603 SELECT CASE ( cp_cfg )604 CASE ( 'gyre' ) ; CALL ctl_warn( ' The Benchmark is activated ' )605 CASE DEFAULT ; CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:', &606 & ' cp_cfg = "gyre" in namelist &namcfg or set nbench = 0' )607 END SELECT608 ENDIF609 !610 397 IF( 1_wp /= SIGN(1._wp,-0._wp) ) CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows ', & 611 398 & 'f2003 standard. ' , & … … 653 440 !! ** Method : 654 441 !!---------------------------------------------------------------------- 655 USE diawri , ONLY: dia_wri_alloc656 442 USE dom_oce , ONLY: dom_oce_alloc 657 USE trc_oce , ONLY: trc_oce_alloc658 #if defined key_diadct659 USE diadct , ONLY: diadct_alloc660 #endif661 #if defined key_bdy662 USE bdy_oce , ONLY: bdy_oce_alloc663 #endif664 443 ! 665 444 INTEGER :: ierr … … 667 446 ! 668 447 ierr = oce_alloc () ! ocean 669 ierr = ierr + dia_wri_alloc ()670 448 ierr = ierr + dom_oce_alloc () ! ocean domain 671 ierr = ierr + zdf_oce_alloc () ! ocean vertical physics672 !673 ierr = ierr + trc_oce_alloc () ! shared TRC / TRA arrays674 !675 #if defined key_diadct676 ierr = ierr + diadct_alloc () !677 #endif678 #if defined key_bdy679 ierr = ierr + bdy_oce_alloc () ! bdy masks (incl. initialization)680 #endif681 449 ! 682 450 IF( lk_mpp ) CALL mpp_sum( ierr ) -
branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r6140 r6827 9 9 USE oce ! ocean dynamics and tracers variables 10 10 USE dom_oce ! ocean space and time domain variables 11 USE zdf_oce ! ocean vertical physics variables12 11 13 12 USE daymod ! calendar (day routine) 14 13 15 USE sbc_oce ! surface boundary condition: ocean16 USE sbcmod ! surface boundary condition (sbc routine)17 USE sbcrnf ! surface boundary condition: runoff variables18 USE sbccpl ! surface boundary condition: coupled formulation (call send at end of step)19 USE sbcapr ! surface boundary condition: atmospheric pressure20 USE sbctide ! Tide initialisation21 22 USE traqsr ! solar radiation penetration (tra_qsr routine)23 USE trasbc ! surface boundary condition (tra_sbc routine)24 USE trabbc ! bottom boundary condition (tra_bbc routine)25 USE trabbl ! bottom boundary layer (tra_bbl routine)26 USE tradmp ! internal damping (tra_dmp routine)27 USE traadv ! advection scheme control (tra_adv_ctl routine)28 USE traldf ! lateral mixing (tra_ldf routine)29 USE trazdf ! vertical mixing (tra_zdf routine)30 USE tranxt ! time-stepping (tra_nxt routine)31 USE tranpc ! non-penetrative convection (tra_npc routine)32 33 14 USE eosbn2 ! equation of state (eos_bn2 routine) 34 15 35 USE divhor ! horizontal divergence (div_hor routine)36 USE dynadv ! advection (dyn_adv routine)37 USE dynbfr ! Bottom friction terms (dyn_bfr routine)38 USE dynvor ! vorticity term (dyn_vor routine)39 USE dynhpg ! hydrostatic pressure grad. (dyn_hpg routine)40 USE dynldf ! lateral momentum diffusion (dyn_ldf routine)41 USE dynzdf ! vertical diffusion (dyn_zdf routine)42 USE dynspg ! surface pressure gradient (dyn_spg routine)43 44 USE dynnxt ! time-stepping (dyn_nxt routine)45 46 USE stopar ! Stochastic parametrization (sto_par routine)47 USE stopts48 49 USE bdy_par ! for lk_bdy50 USE bdy_oce ! for dmp logical51 USE bdydta ! open boundary condition data (bdy_dta routine)52 USE bdytra ! bdy cond. for tracers (bdy_tra routine)53 USE bdydyn3d ! bdy cond. for baroclinic vel. (bdy_dyn3d routine)54 55 USE sshwzv ! vertical velocity and ssh (ssh_nxt routine)56 ! (ssh_swp routine)57 ! (wzv routine)58 USE domvvl ! variable vertical scale factors (dom_vvl_sf_nxt routine)59 ! (dom_vvl_sf_swp routine)60 61 USE ldfslp ! iso-neutral slopes (ldf_slp routine)62 USE ldfdyn ! lateral eddy viscosity coef. (ldf_dyn routine)63 USE ldftra ! lateral eddy diffusive coef. (ldf_tra routine)64 65 USE zdftmx ! tide-induced vertical mixing (zdf_tmx routine)66 USE zdfbfr ! bottom friction (zdf_bfr routine)67 USE zdftke ! TKE vertical mixing (zdf_tke routine)68 USE zdfgls ! GLS vertical mixing (zdf_gls routine)69 USE zdfddm ! double diffusion mixing (zdf_ddm routine)70 USE zdfevd ! enhanced vertical diffusion (zdf_evd routine)71 USE zdfric ! Richardson vertical mixing (zdf_ric routine)72 USE zdfmxl ! Mixed-layer depth (zdf_mxl routine)73 74 USE step_diu ! Time stepping for diurnal sst75 USE diurnal_bulk ! diurnal SST bulk routines (diurnal_sst_takaya routine)76 USE cool_skin ! diurnal cool skin correction (diurnal_sst_coolskin routine)77 USE sbc_oce ! surface fluxes78 79 USE zpshde ! partial step: hor. derivative (zps_hde routine)80 81 USE diawri ! Standard run outputs (dia_wri routine)82 USE diaptr ! poleward transports (dia_ptr routine)83 USE diadct ! sections transports (dia_dct routine)84 USE diaar5 ! AR5 diagnosics (dia_ar5 routine)85 USE diahth ! thermocline depth (dia_hth routine)86 USE diafwb ! freshwater budget (dia_fwb routine)87 USE diahsb ! heat, salt and volume budgets (dia_hsb routine)88 USE diaharm89 USE diacfl90 USE flo_oce ! floats variables91 USE floats ! floats computation (flo_stp routine)92 93 USE crsfld ! Standard output on coarse grid (crs_fld routine)94 95 USE asminc ! assimilation increments (tra_asm_inc routine)96 ! (dyn_asm_inc routine)97 USE asmbkg98 USE stpctl ! time stepping control (stp_ctl routine)99 USE restart ! ocean restart (rst_wri routine)100 16 USE prtctl ! Print control (prt_ctl routine) 101 102 USE diaobs ! Observation operator103 17 104 18 USE in_out_manager ! I/O manager … … 110 24 USE xios 111 25 #endif 112 #if defined key_agrif113 USE agrif_opa_sponge ! Momemtum and tracers sponges114 USE agrif_opa_update ! Update (2-way nesting)115 #endif116 #if defined key_top117 USE trcstp ! passive tracer time-stepping (trc_stp routine)118 #endif119 26 !!---------------------------------------------------------------------- 120 27 !! NEMO/OPA 3.7 , NEMO Consortium (2014)
Note: See TracChangeset
for help on using the changeset viewer.