Changeset 6736 for branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM
- Timestamp:
- 2016-06-24T09:50:27+02:00 (8 years ago)
- Location:
- branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 12 edited
-
closea.F90 (modified) (10 diffs)
-
daymod.F90 (modified) (2 diffs)
-
dom_oce.F90 (modified) (8 diffs)
-
domain.F90 (modified) (5 diffs)
-
domhgr.F90 (modified) (1 diff)
-
domvvl.F90 (modified) (1 diff)
-
domwri.F90 (modified) (2 diffs)
-
domzgr.F90 (modified) (15 diffs)
-
domzgr_substitute.h90 (modified) (5 diffs)
-
dtatsd.F90 (modified) (6 diffs)
-
istate.F90 (modified) (10 diffs)
-
phycst.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r3632 r6736 19 19 USE oce ! dynamics and tracers 20 20 USE dom_oce ! ocean space and time domain 21 USE phycst ! physical constants22 21 USE in_out_manager ! I/O manager 23 22 USE sbc_oce ! ocean surface boundary conditions … … 185 184 !! put as run-off in open ocean. 186 185 !! 187 !! ** Action : emp updated surface freshwater fluxes and associated heat contentat kt186 !! ** Action : emp, emps updated surface freshwater fluxes at kt 188 187 !!---------------------------------------------------------------------- 189 188 INTEGER, INTENT(in) :: kt ! ocean model time step … … 192 191 REAL(wp), PARAMETER :: rsmall = 1.e-20_wp ! Closed sea correction epsilon 193 192 REAL(wp) :: zze2, ztmp, zcorr ! 194 REAL(wp) :: zcoef, zcoef1 !195 193 COMPLEX(wp) :: ctmp 196 194 REAL(wp), DIMENSION(jpncs) :: zfwf ! 1D workspace … … 245 243 ENDIF 246 244 ! !--------------------! 247 ! ! update emp 245 ! ! update emp, emps ! 248 246 zfwf = 0.e0_wp !--------------------! 249 247 IF( lk_mpp_rep ) THEN ! MPP reproductible calculation … … 284 282 ! 285 283 IF( ncstt(jc) == 0 ) THEN ! water/evap excess is shared by all open ocean 286 zcoef = zfwf(jc) / surf(jpncs+1) 287 zcoef1 = rcp * zcoef 288 emp(:,:) = emp(:,:) + zcoef 289 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 284 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 285 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 290 286 ! accumulate closed seas correction 291 zcorr = zcorr + zcoef287 zcorr = zcorr + zfwf(jc) / surf(jpncs+1) 292 288 ! 293 289 ELSEIF( ncstt(jc) == 1 ) THEN ! Excess water in open sea, at outflow location, excess evap shared … … 298 294 IF ( ji > 1 .AND. ji < jpi & 299 295 .AND. jj > 1 .AND. jj < jpj ) THEN 300 zcoef = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 301 zcoef1 = rcp * zcoef 302 emp(ji,jj) = emp(ji,jj) + zcoef 303 qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 296 emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 297 emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 304 298 ENDIF 305 299 END DO 306 300 ELSE 307 zcoef = zfwf(jc) / surf(jpncs+1) 308 zcoef1 = rcp * zcoef 309 emp(:,:) = emp(:,:) + zcoef 310 qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:) 301 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 302 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 311 303 ! accumulate closed seas correction 312 zcorr = zcorr + zcoef304 zcorr = zcorr + zfwf(jc) / surf(jpncs+1) 313 305 ENDIF 314 306 ELSEIF( ncstt(jc) == 2 ) THEN ! Excess e-p-r (either sign) goes to open ocean, at outflow location … … 318 310 IF( ji > 1 .AND. ji < jpi & 319 311 .AND. jj > 1 .AND. jj < jpj ) THEN 320 zcoef = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 321 zcoef1 = rcp * zcoef 322 emp(ji,jj) = emp(ji,jj) + zcoef 323 qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj) 312 emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 313 emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 324 314 ENDIF 325 315 END DO … … 328 318 DO jj = ncsj1(jc), ncsj2(jc) 329 319 DO ji = ncsi1(jc), ncsi2(jc) 330 zcoef = zfwf(jc) / surf(jc) 331 zcoef1 = rcp * zcoef 332 emp(ji,jj) = emp(ji,jj) - zcoef 333 qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj) 320 emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc) 321 emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc) 334 322 END DO 335 323 END DO … … 342 330 DO jj = ncsj1(jc), ncsj2(jc) 343 331 DO ji = ncsi1(jc), ncsi2(jc) 344 emp (ji,jj) = emp(ji,jj) - zcorr345 qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)332 emp (ji,jj) = emp (ji,jj) - zcorr 333 emps(ji,jj) = emps(ji,jj) - zcorr 346 334 END DO 347 335 END DO … … 350 338 ! 351 339 emp (:,:) = emp (:,:) * tmask(:,:,1) 340 emps(:,:) = emps(:,:) * tmask(:,:,1) 352 341 ! 353 342 CALL lbc_lnk( emp , 'T', 1._wp ) 343 CALL lbc_lnk( emps, 'T', 1._wp ) 354 344 ! 355 345 IF( nn_timing == 1 ) CALL timing_stop('sbc_clo') -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r3851 r6736 32 32 USE ioipsl, ONLY : ymds2ju ! for calendar 33 33 USE prtctl ! Print control 34 USE restart ! 34 35 USE trc_oce, ONLY : lk_offline ! offline flag 35 36 USE timing ! Timing 36 USE restart ! restart37 37 38 38 IMPLICIT NONE … … 153 153 IF ( nleapy == 1 ) THEN ! we are using calandar with leap years 154 154 IF ( MOD(nyear-1, 4) == 0 .AND. ( MOD(nyear-1, 400) == 0 .OR. MOD(nyear-1, 100) /= 0 ) ) THEN 155 nyear_len(0) 156 ENDIF 157 IF ( MOD(nyear , 4) == 0 .AND. ( MOD(nyear , 400) == 0 .OR. MOD(nyear, 100) /= 0 ) ) THEN155 nyear_len(0) = 366 156 ENDIF 157 IF ( MOD(nyear, 4) == 0 .AND. ( MOD(nyear, 400) == 0 .OR. MOD(nyear, 100) /= 0 ) ) THEN 158 158 nmonth_len(2) = 29 159 nyear_len(1) = 366 160 ENDIF 161 IF ( MOD(nyear+1, 4) == 0 .AND. ( MOD(nyear+1, 400) == 0 .OR. MOD(nyear+1, 100) /= 0 ) ) THEN 162 nyear_len(2) = 366 159 nyear_len(1) = 366 163 160 ENDIF 164 161 ENDIF -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r3851 r6736 8 8 !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 9 9 !! 4.0 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation 10 !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Add arrays associated11 !! to the optimization of BDY communications12 10 !!---------------------------------------------------------------------- 13 11 … … 82 80 INTEGER, PUBLIC :: narea !: number for local area 83 81 INTEGER, PUBLIC :: nbondi, nbondj !: mark of i- and j-direction local boundaries 84 INTEGER, ALLOCATABLE, PUBLIC :: nbondi_bdy(:) !: mark i-direction local boundaries for BDY open boundaries85 INTEGER, ALLOCATABLE, PUBLIC :: nbondj_bdy(:) !: mark j-direction local boundaries for BDY open boundaries86 INTEGER, ALLOCATABLE, PUBLIC :: nbondi_bdy_b(:) !: mark i-direction of neighbours local boundaries for BDY open boundaries87 INTEGER, ALLOCATABLE, PUBLIC :: nbondj_bdy_b(:) !: mark j-direction of neighbours local boundaries for BDY open boundaries88 89 82 INTEGER, PUBLIC :: npolj !: north fold mark (0, 3 or 4) 90 83 INTEGER, PUBLIC :: nlci, nldi, nlei !: i-dimensions of the local subdomain and its first and last indoor indices … … 131 124 LOGICAL, PUBLIC :: ln_zps = .FALSE. !: z-coordinate - partial step 132 125 LOGICAL, PUBLIC :: ln_sco = .FALSE. !: s-coordinate or hybrid z-s coordinate 133 126 LOGICAL, PUBLIC :: ln_s_sigma = .FALSE. ! use hybrid s-sigma -coordinate & stretching function 127 LOGICAL, PUBLIC :: ln_hyb = .FALSE. !: MANE1 s-coordinate or hybrid z-s coordinate 134 128 !! All coordinates 135 129 !! --------------- … … 172 166 !! =----------------======--------------- 173 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsigt, gsigw !: model level depth coefficient at t-, w-levels (analytic) 168 #if defined key_smsh 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsigt3, gsigw3 !: model level depth coefficient for sigma_s levels 170 #endif 174 171 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsi3w !: model level depth coefficient at w-level (sum of gsigw) 172 #if defined key_smsh 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsi3w3 !: model level depth coefficient for sigma_s levels 174 #endif 175 175 176 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: esigt, esigw !: vertical scale factor coef. at t-, w-levels 176 177 178 #if defined key_smsh 179 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: esigt3, esigw3 !: vertical scale factor coef. for sigma_S levels 180 #endif 177 181 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatv , hbatf !: ocean depth at the vertical of V--F 178 182 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatt , hbatu !: T--U points (m) … … 181 185 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at V--F 182 186 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing T--U points (m) 183 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio184 187 185 188 !!---------------------------------------------------------------------- … … 218 221 REAL(wp), PUBLIC :: adatrj !: number of elapsed days since the begining of the whole simulation 219 222 ! !: (cumulative duration of previous runs that may have used different time-step size) 220 INTEGER , PUBLIC, DIMENSION(0: 2) :: nyear_len !: length in days of the previous/current/next year223 INTEGER , PUBLIC, DIMENSION(0: 1) :: nyear_len !: length in days of the previous/current year 221 224 INTEGER , PUBLIC, DIMENSION(0:13) :: nmonth_len !: length in days of the months of the current year 222 225 INTEGER , PUBLIC, DIMENSION(0:13) :: nmonth_half !: second since Jan 1st 0h of the current year and the half of the months … … 293 296 & hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , STAT=ierr(6) ) 294 297 ! 295 ALLOCATE( gdept_0(jpk) , gdepw_0(jpk) , & 296 & e3t_0 (jpk) , e3w_0 (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 297 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 298 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 298 ALLOCATE( gdept_0(jpk) , gdepw_0(jpk) , & 299 & e3t_0 (jpk) , e3w_0 (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 300 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 301 #if defined key_smsh 302 & gsigt3 (jpi,jpj,jpk) , gsigw3 (jpi,jpj,jpk) , & 303 & esigt3 (jpi,jpj,jpk) , esigw3 (jpi,jpj,jpk) , & 304 & gsi3w3 (jpi,jpj,jpk) , & 305 #endif 306 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 299 307 ! 300 308 ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) , & … … 302 310 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 303 311 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 304 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) ,STAT=ierr(8) )312 & hift (jpi,jpj) , hifu (jpi,jpj) , STAT=ierr(8) ) 305 313 306 314 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r3764 r6736 36 36 USE dyncor_c1d ! Coriolis term (c1d case) (cor_c1d routine) 37 37 USE timing ! Timing 38 USE lbclnk ! ocean lateral boundary condition (or mpp link)39 38 40 39 IMPLICIT NONE … … 85 84 CALL dom_zgr ! Vertical mesh and bathymetry 86 85 CALL dom_msk ! Masks 87 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency88 86 IF( lk_vvl ) CALL dom_vvl ! Vertical variable mesh 89 87 ! … … 123 121 NAMELIST/namrun/ nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, & 124 122 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 125 & nn_write, ln_dimgnnn, ln_mskland , ln_clobber , nn_chunksz 123 & nn_write, ln_dimgnnn, ln_mskland , ln_clobber , nn_chunksz, ln_fse3t_b 126 124 NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh , rn_hmin, & 127 125 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin , & … … 156 154 WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber 157 155 WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz 156 WRITE(numout,*) ' fse3t_b in restart? ln_fse3t_b = ', ln_fse3t_b 158 157 ENDIF 159 158 … … 320 319 END SUBROUTINE dom_ctl 321 320 322 SUBROUTINE dom_stiff323 !!----------------------------------------------------------------------324 !! *** ROUTINE dom_stiff ***325 !!326 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency327 !!328 !! ** Method : Compute Haney (1991) hydrostatic condition ratio329 !! Save the maximum in the vertical direction330 !! (this number is only relevant in s-coordinates)331 !!332 !! Haney, R. L., 1991: On the pressure gradient force333 !! over steep topography in sigma coordinate ocean models.334 !! J. Phys. Oceanogr., 21, 610???619.335 !!----------------------------------------------------------------------336 INTEGER :: ji, jj, jk337 REAL(wp) :: zrxmax338 REAL(wp), DIMENSION(4) :: zr1339 !!----------------------------------------------------------------------340 rx1(:,:) = 0.e0341 zrxmax = 0.e0342 zr1(:) = 0.e0343 344 DO ji = 2, jpim1345 DO jj = 2, jpjm1346 DO jk = 1, jpkm1347 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji-1,jj ,jk ) &348 & +gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1)) &349 & /(gdepw(ji ,jj ,jk )+gdepw(ji-1,jj ,jk ) &350 & -gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1) + rsmall) )351 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw(ji+1,jj ,jk )-gdepw(ji ,jj ,jk ) &352 & +gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1)) &353 & /(gdepw(ji+1,jj ,jk )+gdepw(ji ,jj ,jk ) &354 & -gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )355 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw(ji ,jj+1,jk )-gdepw(ji ,jj ,jk ) &356 & +gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1)) &357 & /(gdepw(ji ,jj+1,jk )+gdepw(ji ,jj ,jk ) &358 & -gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )359 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji ,jj-1,jk ) &360 & +gdepw(ji ,jj ,jk+1)-gdepw(ji ,jj-1,jk+1)) &361 & /(gdepw(ji ,jj ,jk )+gdepw(ji ,jj-1,jk ) &362 & -gdepw(ji, jj ,jk+1)-gdepw(ji ,jj-1,jk+1) + rsmall) )363 zrxmax = MAXVAL(zr1(1:4))364 rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)365 END DO366 END DO367 END DO368 369 CALL lbc_lnk( rx1, 'T', 1. )370 371 zrxmax = MAXVAL(rx1)372 373 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain374 375 IF(lwp) THEN376 WRITE(numout,*)377 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax378 WRITE(numout,*) '~~~~~~~~~'379 ENDIF380 381 END SUBROUTINE dom_stiff382 383 384 385 321 !!====================================================================== 386 322 END MODULE domain -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r3294 r6736 26 26 USE lib_mpp ! MPP library 27 27 USE timing ! Timing 28 28 !! test - delete this line 29 29 IMPLICIT NONE 30 30 PRIVATE -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r3294 r6736 443 443 ! End of individual corrections to scale factors 444 444 445 #if ! defined key_melange 445 446 IF( ln_zps ) THEN ! minimum of the e3t at partial cell level 447 #endif 446 448 DO jj = 2, jpjm1 447 449 DO ji = fs_2, fs_jpim1 448 450 iku = mbku(ji,jj) 451 #if defined key_melange 452 IF(iku>39) THEN 453 #endif 454 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 455 #if defined key_melange 456 ENDIF 457 #endif 449 458 ikv = mbkv(ji,jj) 450 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 459 #if defined key_melange 460 IF(ikv>39) THEN 461 #endif 451 462 pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji ,jj+1,ikv) ) 452 END DO 453 END DO 454 ENDIF 463 #if defined key_melange 464 ENDIF 465 #endif 466 END DO 467 END DO 468 #if ! defined key_melange 469 ENDIF 470 #endif 455 471 456 472 pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:) ! anomaly to avoid zero along closed boundary/extra halos -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r3680 r6736 172 172 173 173 IF( ln_sco ) THEN ! s-coordinate 174 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) 175 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 174 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) ! ! depth 175 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 176 176 CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 177 177 CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 178 178 ! 179 CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef. 180 CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 181 CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w ) 182 CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) 183 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 179 #if defined key_smsh 180 CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt3 ) ! ! scaling coef. 181 CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw3 ) 182 CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w3 ) 183 CALL iom_rstput( 0, 0, inum4, 'esigt', esigt3 ) 184 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw3 ) 185 #else 186 CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef. 187 CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) 188 CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w ) 189 CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) 190 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 191 #endif 192 184 193 ! 185 194 CALL iom_rstput( 0, 0, inum4, 'e3t', e3t ) ! ! scale factors … … 187 196 CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 188 197 CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 189 CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 ) ! ! Max. grid stiffness ratio 190 ! 191 CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept ) ! ! stretched system 192 CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw ) 198 ! 199 #if defined key_smsh 200 CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 ) 201 DO jk = 1,jpk 202 DO jj = 1, jpjm1 203 DO ji = 1, fs_jpim1 ! vector opt. 204 zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj ,jk) ) 205 zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji ,jj+1,jk) ) 206 END DO 207 END DO 208 END DO 209 CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) 210 CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 ) 211 CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 ) 212 CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 ) 213 DO jj = 1,jpj 214 DO ji = 1,jpi 215 zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj) ) * tmask(ji,jj,1) 216 zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1) 217 END DO 218 END DO 219 CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 ) 220 CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 221 #endif 222 CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 ) ! ! stretched system 223 CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 193 224 ENDIF 194 225 -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3764 r6736 15 15 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 16 16 !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function18 17 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 18 !!---------------------------------------------------------------------- 20 19 21 20 !!---------------------------------------------------------------------- … … 29 28 !! zgr_zps : z-coordinate with partial steps 30 29 !! zgr_sco : s-coordinate 31 !! fssig : tanh stretch function 32 !! fssig1 : Song and Haidvogel 1994 stretch function 33 !! fgamma : Siddorn and Furner 2012 stretching function 30 !! fssig : sigma coordinate non-dimensional function 31 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) 34 32 !!--------------------------------------------------------------------- 35 33 USE oce ! ocean variables … … 50 48 51 49 ! !!* Namelist namzgr_sco * 52 LOGICAL :: ln_s_sh94 = .false. ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)53 LOGICAL :: ln_s_sf12 = .true. ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)54 !55 50 REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) 56 51 REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 57 REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1)58 REAL(wp) :: rn_hc = 150._wp ! Critical depth for transition from sigma to stretched coordinates59 ! Song and Haidvogel 1994 stretching parameters60 52 REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20) 61 53 REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1) 62 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter 54 REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1) 55 ! LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 56 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter for song and haidvogel stretching 63 57 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 64 ! Siddorn and Furner stretching parameters 65 LOGICAL :: ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 66 REAL(wp) :: rn_alpha = 4.4_wp ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 67 REAL(wp) :: rn_efold = 0.0_wp ! efold length scale for transition to stretched coord 68 REAL(wp) :: rn_zs = 1.0_wp ! depth of surface grid box 69 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 70 REAL(wp) :: rn_zb_a = 0.024_wp ! bathymetry scaling factor for calculating Zb 71 REAL(wp) :: rn_zb_b = -0.2_wp ! offset for calculating Zb 58 REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates 59 REAL(wp) :: rn_zsigma = 300._wp ! Maximum depth for s-sigma layer 60 INTEGER :: nn_sig_lev = 10 ! Maximum number of levels of s-sigma layer 61 REAL(wp) :: rn_kth = 15._wp ! Approximate layer number, beyond which streching will be maximum 62 REAL(wp) :: rn_acr = 9.00_wp ! 72 63 73 64 !! * Substitutions … … 100 91 INTEGER :: ioptio, ibat ! local integer 101 92 ! 102 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 93 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco , ln_hyb 103 94 !!---------------------------------------------------------------------- 104 95 ! … … 116 107 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 117 108 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 109 WRITE(numout,*) ' hybrid s-z-coordinates,s at shelf ln_hyb = ', ln_hyb 110 118 111 ENDIF 119 112 … … 131 124 IF( ln_zco ) CALL zgr_zco ! z-coordinate 132 125 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 133 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 126 IF( ln_sco.AND. .NOT. ln_hyb ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate (z at upper levels ) 127 IF( ln_sco .AND. ln_hyb ) CALL zgr_hyb ! hybrid s-sigma z ( s- at shel 134 128 ! 135 129 ! final adjustment of mbathy & check … … 520 514 ENDIF 521 515 ! 516 ! 522 517 CALL wrk_dealloc( jpidta, jpjdta, idta ) 523 518 CALL wrk_dealloc( jpidta, jpjdta, zdta ) … … 639 634 END DO 640 635 END DO 636 IF( lk_mpp ) CALL mpp_sum( icompt ) 641 637 IF( icompt == 0 ) THEN 642 638 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' … … 1053 1049 END SUBROUTINE zgr_zps 1054 1050 1051 1052 FUNCTION fssig( pk ) RESULT( pf ) 1053 !!---------------------------------------------------------------------- 1054 !! *** ROUTINE eos_init *** 1055 !! 1056 !! ** Purpose : provide the analytical function in s-coordinate 1057 !! 1058 !! ** Method : the function provide the non-dimensional position of 1059 !! T and W (i.e. between 0 and 1) 1060 !! T-points at integer values (between 1 and jpk) 1061 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1062 !!---------------------------------------------------------------------- 1063 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 1064 REAL(wp) :: pf ! sigma value 1065 !!---------------------------------------------------------------------- 1066 ! 1067 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 1068 & - TANH( rn_thetb * rn_theta ) ) & 1069 & * ( COSH( rn_theta ) & 1070 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1071 & / ( 2._wp * SINH( rn_theta ) ) 1072 ! 1073 END FUNCTION fssig 1074 1075 1076 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 1077 !!---------------------------------------------------------------------- 1078 !! *** ROUTINE eos_init *** 1079 !! 1080 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 1081 !! 1082 !! ** Method : the function provides the non-dimensional position of 1083 !! T and W (i.e. between 0 and 1) 1084 !! T-points at integer values (between 1 and jpk) 1085 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1086 !!---------------------------------------------------------------------- 1087 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 1088 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 1089 REAL(wp) :: pf1 ! sigma value 1090 !!---------------------------------------------------------------------- 1091 ! 1092 IF ( rn_theta == 0._wp ) then ! uniform sigma 1093 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 1094 ELSE ! stretched sigma 1095 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 1096 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1097 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1098 ENDIF 1099 ! 1100 END FUNCTION fssig1 1101 FUNCTION fssig2 ( pk, kmax ) RESULT( pf2 ) 1102 !!---------------------------------------------------------------------- 1103 !! *** ROUTINE eos_init *** 1104 !! 1105 !! ** Purpose : provide the analytical function in s-coordinate 1106 !! 1107 !! ** Method : the function provide the non-dimensional position of 1108 !! T and W (i.e. between 0 and 1) 1109 !! T-points at integer values (between 1 and kmax ) 1110 !! W-points at integer values - 1/2 (between 0.5 and kmax-0.5) 1111 !! 1112 !! Reference : ??? 1113 !!---------------------------------------------------------------------- 1114 REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate 1115 REAL(wp) :: pf2 ! sigma value 1116 INTEGER, INTENT (in) :: kmax ! max of sigma)level 1117 !!---------------------------------------------------------------------- 1118 ! 1119 pf2 = ( TANH( rn_theta * ( -(pk-0.5) / REAL(kmax-1,wp) + rn_thetb ) ) & 1120 & - TANH( rn_thetb * rn_theta ) ) & 1121 & * ( COSH( rn_theta ) & 1122 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1123 & / ( 2._wp * SINH( rn_theta ) ) 1124 ! 1125 END FUNCTION fssig2 1126 1127 FUNCTION fssig3( pk1, pbb ,kmax ) RESULT( pf3 ) 1128 !!---------------------------------------------------------------------- 1129 !! *** ROUTINE eos_init *** 1130 !! 1131 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 1132 !! 1133 !! ** Method : the function provides the non-dimensional position of 1134 !! T and W (i.e. between 0 and 1) 1135 !! T-points at integer values (between 1 and jpk) 1136 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1137 !! 1138 !! Reference : ??? 1139 !!---------------------------------------------------------------------- 1140 REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate 1141 REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient 1142 REAL(wp) :: pf3 ! sigma value 1143 INTEGER, INTENT (in) :: kmax ! max number of s-sigma levels 1144 !!---------------------------------------------------------------------- 1145 ! 1146 IF ( rn_theta == 0 ) then ! uniform sigma 1147 pf3 = -(pk1-0.5_wp) / REAL( kmax-1,wp ) 1148 ELSE ! stretched sigma 1149 pf3 = (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5_wp)/REAL(kmax-1,wp)) ) ) / sinh(rn_theta) + & 1150 & pbb * ( (tanh( rn_theta*( (-(pk1-0.5_wp)/REAL(kmax-1,wp)) + 0.5_wp) ) - tanh(0.5*rn_theta) ) / & 1151 & (2._wp*tanh(0.5_wp*rn_theta) ) ) 1152 ENDIF 1153 END FUNCTION fssig3 1154 1155 SUBROUTINE fszref (zkth, zdzmin, zacr, zhmax,jpup,zhsigm ) 1156 INTEGER :: jk ! dummy loop indices 1157 REAL(wp) :: zt, zw ! temporary scalars 1158 REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in 1159 REAL(wp) :: zacr, zdzmin, zhmax, zhmax_r ! read from namelist or par_XXX.h90 1160 REAL(wp) :: zhsigm ! depth of sigma layer 1161 INTEGER :: jpup, jpkmax ! the last sigma level and number of z-levels 1162 !!---------------------------------------------------------------------- 1163 ! compute reference depth leveles 1164 ! Set variables from parameters 1165 ! ------------------------------ 1166 ! zkth = rn_kth ; zacr = rn_acr 1167 ! zdzmin = rn_dzmin ; zhmax_r = rn_hmax 1168 1169 ! za0, za1, zsur are computed from zdzmin , zhmax, zkth, zacr 1170 ! 1171 jpkmax= jpk - jpup 1172 zhmax_r = zhmax - zhsigm 1173 1174 za1 = ( zdzmin - zhmax_r / REAL(jpkmax,wp) ) & 1175 & / ( TANH((1-zkth)/zacr) - zacr/REAL(jpkmax,wp) * ( LOG( COSH( (jpkmax + 1 - zkth) / zacr) ) & 1176 & - LOG( COSH( ( 1 - zkth) / zacr) ) ) ) 1177 za0 = zdzmin - za1 * TANH( (1-zkth) / zacr ) 1178 zsur = - za0 - za1 * zacr * LOG( COSH( (1-zkth) / zacr ) ) 1179 ! Reference z-coordinate (depth - scale factor at T- and W-points) 1180 ! ====================== 1181 IF( zkth == 0.e0 ) THEN ! uniform vertical grid 1182 za1 = zhmax_r / REAL(jpkmax-1,wp) 1183 DO jk = 1, jpkmax+1 1184 zw = REAL( jk,wp ) 1185 zt = REAL( jk,wp ) + 0.5_wp 1186 gdepw_0(jk+jpup-1 ) = ( zw - 1 ) * za1 + zhsigm 1187 gdept_0(jk+jpup-1 ) = ( zt - 1 ) * za1 + zhsigm 1188 e3w_0 (jk+jpup-1 ) = za1 1189 e3t_0 (jk+jpup-1 ) = za1 1190 1191 END DO 1192 ELSE ! Madec & Imbard 1996 function 1193 DO jk = 1, jpkmax+1 1194 zw = REAL( jk,wp) 1195 zt = REAL( jk,wp ) + 0.5_wp 1196 gdepw_0(jk+jpup-1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )+zhsigm 1197 gdept_0(jk+jpup-1) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )+zhsigm 1198 e3w_0 (jk+jpup-1) = za0 + za1 * TANH( (zw-zkth) / zacr ) 1199 e3t_0 (jk+jpup-1) = za0 + za1 * TANH( (zt-zkth) / zacr ) 1200 1201 END DO 1202 gdepw_0(jpup) = zhsigm ! force first w-level to be exactly at zhsigm 1203 ENDIF 1204 IF(lwp) WRITE (numout,*) " max and min z-vertical level",jpkmax+1,jpup 1205 1206 END SUBROUTINE fszref 1207 1208 1055 1209 SUBROUTINE zgr_sco 1056 1210 !!---------------------------------------------------------------------- … … 1071 1225 !! hbatv = mj( hbatt ) 1072 1226 !! hbatf = mi( mj( hbatt ) ) 1073 !! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical1227 !! - Compute gsigt, gsigw, esigt, esigw from an analytical 1074 1228 !! function and its derivative given as function. 1075 !! z_gsigt(k) = fssig (k )1076 !! z_gsigw(k) = fssig (k-0.5)1077 !! z_esigt(k) = fsdsig(k )1078 !! z_esigw(k) = fsdsig(k-0.5)1079 !! Th ree options for stretching are give, and they canbe modified1080 !! following the user s requirements. Nevertheless, the output as1229 !! gsigt(k) = fssig (k ) 1230 !! gsigw(k) = fssig (k-0.5) 1231 !! esigt(k) = fsdsig(k ) 1232 !! esigw(k) = fsdsig(k-0.5) 1233 !! This routine is given as an example, it must be modified 1234 !! following the user s desiderata. nevertheless, the output as 1081 1235 !! well as the way to compute the model levels and scale factors 1082 !! must be respected in order to insure second order a ccuracy1236 !! must be respected in order to insure second order a!!uracy 1083 1237 !! schemes. 1084 1238 !! 1085 !! The three methods for stretching available are: 1086 !! 1087 !! s_sh94 (Song and Haidvogel 1994) 1088 !! a sinh/tanh function that allows sigma and stretched sigma 1089 !! 1090 !! s_sf12 (Siddorn and Furner 2012?) 1091 !! allows the maintenance of fixed surface and or 1092 !! bottom cell resolutions (cf. geopotential coordinates) 1093 !! within an analytically derived stretched S-coordinate framework. 1094 !! 1095 !! s_tanh (Madec et al 1996) 1096 !! a cosh/tanh function that gives stretched coordinates 1097 !! 1239 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1098 1240 !!---------------------------------------------------------------------- 1099 1241 ! 1100 1242 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1101 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1102 REAL(wp) :: zrmax, ztaper ! temporary scalars 1103 ! 1104 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1105 1106 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1107 rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1108 !!---------------------------------------------------------------------- 1243 INTEGER :: inum ! temporary logical unit 1244 INTEGER :: iip1, ijp1, iim1, ijm1, kdep ! temporary integers 1245 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper, maxzenv ! temporary scalars 1246 #if defined key_melange 1247 REAL(wp) :: rn_hc_bak ! temporary scalars 1248 #endif 1249 REAL(wp) :: zrfact ! temporary scalars 1250 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1251 ! 1252 #if defined key_fudge 1253 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, zri, zrj, zhbat, fenv 1254 #else 1255 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, zri, zrj, zhbat 1256 #endif 1257 #if defined key_smsh 1258 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1259 #else 1260 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1261 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1262 #endif 1263 NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 1264 #if defined key_melange 1265 NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc, nn_sig_lev 1266 #endif 1267 !!---------------------------------------------------------------------- 1109 1268 ! 1110 1269 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1111 1270 ! 1112 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1113 ! 1271 CALL wrk_alloc( jpi, jpj, ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1272 #if defined key_fudge 1273 CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat, fenv ) 1274 #else 1275 CALL wrk_alloc( jpi, jpj, zenv, zri, zrj, zhbat ) 1276 #endif 1277 #if defined key_smsh 1278 CALL wrk_alloc( jpi, jpj, jpk, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1279 #else 1280 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1281 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1282 #endif 1114 1283 REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters 1115 1284 READ ( numnam, namzgr_sco ) … … 1120 1289 WRITE(numout,*) '~~~~~~~~~~~' 1121 1290 WRITE(numout,*) ' Namelist namzgr_sco' 1122 WRITE(numout,*) ' stretching coeffs ' 1123 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max 1124 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min 1125 WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc 1126 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax 1127 WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 1128 WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients' 1129 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta 1130 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb 1131 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 1132 WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 1133 WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit 1134 WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients' 1135 WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha 1136 WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold 1137 WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs 1138 WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a 1139 WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b 1140 WRITE(numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 1141 ENDIF 1291 WRITE(numout,*) ' sigma-stretching coeffs ' 1292 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max 1293 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min 1294 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta 1295 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb 1296 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax 1297 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1298 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb 1299 WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc 1300 ENDIF 1301 1302 #if defined key_melange 1303 CALL zgr_zps ! Partial step z-coordinate 1304 ! Scale factors and depth at U-, V-, UW and VW-points 1305 DO jk = 1, nn_sig_lev ! initialisation to z-scale factors above ln_s_sigma to remove any zps 1306 e3u (:,:,jk) = e3t_0(jk) 1307 e3v (:,:,jk) = e3t_0(jk) 1308 e3uw(:,:,jk) = e3w_0(jk) 1309 e3vw(:,:,jk) = e3w_0(jk) 1310 END DO 1311 #endif 1312 1313 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1314 esigt3 = 0._wp ; esigw3 = 0._wp 1315 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1316 esigwu3 = 0._wp ; esigwv3 = 0._wp 1142 1317 1143 1318 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate … … 1158 1333 ! ! ============================= 1159 1334 ! use r-value to create hybrid coordinates 1335 zenv(:,:) = bathy(:,:) 1336 #if defined key_melange 1337 DO jj = 1, jpj 1338 DO ji = 1, jpi 1339 zenv(ji,jj) = MIN( bathy(ji,jj), gdepw_0(nn_sig_lev + 1) ) 1340 ENDDO 1341 ENDDO 1342 #endif 1343 #if defined key_fudge 1344 CALL iom_open ( 'fudge.nc', inum ) 1345 CALL iom_get ( inum, jpdom_data, 'zenv', fenv ) 1346 CALL iom_close( inum ) 1347 DO jj = 1, jpj 1348 DO ji = 1, jpi 1349 zenv(ji,jj) = MAX( zenv(ji,jj), fenv(ji,jj) ) 1350 ENDDO 1351 ENDDO 1352 #endif 1353 ! 1354 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1355 DO jj = 1, jpj 1356 DO ji = 1, jpi 1357 IF( bathy(ji,jj) == 0._wp ) THEN 1358 iip1 = MIN( ji+1, jpi ) 1359 ijp1 = MIN( jj+1, jpj ) 1360 iim1 = MAX( ji-1, 1 ) 1361 ijm1 = MAX( jj-1, 1 ) 1362 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1363 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1364 zenv(ji,jj) = rn_sbot_min 1365 ENDIF 1366 ENDIF 1367 END DO 1368 END DO 1369 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1370 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 1371 ! 1372 ! smooth the bathymetry (if required) 1373 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1374 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1375 ! 1376 jl = 0 1377 zrmax = 1._wp 1378 ! 1379 ! set scaling factor used in reducing vertical gradients 1380 zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 1381 ! 1382 ! initialise temporary evelope depth arrays 1383 ztmpi1(:,:) = zenv(:,:) 1384 ztmpi2(:,:) = zenv(:,:) 1385 ztmpj1(:,:) = zenv(:,:) 1386 ztmpj2(:,:) = zenv(:,:) 1387 ! 1388 ! initialise temporary r-value arrays 1389 zri(:,:) = 1._wp 1390 zrj(:,:) = 1._wp 1391 ! ! ================ ! 1392 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! 1393 ! ! ================ ! 1394 jl = jl + 1 1395 zrmax = 0._wp 1396 ! we set zrmax from previous r-values (zri and zrj) first 1397 ! if set after current r-value calculation (as previously) 1398 ! we could exit DO WHILE prematurely before checking r-value 1399 ! of current zenv 1400 DO jj = 1, nlcj 1401 DO ji = 1, nlci 1402 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 1403 END DO 1404 END DO 1405 zri(:,:) = 0._wp 1406 zrj(:,:) = 0._wp 1407 DO jj = 1, nlcj 1408 DO ji = 1, nlci 1409 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1410 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1411 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 1412 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1413 END IF 1414 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 1415 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1416 END IF 1417 IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 1418 IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 1419 IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 1420 IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 1421 END DO 1422 END DO 1423 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1424 ! 1425 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 1426 ! 1427 DO jj = 1, nlcj 1428 DO ji = 1, nlci 1429 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 1430 END DO 1431 END DO 1432 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1433 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 1434 ! ! ================ ! 1435 END DO ! End loop ! 1436 ! ! ================ ! 1437 DO jj = 1, jpj 1438 DO ji = 1, jpi 1439 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1440 END DO 1441 END DO 1442 ! 1443 ! Envelope bathymetry saved in hbatt 1444 hbatt(:,:) = zenv(:,:) 1445 1446 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1447 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 1448 DO jj = 1, jpj 1449 DO ji = 1, jpi 1450 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) 1451 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 1452 END DO 1453 END DO 1454 ENDIF 1455 ! 1456 IF(lwp) THEN ! Control print 1457 WRITE(numout,*) 1458 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 1459 WRITE(numout,*) 1460 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout ) 1461 IF( nprint == 1 ) THEN 1462 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 1463 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 1464 ENDIF 1465 ENDIF 1466 1467 ! ! ============================== 1468 ! ! hbatu, hbatv, hbatf fields 1469 ! ! ============================== 1470 IF(lwp) THEN 1471 WRITE(numout,*) 1472 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 1473 ENDIF 1474 hbatu(:,:) = rn_sbot_min 1475 hbatv(:,:) = rn_sbot_min 1476 hbatf(:,:) = rn_sbot_min 1477 DO jj = 1, jpjm1 1478 DO ji = 1, jpim1 ! NO vector opt. 1479 hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 1480 hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 1481 hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 1482 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 1483 END DO 1484 END DO 1485 ! 1486 ! Apply lateral boundary condition 1487 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 1488 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) 1489 DO jj = 1, jpj 1490 DO ji = 1, jpi 1491 IF( hbatu(ji,jj) == 0._wp ) THEN 1492 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 1493 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) 1494 ENDIF 1495 END DO 1496 END DO 1497 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) 1498 DO jj = 1, jpj 1499 DO ji = 1, jpi 1500 IF( hbatv(ji,jj) == 0._wp ) THEN 1501 IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min 1502 IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) 1503 ENDIF 1504 END DO 1505 END DO 1506 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) 1507 DO jj = 1, jpj 1508 DO ji = 1, jpi 1509 IF( hbatf(ji,jj) == 0._wp ) THEN 1510 IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min 1511 IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) 1512 ENDIF 1513 END DO 1514 END DO 1515 1516 !!bug: key_helsinki a verifer 1517 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 1518 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 1519 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 1520 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 1521 1522 IF( nprint == 1 .AND. lwp ) THEN 1523 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 1524 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 1525 WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & 1526 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 1527 WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 1528 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 1529 WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 1530 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 1531 ENDIF 1532 !! helsinki 1533 1534 ! ! ======================= 1535 ! ! s-ccordinate fields (gdep., e3.) 1536 ! ! ======================= 1537 ! 1538 ! non-dimensional "sigma" for model level depth at w- and t-levels 1539 1540 IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths 1541 ! ! below rn_hc, with uniform sigma in shallower waters 1542 DO ji = 1, jpi 1543 DO jj = 1, jpj 1544 1545 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1546 DO jk = 1, jpk 1547 #if defined key_melange 1548 gsigw3(ji,jj,jk) = gdepw_0(jk)/gdepw_0(nn_sig_lev + 1) 1549 gsigt3(ji,jj,jk) = gdept_0(jk)/gdepw_0(nn_sig_lev + 1) 1550 #else 1551 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1552 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1553 #endif 1554 END DO 1555 ELSE ! shallow water, uniform sigma 1556 DO jk = 1, jpk 1557 #if defined key_melange 1558 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(nn_sig_lev,wp) 1559 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(nn_sig_lev,wp) 1560 #else 1561 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1562 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1563 #endif 1564 END DO 1565 ENDIF 1566 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1567 ! 1568 DO jk = 1, jpkm1 1569 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1570 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1571 END DO 1572 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1573 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1574 ! 1575 ! Coefficients for vertical depth as the sum of e3w scale factors 1576 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1577 DO jk = 2, jpk 1578 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1579 END DO 1580 ! 1581 #if defined key_melange 1582 DO jk = 1, nn_sig_lev+1 1583 ! DO jk = 1, jpk 1584 IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt? 1585 #else 1586 DO jk = 1, jpk 1587 #endif 1588 #if defined key_melange 1589 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(nn_sig_lev,wp) 1590 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(nn_sig_lev,wp) 1591 ! zcoeft = ( REAL(MIN(jk,nn_sig_lev),wp) - 0.5_wp ) / REAL(nn_sig_lev-1,wp) 1592 ! zcoefw = ( REAL(MIN(jk,nn_sig_lev),wp) - 1.0_wp ) / REAL(nn_sig_lev-1,wp) 1593 #else 1594 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1595 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1596 #endif 1597 #if defined key_melange 1598 rn_hc_bak = rn_hc 1599 rn_hc = MIN( MAX ( & 1600 & (hbatt(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc)) & 1601 & ,0._wp) ,rn_hc) 1602 #endif 1603 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1604 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1605 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1606 #if defined key_melange 1607 rn_hc = rn_hc_bak 1608 #endif 1609 IF( gdepw(ji,jj,jk) < 0._wp ) THEN 1610 WRITE(*,*) 'zgr_sco : gdepw at point (i,j,k)= ', ji, jj, jk, (gsigw3(ji,jj,jk)*10000._wp-zcoefw*10000._wp) 1611 ENDIF 1612 #if defined key_melange 1613 ENDIF 1614 #endif 1615 END DO 1616 ! 1617 END DO ! for all jj's 1618 END DO ! for all ji's 1619 1620 DO ji = 1, jpim1 1621 DO jj = 1, jpjm1 1622 #if defined key_melange 1623 IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be bathy or hbatt? 1624 DO jk = 1, nn_sig_lev+1 ! scale factors should be the same in both zps and sco when H > Hcrit?? 1625 ! DO jk = 1, jpk 1626 #else 1627 DO jk = 1, jpk 1628 #endif 1629 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1630 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1631 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1632 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1633 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1634 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1635 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1636 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1637 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1638 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1639 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1640 ! 1641 #if defined key_melange 1642 rn_hc_bak = rn_hc 1643 rn_hc = MIN( MAX( & 1644 & (hbatt(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc)) & 1645 & ,0._wp) ,rn_hc) 1646 ! e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1647 ! e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1648 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1649 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1650 #else 1651 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1652 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1653 #endif 1654 #if defined key_melange 1655 rn_hc = MIN( MAX( & 1656 & (hbatu(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) & 1657 & ,0._wp) ,rn_hc_bak) 1658 ! e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1659 ! e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1660 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1661 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1662 #else 1663 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1664 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1665 #endif 1666 #if defined key_melange 1667 rn_hc = MIN( MAX( & 1668 & (hbatv(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) & 1669 & ,0._wp) ,rn_hc_bak) 1670 ! e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1671 ! e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1672 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1673 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1674 #else 1675 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1676 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1677 #endif 1678 #if defined key_melange 1679 rn_hc = MIN( MAX( & 1680 & (hbatf(ji,jj)-gdepw_0(nn_sig_lev + 1)) / (1._wp - (gdepw_0(nn_sig_lev + 1)/rn_hc_bak)) & 1681 & ,0._wp), rn_hc_bak) 1682 ! e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev - 1,wp) ) 1683 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(nn_sig_lev ,wp) ) 1684 #else 1685 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp)) 1686 #endif 1687 ! 1688 #if defined key_melange 1689 rn_hc = rn_hc_bak 1690 #endif 1691 END DO 1692 #if defined key_melange 1693 ENDIF 1694 #endif 1695 END DO 1696 END DO 1697 1698 CALL lbc_lnk( e3t , 'T', 1._wp ) 1699 CALL lbc_lnk( e3u , 'U', 1._wp ) 1700 CALL lbc_lnk( e3v , 'V', 1._wp ) 1701 CALL lbc_lnk( e3f , 'F', 1._wp ) 1702 CALL lbc_lnk( e3w , 'W', 1._wp ) 1703 CALL lbc_lnk( e3uw, 'U', 1._wp ) 1704 CALL lbc_lnk( e3vw, 'V', 1._wp ) 1705 1706 ! 1707 ELSE ! not ln_s_sigma 1708 ! 1709 DO jk = 1, jpk 1710 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1711 gsigt(jk) = -fssig( REAL(jk,wp) ) 1712 END DO 1713 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1714 ! 1715 ! Coefficients for vertical scale factors at w-, t- levels 1716 !!gm bug : define it from analytical function, not like juste bellow.... 1717 !!gm or betteroffer the 2 possibilities.... 1718 DO jk = 1, jpkm1 1719 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1720 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1721 END DO 1722 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1723 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1724 1725 !!gm original form 1726 !!org DO jk = 1, jpk 1727 !!org esigt(jk)=fsdsig( FLOAT(jk) ) 1728 !!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 1729 !!org END DO 1730 !!gm 1731 ! 1732 ! Coefficients for vertical depth as the sum of e3w scale factors 1733 gsi3w(1) = 0.5_wp * esigw(1) 1734 DO jk = 2, jpk 1735 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1736 END DO 1737 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1738 DO jk = 1, jpk 1739 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1740 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1741 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 1742 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 1743 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 1744 END DO 1745 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1746 DO jj = 1, jpj 1747 DO ji = 1, jpi 1748 DO jk = 1, jpk 1749 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1750 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1751 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1752 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1753 ! 1754 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1755 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1756 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1757 END DO 1758 END DO 1759 END DO 1760 ! 1761 ENDIF ! ln_s_sigma 1762 1763 where (e3t (:,:,:).eq.0._wp) e3t(:,:,:) = 1._wp 1764 where (e3u (:,:,:).eq.0._wp) e3u(:,:,:) = 1._wp 1765 where (e3v (:,:,:).eq.0._wp) e3v(:,:,:) = 1._wp 1766 where (e3f (:,:,:).eq.0._wp) e3f(:,:,:) = 1._wp 1767 where (e3w (:,:,:).eq.0._wp) e3w(:,:,:) = 1._wp 1768 where (e3uw (:,:,:).eq.0._wp) e3uw(:,:,:) = 1._wp 1769 where (e3vw (:,:,:).eq.0._wp) e3vw(:,:,:) = 1._wp 1770 1771 1772 fsdept(:,:,:) = gdept (:,:,:) 1773 fsdepw(:,:,:) = gdepw (:,:,:) 1774 fsde3w(:,:,:) = gdep3w(:,:,:) 1775 fse3t (:,:,:) = e3t (:,:,:) 1776 fse3u (:,:,:) = e3u (:,:,:) 1777 fse3v (:,:,:) = e3v (:,:,:) 1778 fse3f (:,:,:) = e3f (:,:,:) 1779 fse3w (:,:,:) = e3w (:,:,:) 1780 fse3uw(:,:,:) = e3uw (:,:,:) 1781 fse3vw(:,:,:) = e3vw (:,:,:) 1782 !! 1783 ! HYBRID : 1784 DO jj = 1, jpj 1785 DO ji = 1, jpi 1786 #if defined key_melange 1787 IF( bathy(ji,jj) < gdepw_0(nn_sig_lev + 1) ) THEN ! should this be hbatt or bathy 1788 DO jk = 1, nn_sig_lev 1789 ! DO jk = 1, jpkm1 1790 #else 1791 DO jk = 1, jpkm1 1792 #endif 1793 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1794 END DO 1795 #if defined key_melange 1796 ENDIF 1797 #endif 1798 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1799 END DO 1800 END DO 1801 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1802 & ' MAX ', MAXVAL( mbathy(:,:) ) 1803 1804 ! ! ============= 1805 IF(lwp) THEN ! Control print 1806 ! ! ============= 1807 WRITE(numout,*) 1808 WRITE(numout,*) ' domzgr: vertical coefficients for model level' 1809 WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" ) 1810 WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk ) 1811 ENDIF 1812 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 1813 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1814 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & 1815 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) ) 1816 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), & 1817 & ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), & 1818 & ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), & 1819 & ' w ', MINVAL( fse3w (:,:,:) ) 1820 1821 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & 1822 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) ) 1823 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), & 1824 & ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), & 1825 & ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), & 1826 & ' w ', MAXVAL( fse3w (:,:,:) ) 1827 ENDIF 1828 ! 1829 IF(lwp) THEN ! selected vertical profiles 1830 WRITE(numout,*) 1831 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 1832 WRITE(numout,*) ' ~~~~~~ --------------------' 1833 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1834 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), & 1835 & fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk ) 1836 DO jj = mj0(20), mj1(20) 1837 DO ji = mi0(20), mi1(20) 1838 WRITE(numout,*) 1839 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1840 WRITE(numout,*) ' ~~~~~~ --------------------' 1841 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1842 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1843 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1844 END DO 1845 END DO 1846 DO jj = mj0(74), mj1(74) 1847 DO ji = mi0(100), mi1(100) 1848 WRITE(numout,*) 1849 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1850 WRITE(numout,*) ' ~~~~~~ --------------------' 1851 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1852 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1853 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1854 END DO 1855 END DO 1856 ENDIF 1857 1858 !!gm bug? no more necessary? if ! defined key_helsinki 1859 DO jk = 1, jpk 1860 DO jj = 1, jpj 1861 DO ji = 1, jpi 1862 IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1863 WRITE(*,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk, fse3w(ji,jj,jk), fse3t(ji,jj,jk) 1864 ! WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1865 ! CALL ctl_stop( ctmp1 ) 1866 ENDIF 1867 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1868 WRITE(*,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk, fsdepw(ji,jj,jk), fsdept(ji,jj,jk) 1869 ! WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1870 ! CALL ctl_stop( ctmp1 ) 1871 ENDIF 1872 END DO 1873 END DO 1874 END DO 1875 !!gm bug #endif 1876 ! 1877 1878 CALL wrk_dealloc( jpi, jpj, zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat ) 1879 1880 #if defined key_smsh 1881 CALL wrk_dealloc( jpi, jpj, jpk, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1882 #else 1883 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1884 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1885 #endif 1886 ! 1887 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1888 ! 1889 END SUBROUTINE zgr_sco 1890 SUBROUTINE zgr_hyb 1891 !!---------------------------------------------------------------------- 1892 !! *** ROUTINE zgr_sco *** 1893 !! Combination of zgr_sco in upper layers ( shelf ) and zgr_zps in abyss !! 1894 !! ** Purpose : define the s-z coordinate system 1895 !! 1896 !! ** Method : s-coordinate in upper layers and z-coordinates below 1897 !! The depth of model levels is defined as the product of an 1898 !! analytical function by the local bathymetry, while the vertical 1899 !! scale factors are defined as the product of the first derivative 1900 !! of the analytical function by the bathymetry. 1901 !! (this solution save memory as depth and scale factors are not 1902 !! 3d fields) 1903 !! - Read bathymetry (in meters) at t-point and compute the 1904 !! bathymetry at u-, v-, and f-points. 1905 !! hbatu = mi( hbatt ) 1906 !! hbatv = mj( hbatt ) 1907 !! hbatf = mi( mj( hbatt ) ) 1908 !! - Compute gsigt, gsigw, esigt, esigw from an analytical 1909 !! function 1910 !! gsigt(k) = fssig (k ) 1911 !! gsigw(k) = fssig (k-0.5) 1912 !! This routine is given as an example, it must be modified 1913 !! following the user s desiderata. nevertheless, the output as 1914 !! well as the way to compute the model levels and scale factors 1915 !! must be respected in order to insure second order a!!uracy 1916 !! schemes. 1917 !! 1918 1919 !!====================================================================== 1920 INTEGER :: ji, jj, jk, jl, ik ! dummy loop argument 1921 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1922 INTEGER :: jpksigm ! temporary integer for maxnumber of s-levels 1923 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper,zrmin,e3t_t,e3w_t ! temporary scalars 1924 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1925 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 1926 REAL(wp) :: zmax, zmin ,zsigma ! Maximum and minimum depth and depth of sigma layer 1927 REAL(wp) :: zacr , zkth ,za1 ! parameters for z- layer (as ppacr , ppzkth ) 1928 1929 1930 ! 1931 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat ,zrpt 1932 1933 NAMELIST/namzgr_hyb/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, & 1934 ln_s_sigma, rn_bb, rn_hc,rn_zsigma , nn_sig_lev , rn_kth , rn_acr 1935 1936 1937 !!---------------------------------------------------------------------- 1938 ! 1939 IF( nn_timing == 1 ) CALL timing_start('zgr_hyb') 1940 ! 1941 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1942 1943 ! CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3 ) 1944 ! 1945 REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters 1946 READ ( numnam, namzgr_hyb ) 1947 IF(lwp) THEN ! control print 1948 WRITE(numout,*) 1949 WRITE(numout,*) 'dom:zgr_hyb : s-coordinate or hybrid z-s-coordinate' 1950 WRITE(numout,*) '~~~~~~~~~~~' 1951 WRITE(numout,*) ' Namelist namzgr_hyb' 1952 WRITE(numout,*) ' sigma-stretching coeffs ' 1953 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max 1954 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min 1955 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta 1956 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb 1957 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax 1958 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1959 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb 1960 WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc 1961 WRITE(numout,*) ' Sigma depth rn_zsigma = ', rn_zsigma 1962 WRITE(numout,*) ' The same as pp_arc rn_arc = ', rn_acr 1963 WRITE(numout,*) ' Number of sigma levels rn_arc = ', nn_sig_lev 1964 WRITE(numout,*) ' Number of levels for stretching z-coord rn_kth = ', rn_kth 1965 1966 1967 ENDIF 1968 zsigma = rn_zsigma 1969 jpksigm = nn_sig_lev 1970 zmax = rn_sbot_max 1971 zacr = rn_acr 1972 zkth = rn_kth 1973 e3t(:,:,:) = 1._wp 1974 e3w(:,:,:) = 1._wp 1975 e3u(:,:,:) = 1._wp 1976 e3v(:,:,:) = 1._wp 1977 e3f(:,:,:) = 1._wp 1978 e3uw(:,:,:)= 1._wp 1979 e3vw(:,:,:)= 1._wp 1980 1981 1982 1983 1984 1985 DO jj = 1, jpj 1986 DO ji= 1, jpi 1987 IF( bathy(ji,jj) <= 0._wp ) THEN ; bathy(ji,jj) = 0.e0_wp 1988 ELSE ; bathy(ji,jj) = MIN( rn_sbot_max, MAX( bathy(ji,jj),rn_sbot_min ) ) 1989 ENDIF 1990 END DO 1991 END DO 1992 1993 ! create bathymetry for enveloping 1160 1994 DO jj = 1, jpj 1161 1995 DO ji = 1, jpi 1162 1996 zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 1163 END DO1164 END DO1165 !1166 ! Smooth the bathymetry (if required)1997 zenv(ji,jj) = MIN (zenv(ji,jj), zsigma ) 1998 hbatt(ji,jj) = zenv(ji,jj) 1999 END DO 2000 END DO 1167 2001 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1168 2002 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1169 !1170 2003 jl = 0 1171 2004 zrmax = 1._wp … … 1197 2030 END DO 1198 2031 END DO 1199 ! 1200 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 2032 IF(lwp)WRITE(numout,*) 'zgr_hyb : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1201 2033 ! 1202 2034 DO jj = 1, nlcj … … 1222 2054 DO jj = 1, nlcj 1223 2055 DO ji = 1, nlci 1224 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )2056 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), hbatt(ji,jj) ) 1225 2057 END DO 1226 2058 END DO 1227 2059 ! 1228 ! Apply lateral boundary condition CAUTION: ke epthe value when the lbc field is zero2060 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 1229 2061 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1230 2062 DO jj = 1, nlcj … … 1237 2069 ! ! ================ ! 1238 2070 ! 1239 ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions 1240 DO ji = nlci+1, jpi 1241 zenv(ji,1:nlcj) = zenv(nlci,1:nlcj) 1242 END DO 1243 ! 1244 DO jj = nlcj+1, jpj 1245 zenv(:,jj) = zenv(:,nlcj) 1246 END DO 1247 ! 1248 ! Envelope bathymetry saved in hbatt 2071 ! ! envelop bathymetry saved in hbatt 1249 2072 hbatt(:,:) = zenv(:,:) 1250 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1251 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 1252 DO jj = 1, jpj 1253 DO ji = 1, jpi 1254 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 ) 1255 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 1256 END DO 1257 END DO 1258 ENDIF 1259 ! 1260 IF(lwp) THEN ! Control print 2073 IF( lk_mpp ) CALL mpp_max( nstop ) 2074 IF (lwp) write(numout,*)"after envelope", nstop 2075 2076 ! define new reference levels 2077 ! for s- levels, allows stretching at surface and bottom layer 2078 IF( jpksigm > 1 )THEN 2079 IF(ln_s_sigma)THEN 2080 DO jk = 1, jpksigm 2081 gsigw(jk) = -fssig3( REAL(jk,wp)-0.5_wp , rn_bb,jpksigm ) 2082 gsigt(jk) = -fssig3( REAL(jk,wp) , rn_bb,jpksigm ) 2083 END DO 2084 ELSE 2085 DO jk = 1, jpksigm 2086 gsigw(jk) = -fssig2( REAL(jk,wp)-0.5_wp ,jpksigm ) 2087 gsigt(jk) = -fssig2( REAL(jk,wp) ,jpksigm ) 2088 2089 END DO 2090 ENDIF 2091 gsigw(1)=0._wp ! set gsigw exactly to zero 2092 IF( lk_mpp ) CALL mpp_max( nstop ) 2093 IF (lwp) THEN 2094 write(numout,*)"after fssig", nstop,"gsigw,gsigt=" 2095 do jk=1,jpksigm 2096 write(numout,*)jk,gsigw(jk),gsigt(jk) 2097 enddo 2098 ENDIF 2099 2100 DO jk=1,jpksigm 2101 DO jj=1,jpj 2102 DO ji=1,jpi 2103 zrmin= min( hbatt(ji,jj), zsigma ) 2104 IF(hbatt(ji,jj).lt.rn_hc)THEN 2105 zcoefw=REAL(jk-1,wp) / REAL(jpksigm-1,wp) 2106 zcoeft=(REAL(jk-1,wp)+0.5)/ REAL(jpksigm-1,wp) 2107 ELSE 2108 zcoefw=gsigw(jk) 2109 zcoeft=gsigt(jk) 2110 2111 ENDIF 2112 2113 gdept(ji,jj,jk)=scosrf(ji,jj)+(zrmin-rn_hc)*zcoeft & 2114 +rn_hc* (REAL(jk,wp)- 0.5_wp) / REAL(jpksigm-1,wp) 2115 gdepw(ji,jj,jk)=scosrf(ji,jj)+(zrmin-rn_hc)*zcoefw & 2116 +rn_hc*( REAL(jk,wp)- 1.0_wp ) / REAL(jpksigm-1,wp) 2117 2118 ENDDO 2119 ENDDO 2120 ENDDO 2121 gdepw(:,:,1)= scosrf(:,:) 2122 ! redefine gdept_0, gdepw_0 which will be used in diawri.F90 2123 DO jk=1,jpksigm 2124 IF(zsigma.lt.rn_hc)THEN 2125 zcoefw=REAL(jk-1,wp) / REAL(jpksigm-1,wp) 2126 zcoeft=(REAL(jk-1,wp)+0.5)/ REAL(jpksigm-1,wp) 2127 ELSE 2128 zcoefw=gsigw(jk) 2129 zcoeft=gsigt(jk) 2130 ENDIF 2131 2132 gdept_0(jk)= zcoeft * (zsigma-rn_hc)+rn_hc* (REAL(jk,wp)- 0.5_wp )/REAL(jpksigm-1,wp) 2133 gdepw_0(jk)= zcoefw * (zsigma-rn_hc)+rn_hc* (REAL(jk,wp)- 1.0_wp )/REAL(jpksigm-1,wp) 2134 ENDDO 2135 2136 DO jk=1,jpksigm-1 2137 e3t_0(jk) = gdepw_0(jk+1)-gdepw_0(jk) 2138 e3w_0(jk+1)= gdept_0(jk+1)-gdept_0(jk) 2139 ENDDO 2140 e3w_0(1) = 2._wp * ( gdept_0(1 ) - gdepw_0(1 ) ) 2141 2142 2143 ! now for lower z-levels : 2144 zmin = e3t_0 (jpksigm -1 ) ! min layer width in z- zone is the same as lowest in s- layer 2145 ELSE 2146 zsigma = 0._wp 2147 hbatt(:,:) = zsigma 2148 zmin = 5._wp 2149 ENDIF 2150 IF(lwp) write(numout,*) ": last vertical level of sigma-coordinates",zmin 2151 CALL fszref ( zkth, zmin, zacr, zmax, jpksigm , zsigma ) 2152 2153 2154 DO jk=jpksigm,jpkm1 2155 2156 e3t_0(jk) = gdepw_0(jk+1)-gdepw_0(jk) 2157 e3w_0(jk+1)= gdept_0(jk+1)-gdept_0(jk) 2158 ENDDO 2159 e3w_0(1) = 2._wp * ( gdept_0(1 ) - gdepw_0(1 ) ) 2160 e3t_0(jpk) = 2._wp * ( gdept_0(jpk) - gdepw_0(jpk) ) 2161 2162 IF( lk_mpp ) CALL mpp_max( nstop ) 2163 IF (lwp) write(numout,*)"e3t0" ,nstop 2164 2165 IF(lwp) THEN ! control print 1261 2166 WRITE(numout,*) 1262 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 2167 WRITE(numout,*) ' zhyb Reference z-coordinate depth and scale factors:' 2168 WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) 2169 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk ) 2170 open(333,file='zmesh.dat') 2171 WRITE(333,*)'initial zsigma =', zsigma, jpk 2172 WRITE(333,*)'initial jpksigm =', jpksigm 2173 WRITE(333,*)'rn_bb =', rn_bb,'rn_theta=',rn_theta 2174 do jk=1,jpk 2175 WRITE(333,'(i4,1x,4(1x,e13.6))') jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk),e3w_0(jk) 2176 enddo 2177 close(333) 2178 ENDIF 2179 2180 DO jk=jpksigm,jpk 2181 DO jj=1,jpj 2182 DO ji=1,jpi 2183 IF(jpksigm>1)THEN 2184 gdept(ji,jj,jk)=scosrf(ji,jj) + gdept_0(jk) *hbatt(ji,jj)/zsigma ! differ from gdept0+scorf only at land 2185 gdepw(ji,jj,jk)=scosrf(ji,jj) + gdepw_0(jk) *hbatt(ji,jj)/zsigma ! as hbatt=zsigma over deep part of basin 2186 ELSE 2187 gdept(ji,jj,jk)=scosrf(ji,jj) + gdept_0(jk) 2188 gdepw(ji,jj,jk)=scosrf(ji,jj) + gdepw_0(jk) 2189 2190 ENDIF 2191 ENDDO 2192 ENDDO 2193 ENDDO 2194 2195 ! define e3t, e3w for general levels 2196 DO jk=1,jpkm1 2197 e3t(:,:,jk) = gdepw(:,:,jk+1)-gdepw(:,:,jk) 2198 e3w(:,:,jk+1)= gdept(:,:,jk+1)-gdept(:,:,jk) 2199 ENDDO 2200 e3w(:,:,1) = 2._wp * ( gdept(:,:,1 ) - gdepw(:,:,1 ) ) 2201 e3t(:,:,jpk) = 2._wp * ( gdept(:,:,jpk) - gdepw(:,:,jpk) ) 2202 2203 ! and surface : 2204 ! ! HYBRID mbathy : 2205 2206 IF(lwp) THEN ! control print 1263 2207 WRITE(numout,*) 1264 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout ) 1265 IF( nprint == 1 ) THEN 1266 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 1267 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 1268 ENDIF 1269 ENDIF 1270 1271 ! ! ============================== 1272 ! ! hbatu, hbatv, hbatf fields 1273 ! ! ============================== 1274 IF(lwp) THEN 1275 WRITE(numout,*) 1276 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 1277 ENDIF 1278 hbatu(:,:) = rn_sbot_min 1279 hbatv(:,:) = rn_sbot_min 1280 hbatf(:,:) = rn_sbot_min 1281 DO jj = 1, jpjm1 1282 DO ji = 1, jpim1 ! NO vector opt. 1283 hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 1284 hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 1285 hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 1286 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 1287 END DO 1288 END DO 1289 ! 1290 ! Apply lateral boundary condition 1291 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 1292 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) 2208 WRITE(numout,*) ' zhyb centre of basin s-z-coordinate depth and scale factors:' 2209 WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) 2210 write(numout,*)"bathy" ,"min e3t" 2211 do jk=1,jpk 2212 WRITE(numout, "(10x, i4, 4f9.2)" ) jk, gdept(20,20,jk), gdepw(20,20,jk), & 2213 & e3t(20,20,jk), e3w(20,20,jk) 2214 enddo 2215 ENDIF 2216 2217 mbathy(:,:)=0 2218 ! WHERE( 0._wp < bathy(:,:)) mbathy(:,:)=jpkm1 1293 2219 DO jj = 1, jpj 1294 2220 DO ji = 1, jpi 1295 IF( hbatu(ji,jj) == 0._wp ) THEN 1296 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 1297 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) 2221 DO jk = 1, jpkm1 2222 IF( bathy(ji,jj) >= gdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2223 2224 END DO 2225 END DO 2226 END DO 2227 2228 ! DO jk = jpkm1, jpksigm+1, -1 2229 ! zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) 2230 ! WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 2231 ! END DO 2232 2233 2234 2235 ! z-partial steps :goto 20 2236 DO jj = 1, jpj 2237 DO ji = 1, jpi 2238 ik = mbathy(ji,jj) 2239 IF( ik > jpksigm ) THEN ! ocean point only 2240 ! max ocean level case 2241 IF( ik == jpkm1 ) THEN 2242 zdepwp = bathy(ji,jj) 2243 ze3tp = bathy(ji,jj) - gdepw_0(ik) 2244 ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) ) 2245 e3t(ji,jj,ik ) = ze3tp 2246 e3t(ji,jj,ik+1) = ze3tp 2247 e3w(ji,jj,ik ) = ze3wp 2248 e3w(ji,jj,ik+1) = ze3tp 2249 gdepw(ji,jj,ik+1) = zdepwp 2250 gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp 2251 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 2252 ! 2253 ELSE ! standard case 2254 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj) 2255 ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 2256 ENDIF 2257 !gm Bug? check the gdepw_0 2258 ! ... on ik 2259 gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 2260 & * ((gdept_0( ik ) - gdepw_0(ik) ) & 2261 & / ( gdepw_0( ik+1) - gdepw_0(ik) )) 2262 e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 2263 & / ( gdepw_0( ik+1) - gdepw_0(ik) ) 2264 e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) & 2265 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 2266 ! ... on ik+1 2267 e3w (ji,jj,ik+1) = e3t (ji,jj,ik) 2268 e3t (ji,jj,ik+1) = e3t (ji,jj,ik) 2269 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 2270 ENDIF 1298 2271 ENDIF 1299 2272 END DO 1300 2273 END DO 1301 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) 2274 ! 2275 jl = 0 1302 2276 DO jj = 1, jpj 1303 2277 DO ji = 1, jpi 1304 IF( hbatv(ji,jj) == 0._wp ) THEN 1305 IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min 1306 IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) 2278 ik = mbathy(ji,jj) 2279 IF( ik > jpksigm ) THEN ! ocean point only 2280 e3tp (ji,jj) = e3t(ji,jj,ik ) 2281 e3wp (ji,jj) = e3w(ji,jj,ik ) 2282 ! test 2283 zmin= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) 2284 IF( zmin <= 0._wp .AND. lwp ) THEN 2285 jl = jl + 1 2286 WRITE(numout,*) ' it = ', jl, ' ik = ', ik, ' (i,j) = ', ji, jj 2287 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 2288 WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zmin 2289 WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik ) 2290 ENDIF 1307 2291 ENDIF 1308 2292 END DO 1309 2293 END DO 1310 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) 1311 DO jj = 1, jpj 1312 DO ji = 1, jpi 1313 IF( hbatf(ji,jj) == 0._wp ) THEN 1314 IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min 1315 IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) 1316 ENDIF 1317 END DO 1318 END DO 1319 1320 !!bug: key_helsinki a verifer 1321 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 1322 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 1323 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 1324 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 1325 1326 IF( nprint == 1 .AND. lwp ) THEN 1327 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 1328 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 1329 WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & 1330 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 1331 WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 1332 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 1333 WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 1334 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 1335 ENDIF 1336 !! helsinki 1337 1338 ! ! ======================= 1339 ! ! s-ccordinate fields (gdep., e3.) 1340 ! ! ======================= 1341 ! 1342 ! non-dimensional "sigma" for model level depth at w- and t-levels 1343 1344 1345 !======================================================================== 1346 ! Song and Haidvogel 1994 (ln_s_sh94=T) 1347 ! Siddorn and Furner 2012 (ln_sf12=T) 1348 ! or tanh function (both false) 1349 !======================================================================== 1350 IF ( ln_s_sh94 ) THEN 1351 CALL s_sh94() 1352 ELSE IF ( ln_s_sf12 ) THEN 1353 CALL s_sf12() 1354 ELSE 1355 CALL s_tanh() 1356 ENDIF 1357 1358 CALL lbc_lnk( e3t , 'T', 1._wp ) 1359 CALL lbc_lnk( e3u , 'U', 1._wp ) 1360 CALL lbc_lnk( e3v , 'V', 1._wp ) 1361 CALL lbc_lnk( e3f , 'F', 1._wp ) 1362 CALL lbc_lnk( e3w , 'W', 1._wp ) 1363 CALL lbc_lnk( e3uw, 'U', 1._wp ) 1364 CALL lbc_lnk( e3vw, 'V', 1._wp ) 1365 1366 fsdepw(:,:,:) = gdepw (:,:,:) 1367 fsde3w(:,:,:) = gdep3w(:,:,:) 1368 ! 1369 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0 1370 where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.0 1371 where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.0 1372 where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.0 1373 where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.0 1374 where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.0 1375 where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.0 1376 1377 1378 fsdept(:,:,:) = gdept (:,:,:) 1379 fsdepw(:,:,:) = gdepw (:,:,:) 1380 fsde3w(:,:,:) = gdep3w(:,:,:) 1381 fse3t (:,:,:) = e3t (:,:,:) 1382 fse3u (:,:,:) = e3u (:,:,:) 1383 fse3v (:,:,:) = e3v (:,:,:) 1384 fse3f (:,:,:) = e3f (:,:,:) 1385 fse3w (:,:,:) = e3w (:,:,:) 1386 fse3uw(:,:,:) = e3uw (:,:,:) 1387 fse3vw(:,:,:) = e3vw (:,:,:) 1388 !! 1389 ! HYBRID : 1390 DO jj = 1, jpj 1391 DO ji = 1, jpi 1392 DO jk = 1, jpkm1 1393 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1394 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2294 2295 ! Scale factors and depth at U-, V-, UW and VW-points 2296 DO jk = 1, jpk ! initialisation to z-scale factors 2297 e3u (:,:,jk) = e3t(:,:,jk) 2298 e3v (:,:,jk) = e3t(:,:,jk) 2299 e3uw(:,:,jk) = e3w(:,:,jk) 2300 e3vw(:,:,jk) = e3w(:,:,jk) 2301 e3f (:,:,jk) = e3t(:,:,jk) 2302 END DO 2303 DO jk = 1, jpksigm-1 2304 DO jj = 1, jpjm1 2305 DO ji = 1, fs_jpim1 2306 e3u(ji,jj,jk)=( REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk) + & 2307 REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3t(ji+1,jj,jk) ) & 2308 /MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) 2309 2310 e3uw(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3w(ji,jj,jk) + & 2311 REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3w(ji+1,jj,jk) ) & 2312 /REAL(MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj))),wp) 2313 2314 e3v(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk) + & 2315 REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3t(ji,jj+1,jk) ) & 2316 /REAL (MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))),wp) 2317 2318 e3vw(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3w(ji,jj,jk) + & 2319 REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3w(ji,jj+1,jk) ) & 2320 /REAL(MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1))),wp) 2321 2322 e3f(ji,jj,jk)=(REAL(MIN(1,mbathy(ji,jj)),wp)* e3t(ji,jj,jk) + & 2323 REAL(MIN(1,mbathy(ji+1,jj)),wp)*e3t(ji+1,jj,jk) + & 2324 REAL(MIN(1,mbathy(ji+1,jj+1)),wp)* e3t(ji+1,jj+1,jk)+ & 2325 REAL(MIN(1,mbathy(ji,jj+1)),wp)*e3t(ji,jj+1,jk) ) & 2326 /REAL(MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) & 2327 + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1))),wp) 2328 2329 ENDDO 1395 2330 END DO 1396 2331 END DO 1397 END DO 1398 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1399 & ' MAX ', MAXVAL( mbathy(:,:) ) 1400 1401 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 1402 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1403 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & 1404 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) ) 1405 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), & 1406 & ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), & 1407 & ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), & 1408 & ' w ', MINVAL( fse3w (:,:,:) ) 1409 1410 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & 1411 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) ) 1412 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), & 1413 & ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), & 1414 & ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), & 1415 & ' w ', MAXVAL( fse3w (:,:,:) ) 1416 ENDIF 1417 ! END DO 1418 IF(lwp) THEN ! selected vertical profiles 1419 WRITE(numout,*) 1420 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 1421 WRITE(numout,*) ' ~~~~~~ --------------------' 1422 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1423 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), & 1424 & fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk ) 1425 DO jj = mj0(20), mj1(20) 1426 DO ji = mi0(20), mi1(20) 1427 WRITE(numout,*) 1428 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1429 WRITE(numout,*) ' ~~~~~~ --------------------' 1430 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1431 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1432 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 2332 2333 DO jk = jpksigm,jpk ! Computed as the minimum of neighbooring scale factors 2334 DO jj = 1, jpjm1 2335 DO ji = 1, fs_jpim1 ! vector opt. 2336 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) ) 2337 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) ) 2338 e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 2339 e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 1433 2340 END DO 1434 2341 END DO 1435 DO jj = mj0(74), mj1(74) 1436 DO ji = mi0(100), mi1(100) 1437 WRITE(numout,*) 1438 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1439 WRITE(numout,*) ' ~~~~~~ --------------------' 1440 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1441 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1442 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 2342 END DO 2343 2344 CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions 2345 CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp ) 2346 ! 2347 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 2348 WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk) 2349 WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk) 2350 WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk) 2351 WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk) 2352 END DO 2353 2354 DO jk = jpksigm, jpk ! Computed as the minimum of neighbooring V-scale factors 2355 DO jj = 1, jpjm1 2356 DO ji = 1, fs_jpim1 ! vector opt. 2357 e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 1443 2358 END DO 1444 2359 END DO 1445 ENDIF 1446 1447 !================================================================================ 1448 ! check the coordinate makes sense 1449 !================================================================================ 1450 DO ji = 1, jpi 1451 DO jj = 1, jpj 1452 1453 IF( hbatt(ji,jj) > 0._wp) THEN 1454 DO jk = 1, mbathy(ji,jj) 1455 ! check coordinate is monotonically increasing 1456 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1457 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1458 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1459 WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 1460 WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 1461 CALL ctl_stop( ctmp1 ) 1462 ENDIF 1463 ! and check it has never gone negative 1464 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1465 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1466 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1467 WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1468 WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1469 CALL ctl_stop( ctmp1 ) 1470 ENDIF 1471 ! and check it never exceeds the total depth 1472 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 1473 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1474 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1475 WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1476 CALL ctl_stop( ctmp1 ) 1477 ENDIF 1478 END DO 1479 1480 DO jk = 1, mbathy(ji,jj)-1 1481 ! and check it never exceeds the total depth 1482 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 1483 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1484 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1485 WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1486 CALL ctl_stop( ctmp1 ) 1487 ENDIF 1488 END DO 1489 1490 ENDIF 1491 1492 END DO 1493 END DO 1494 ! 1495 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1496 ! 1497 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1498 ! 1499 END SUBROUTINE zgr_sco 1500 1501 !!====================================================================== 1502 SUBROUTINE s_sh94() 1503 1504 !!---------------------------------------------------------------------- 1505 !! *** ROUTINE s_sh94 *** 1506 !! 1507 !! ** Purpose : stretch the s-coordinate system 1508 !! 1509 !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 1510 !! mixed S/sigma coordinate 1511 !! 1512 !! Reference : Song and Haidvogel 1994. 1513 !!---------------------------------------------------------------------- 1514 ! 1515 INTEGER :: ji, jj, jk ! dummy loop argument 1516 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1517 ! 1518 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 1519 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1520 1521 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1522 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1523 1524 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 1525 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 1526 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 1527 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 1528 1529 DO ji = 1, jpi 1530 DO jj = 1, jpj 1531 1532 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1533 DO jk = 1, jpk 1534 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1535 z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1536 END DO 1537 ELSE ! shallow water, uniform sigma 1538 DO jk = 1, jpk 1539 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1540 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1541 END DO 1542 ENDIF 1543 ! 1544 DO jk = 1, jpkm1 1545 z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 1546 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 1547 END DO 1548 z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) ) 1549 z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 1550 ! 1551 ! Coefficients for vertical depth as the sum of e3w scale factors 1552 z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 1553 DO jk = 2, jpk 1554 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 1555 END DO 1556 ! 1557 DO jk = 1, jpk 1558 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1559 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1560 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1561 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1562 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1563 END DO 1564 ! 1565 END DO ! for all jj's 1566 END DO ! for all ji's 1567 1568 DO ji = 1, jpim1 1569 DO jj = 1, jpjm1 1570 DO jk = 1, jpk 1571 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 1572 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1573 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 1574 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1575 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 1576 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 1577 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1578 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 1579 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1580 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 1581 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1582 ! 1583 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1584 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1585 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1586 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1587 ! 1588 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1589 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1590 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1591 END DO 1592 END DO 1593 END DO 1594 1595 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1596 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1597 1598 END SUBROUTINE s_sh94 1599 1600 SUBROUTINE s_sf12 1601 1602 !!---------------------------------------------------------------------- 1603 !! *** ROUTINE s_sf12 *** 1604 !! 1605 !! ** Purpose : stretch the s-coordinate system 1606 !! 1607 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 1608 !! mixed S/sigma/Z coordinate 1609 !! 1610 !! This method allows the maintenance of fixed surface and or 1611 !! bottom cell resolutions (cf. geopotential coordinates) 1612 !! within an analytically derived stretched S-coordinate framework. 1613 !! 1614 !! 1615 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 1616 !!---------------------------------------------------------------------- 1617 ! 1618 INTEGER :: ji, jj, jk ! dummy loop argument 1619 REAL(wp) :: zsmth ! smoothing around critical depth 1620 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 1621 ! 1622 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 1623 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 1624 1625 ! 1626 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1627 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1628 1629 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 1630 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 1631 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 1632 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 1633 1634 DO ji = 1, jpi 1635 DO jj = 1, jpj 1636 1637 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 1638 1639 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 1640 ! could be changed by users but care must be taken to do so carefully 1641 zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 1642 1643 zzs = rn_zs / hbatt(ji,jj) 1644 1645 IF (rn_efold /= 0.0_wp) THEN 1646 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 1647 ELSE 1648 zsmth = 1.0_wp 1649 ENDIF 1650 1651 DO jk = 1, jpk 1652 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1653 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 1654 ENDDO 1655 z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) 1656 z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) 1657 1658 ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 1659 1660 DO jk = 1, jpk 1661 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1662 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 1663 END DO 1664 1665 ELSE ! shallow water, z coordinates 1666 1667 DO jk = 1, jpk 1668 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 1669 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 1670 END DO 1671 1672 ENDIF 1673 1674 DO jk = 1, jpkm1 1675 z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 1676 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 1677 END DO 1678 z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 )) 1679 z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 1680 1681 ! Coefficients for vertical depth as the sum of e3w scale factors 1682 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 1683 DO jk = 2, jpk 1684 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 1685 END DO 1686 1687 DO jk = 1, jpk 1688 gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 1689 gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 1690 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 1691 END DO 1692 1693 ENDDO ! for all jj's 1694 ENDDO ! for all ji's 1695 1696 DO ji=1,jpi-1 1697 DO jj=1,jpj-1 1698 1699 DO jk = 1, jpk 1700 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 1701 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1702 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 1703 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1704 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 1705 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 1706 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1707 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 1708 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1709 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 1710 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1711 1712 e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 1713 e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 1714 e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 1715 e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 1716 ! 1717 e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 1718 e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 1719 e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 1720 END DO 1721 1722 ENDDO 1723 ENDDO 1724 ! 1725 CALL lbc_lnk(e3t ,'T',1.) ; CALL lbc_lnk(e3u ,'T',1.) 1726 CALL lbc_lnk(e3v ,'T',1.) ; CALL lbc_lnk(e3f ,'T',1.) 1727 CALL lbc_lnk(e3w ,'T',1.) 1728 CALL lbc_lnk(e3uw,'T',1.) ; CALL lbc_lnk(e3vw,'T',1.) 1729 ! 1730 ! ! ============= 1731 1732 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 1733 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 1734 1735 END SUBROUTINE s_sf12 1736 1737 SUBROUTINE s_tanh() 1738 1739 !!---------------------------------------------------------------------- 1740 !! *** ROUTINE s_tanh*** 1741 !! 1742 !! ** Purpose : stretch the s-coordinate system 1743 !! 1744 !! ** Method : s-coordinate stretch 1745 !! 1746 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1747 !!---------------------------------------------------------------------- 1748 1749 INTEGER :: ji, jj, jk ! dummy loop argument 1750 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1751 1752 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 1753 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 1754 1755 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 1756 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 1757 1758 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp 1759 z_esigt = 0._wp ; z_esigw = 0._wp 1760 1761 DO jk = 1, jpk 1762 z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1763 z_gsigt(jk) = -fssig( REAL(jk,wp) ) 1764 END DO 1765 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk) 1766 ! 1767 ! Coefficients for vertical scale factors at w-, t- levels 1768 !!gm bug : define it from analytical function, not like juste bellow.... 1769 !!gm or betteroffer the 2 possibilities.... 1770 DO jk = 1, jpkm1 1771 z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk) 1772 z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 1773 END DO 1774 z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) ) 1775 z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 1776 ! 1777 ! Coefficients for vertical depth as the sum of e3w scale factors 1778 z_gsi3w(1) = 0.5_wp * z_esigw(1) 2360 END DO 2361 CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions 2362 ! 2363 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 2364 WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk) 2365 END DO 2366 !!gm bug ? : must be a do loop with mj0,mj1 2367 ! 2368 e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 2369 e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 2370 e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 2371 e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 2372 e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 2373 2374 ! Control of the sign 2375 Do jk=1,jpk 2376 do jj=1,jpj 2377 do ji=1,jpi 2378 IF( ( e3t (ji,jj,jk) ) <= 0._wp )then 2379 write(numout,*)' zgr_hyb : e r r o r e3t <= 0',ji,jj,jk,e3t (ji,jj,jk); endif 2380 IF( ( e3w (ji,jj,jk) ) <= 0._wp )then 2381 write(numout,*)' zgr_hyb : e r r o r e3t <= 0',ji,jj,jk,e3w (ji,jj,jk); endif 2382 2383 2384 IF( ( gdept(ji,jj,jk) ) < 0._wp )THEN 2385 write (numout,*)' zgr_hyb : e r r o r gdept < 0',ji,jj,jk ,gdept(ji,jj,jj); endif 2386 IF( ( gdepw(ji,jj,jk) ) < 0._wp )then 2387 write (numout,*)' zgr_hyb : e r r o r gdepw < 0',ji,jj,jk , gdepw(ji,jj,jj); endif 2388 enddo 2389 enddo 2390 enddo 2391 2392 ! Compute gdep3w (vertical sum of e3w) 2393 gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1) 1779 2394 DO jk = 2, jpk 1780 z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 1781 END DO 1782 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1783 DO jk = 1, jpk 1784 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1785 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1786 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 1787 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 1788 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 1789 END DO 1790 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1791 DO jj = 1, jpj 1792 DO ji = 1, jpi 1793 DO jk = 1, jpk 1794 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1795 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1796 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1797 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1798 ! 1799 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1800 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1801 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1802 END DO 1803 END DO 1804 END DO 1805 1806 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 1807 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 1808 1809 END SUBROUTINE s_tanh 1810 1811 FUNCTION fssig( pk ) RESULT( pf ) 1812 !!---------------------------------------------------------------------- 1813 !! *** ROUTINE fssig *** 1814 !! 1815 !! ** Purpose : provide the analytical function in s-coordinate 1816 !! 1817 !! ** Method : the function provide the non-dimensional position of 1818 !! T and W (i.e. between 0 and 1) 1819 !! T-points at integer values (between 1 and jpk) 1820 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1821 !!---------------------------------------------------------------------- 1822 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 1823 REAL(wp) :: pf ! sigma value 1824 !!---------------------------------------------------------------------- 1825 ! 1826 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 1827 & - TANH( rn_thetb * rn_theta ) ) & 1828 & * ( COSH( rn_theta ) & 1829 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1830 & / ( 2._wp * SINH( rn_theta ) ) 1831 ! 1832 END FUNCTION fssig 1833 1834 1835 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 1836 !!---------------------------------------------------------------------- 1837 !! *** ROUTINE fssig1 *** 1838 !! 1839 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 1840 !! 1841 !! ** Method : the function provides the non-dimensional position of 1842 !! T and W (i.e. between 0 and 1) 1843 !! T-points at integer values (between 1 and jpk) 1844 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1845 !!---------------------------------------------------------------------- 1846 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 1847 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 1848 REAL(wp) :: pf1 ! sigma value 1849 !!---------------------------------------------------------------------- 1850 ! 1851 IF ( rn_theta == 0 ) then ! uniform sigma 1852 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 1853 ELSE ! stretched sigma 1854 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 1855 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1856 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1857 ENDIF 1858 ! 1859 END FUNCTION fssig1 1860 1861 1862 FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) 1863 !!---------------------------------------------------------------------- 1864 !! *** ROUTINE fgamma *** 1865 !! 1866 !! ** Purpose : provide analytical function for the s-coordinate 1867 !! 1868 !! ** Method : the function provides the non-dimensional position of 1869 !! T and W (i.e. between 0 and 1) 1870 !! T-points at integer values (between 1 and jpk) 1871 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1872 !! 1873 !! This method allows the maintenance of fixed surface and or 1874 !! bottom cell resolutions (cf. geopotential coordinates) 1875 !! within an analytically derived stretched S-coordinate framework. 1876 !! 1877 !! Reference : Siddorn and Furner, in prep 1878 !!---------------------------------------------------------------------- 1879 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 1880 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 1881 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 1882 REAL(wp), INTENT(in ) :: pzs ! surface box depth 1883 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 1884 REAL(wp) :: za1,za2,za3 ! local variables 1885 REAL(wp) :: zn1,zn2 ! local variables 1886 REAL(wp) :: za,zb,zx ! local variables 1887 integer :: jk 1888 !!---------------------------------------------------------------------- 1889 ! 1890 1891 zn1 = 1./(jpk-1.) 1892 zn2 = 1. - zn1 1893 1894 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 1895 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 1896 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2395 gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 2396 END DO 2397 2398 2399 IF( lk_mpp ) CALL mpp_max( nstop ) 2400 IF (lwp) write(numout,*)"zpartial" ,nstop 2401 2402 2403 2404 CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions 2405 1897 2406 1898 za = pzb - za3*(pzs-za1)-za2 1899 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 1900 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 1901 zx = 1.0_wp-za/2.0_wp-zb 1902 1903 DO jk = 1, jpk 1904 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 1905 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 1906 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 1907 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 1908 ENDDO 1909 1910 ! 1911 END FUNCTION fgamma 2407 2408 2409 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 2410 ! 2411 IF( nn_timing == 1 ) CALL timing_stop('zgr_hyb') 2412 2413 END SUBROUTINE zgr_hyb 2414 1912 2415 1913 2416 !!====================================================================== -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
r2528 r6736 32 32 # define fse3vw(i,j,k) e3vw_1(i,j,k) 33 33 34 #if defined key_jth_fix 35 # define fsdept_b(i,j,k) (fsdept_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 36 # define fsdepw_b(i,j,k) (fsdepw_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 37 # define fsde3w_b(i,j,k) (fsde3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))-sshb(i,j)) 38 # define fse3t_b(i,j,k) (fse3t_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 39 # define fse3u_b(i,j,k) (fse3u_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 40 # define fse3v_b(i,j,k) (fse3v_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 41 # define fse3f_b(i,j,k) (fse3f_0(i,j,k)*(1.+sshf_b(i,j)*muf(i,j,k))) 42 # define fse3w_b(i,j,k) (fse3w_0(i,j,k)*(1.+sshb(i,j)*mut(i,j,k))) 43 # define fse3uw_b(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 44 # define fse3vw_b(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 45 #else 34 46 # define fse3t_b(i,j,k) e3t_b(i,j,k) 35 47 # define fse3u_b(i,j,k) e3u_b(i,j,k) … … 37 49 # define fse3uw_b(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k))) 38 50 # define fse3vw_b(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k))) 51 #endif 39 52 40 53 # define fsdept_n(i,j,k) (fsdept_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))) … … 51 64 # define fse3t_m(i,j,k) (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k))) 52 65 66 #if defined key_jth_fix 67 # define fsdept_a(i,j,k) (fsdept_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 68 # define fsdepw_a(i,j,k) (fsdepw_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 69 # define fsde3w_a(i,j,k) (fsde3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))-ssha(i,j)) 70 #endif 53 71 # define fse3t_a(i,j,k) (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 54 72 # define fse3u_a(i,j,k) (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 55 73 # define fse3v_a(i,j,k) (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 74 #if defined key_jth_fix 75 # define fse3f_a(i,j,k) (fse3f_0(i,j,k)*(1.+sshf_a(i,j)*muf(i,j,k))) 76 # define fse3w_a(i,j,k) (fse3w_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 77 # define fse3uw_a(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 78 # define fse3vw_a(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 79 #endif 56 80 57 81 #else … … 68 92 # define fse3vw(i,j,k) fse3vw_0(i,j,k) 69 93 94 #if defined key_jth_fix 95 # define fsdept_b(i,j,k) fsdept_0(i,j,k) 96 # define fsdepw_b(i,j,k) fsdepw_0(i,j,k) 97 # define fsde3w_b(i,j,k) fsde3w_0(i,j,k) 98 #endif 70 99 # define fse3t_b(i,j,k) fse3t_0(i,j,k) 71 100 # define fse3u_b(i,j,k) fse3u_0(i,j,k) 72 101 # define fse3v_b(i,j,k) fse3v_0(i,j,k) 102 #if defined key_jth_fix 103 # define fse3f_b(i,j,k) fse3f_0(i,j,k) 104 # define fse3w_b(i,j,k) fse3w_0(i,j,k) 105 #endif 73 106 # define fse3uw_b(i,j,k) fse3uw_0(i,j,k) 74 107 # define fse3vw_b(i,j,k) fse3vw_0(i,j,k) … … 87 120 # define fse3t_m(i,j,k) fse3t_0(i,j,k) 88 121 122 #if defined key_jth_fix 123 # define fsdept_a(i,j,k) fsdept_0(i,j,k) 124 # define fsdepw_a(i,j,k) fsdepw_0(i,j,k) 125 # define fsde3w_a(i,j,k) fsde3w_0(i,j,k) 126 #endif 89 127 # define fse3t_a(i,j,k) fse3t_0(i,j,k) 90 128 # define fse3u_a(i,j,k) fse3u_0(i,j,k) 91 129 # define fse3v_a(i,j,k) fse3v_0(i,j,k) 130 #if defined key_jth_fix 131 # define fse3f_a(i,j,k) fse3f_0(i,j,k) 132 # define fse3w_a(i,j,k) fse3w_0(i,j,k) 133 # define fse3uw_a(i,j,k) fse3uw_0(i,j,k) 134 # define fse3vw_a(i,j,k) fse3vw_0(i,j,k) 135 #endif 92 136 #endif 93 137 !!---------------------------------------------------------------------- -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r3294 r6736 18 18 USE dom_oce ! ocean space and time domain 19 19 USE fldread ! read input fields 20 USE in_out_manager ! I/O manager21 20 USE phycst ! physical constants 22 21 USE lib_mpp ! MPP library 23 22 USE wrk_nemo ! Memory allocation 24 23 USE timing ! Timing 24 USE in_out_manager ! I/O manager 25 USE iom 25 26 26 27 IMPLICIT NONE … … 34 35 35 36 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read) 37 #if defined key_jdha_init 38 REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) :: gdept_init 39 #endif 36 40 37 41 !! * Substitutions … … 146 150 INTEGER , INTENT(in ) :: kt ! ocean time-step 147 151 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: ptsd ! T & S data 152 ! REAL(wp), DIMENSION(jpi,jpj,jpk) :: gdept_init ! T & S data 148 153 ! 149 154 INTEGER :: ji, jj, jk, jl, jkk ! dummy loop indicies … … 156 161 ! 157 162 CALL fld_read( kt, 1, sf_tsd ) !== read T & S data at kt time step ==! 163 164 #if defined key_jdha_init 165 ALLOCATE( gdept_init(jpi,jpj,jpk) ) 166 CALL iom_open ( sf_tsd(jp_tem)%clname, sf_tsd(jp_tem)%num ) 167 CALL iom_get ( sf_tsd(jp_tem)%num, jpdom_data, 'deptht', gdept_init,1) 168 CALL iom_close( sf_tsd(jp_tem)%num ) ! Close the input file 169 #endif 158 170 ! 159 171 ! … … 223 235 DO jk = 1, jpk ! determines the intepolated T-S profiles at each (i,j) points 224 236 zl = fsdept_0(ji,jj,jk) 237 #if defined key_jdha_init 238 IF( zl < gdept_init(ji,jj,1 ) ) THEN ! above the first level of data 239 #else 225 240 IF( zl < gdept_0(1 ) ) THEN ! above the first level of data 241 #endif 226 242 ztp(jk) = ptsd(ji,jj,1 ,jp_tem) 227 243 zsp(jk) = ptsd(ji,jj,1 ,jp_sal) 244 #if defined key_jdha_init 245 ELSEIF( zl > gdept_init(ji,jj,jpk) ) THEN ! below the last level of data 246 #else 228 247 ELSEIF( zl > gdept_0(jpk) ) THEN ! below the last level of data 248 #endif 229 249 ztp(jk) = ptsd(ji,jj,jpkm1,jp_tem) 230 250 zsp(jk) = ptsd(ji,jj,jpkm1,jp_sal) 231 251 ELSE ! inbetween : vertical interpolation between jkk & jkk+1 232 252 DO jkk = 1, jpkm1 ! when gdept(jkk) < zl < gdept(jkk+1) 253 #if defined key_jdha_init 254 IF( (zl-gdept_init(ji,jj,jkk)) * (zl-gdept_init(ji,jj,jkk+1)) <= 0._wp ) THEN 255 zi = ( zl - gdept_init(ji,jj,jkk) ) / (gdept_init(ji,jj,jkk+1)-gdept_init(ji,jj,jkk)) 256 #else 233 257 IF( (zl-gdept_0(jkk)) * (zl-gdept_0(jkk+1)) <= 0._wp ) THEN 234 258 zi = ( zl - gdept_0(jkk) ) / (gdept_0(jkk+1)-gdept_0(jkk)) 259 #endif 235 260 ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 236 261 zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi … … 299 324 IF( sf_tsd(jp_sal)%ln_tint ) DEALLOCATE( sf_tsd(jp_sal)%fdta ) 300 325 DEALLOCATE( sf_tsd ) ! the structure itself 326 #if defined key_jdha_init 327 DEALLOCATE( gdept_init ) 328 #endif 301 329 ENDIF 302 330 ! -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r3764 r6736 32 32 USE phycst ! physical constants 33 33 USE dtatsd ! data temperature and salinity (dta_tsd routine) 34 USE restart ! ocean restart (rst_read routine) 34 35 USE in_out_manager ! I/O manager 35 36 USE iom ! I/O library … … 43 44 USE sol_oce ! ocean solver variables 44 45 USE lib_mpp ! MPP library 45 USE restart ! restart46 46 USE wrk_nemo ! Memory allocation 47 47 USE timing ! Timing … … 70 70 ! - ML - needed for initialization of e3t_b 71 71 INTEGER :: jk ! dummy loop indice 72 INTEGER :: inum ! temporary logical unit 72 73 !!---------------------------------------------------------------------- 73 74 ! … … 91 92 CALL rst_read ! Read the restart file 92 93 ! ! define e3u_b, e3v_b from e3t_b read in restart file 94 #if ! defined key_jth_fix 93 95 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 96 #endif 94 97 CALL day_init ! model calendar (using both namelist and restart infos) 95 98 ELSE … … 114 117 CALL dta_tsd( nit000, tsb ) ! read 3D T and S data at nit000 115 118 tsn(:,:,:,:) = tsb(:,:,:,:) 119 #if defined key_jdha_ssh_init 120 CALL iom_open ( 'initcd_ssh.nc', inum ) 121 CALL iom_get ( inum, jpdom_data, 'sossheig', sshb(:,:)) 122 CALL iom_close( inum ) ! Close the input file 123 sshn(:,:) = sshb(:,:) 124 #endif 116 125 ! 117 126 ELSE ! Initial T-S fields defined analytically … … 126 135 ! 127 136 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 137 #if ! defined key_jth_fix 128 138 IF( lk_vvl ) THEN 129 139 DO jk = 1, jpk … … 133 143 ! ! define e3u_b, e3v_b from e3t_b initialized in domzgr 134 144 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 145 #endif 135 146 ! 136 147 ENDIF … … 164 175 INTEGER :: ji, jj, jk 165 176 REAL(wp) :: zsal = 35.50 177 #if defined key_istate_fixed 178 REAL(wp) :: ztem = 25.50 179 #endif 166 180 !!---------------------------------------------------------------------- 167 181 ! … … 170 184 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)' 171 185 ! 186 #if ! defined key_istate_fixed 172 187 DO jk = 1, jpk 173 188 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) & … … 175 190 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 176 191 END DO 192 #else 193 tsn(:,:,:,jp_tem) = ztem * tmask(:,:,:) 194 tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 195 #endif 177 196 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 178 197 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) -
branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r3625 r6736 27 27 REAL(wp), PUBLIC :: rsmall = 0.5 * EPSILON( 1.e0 ) !: smallest real computer value 28 28 29 REAL(wp), PUBLIC :: rday = 24.*60.*60. !: day [s] 30 REAL(wp), PUBLIC :: rsiyea !: sideral year [s] 31 REAL(wp), PUBLIC :: rsiday !: sideral day [s] 32 REAL(wp), PUBLIC :: raamo = 12._wp !: number of months in one year 33 REAL(wp), PUBLIC :: rjjhh = 24._wp !: number of hours in one day 34 REAL(wp), PUBLIC :: rhhmm = 60._wp !: number of minutes in one hour 35 REAL(wp), PUBLIC :: rmmss = 60._wp !: number of seconds in one minute 36 REAL(wp), PUBLIC :: omega !: earth rotation parameter [s-1] 37 REAL(wp), PUBLIC :: ra = 6371229._wp !: earth radius [m] 38 REAL(wp), PUBLIC :: grav = 9.80665_wp !: gravity [m/s2] 29 REAL(wp), PUBLIC :: rday = 24.*60.*60. !: day (s) 30 REAL(wp), PUBLIC :: rsiyea !: sideral year (s) 31 REAL(wp), PUBLIC :: rsiday !: sideral day (s) 32 REAL(wp), PUBLIC :: raamo = 12._wp !: number of months in one year 33 REAL(wp), PUBLIC :: rjjhh = 24._wp !: number of hours in one day 34 REAL(wp), PUBLIC :: rhhmm = 60._wp !: number of minutes in one hour 35 REAL(wp), PUBLIC :: rmmss = 60._wp !: number of seconds in one minute 36 !! REAL(wp), PUBLIC :: omega = 7.292115083046061e-5_wp , & !: change the last digit! 37 REAL(wp), PUBLIC :: omega !: earth rotation parameter 38 REAL(wp), PUBLIC :: ra = 6371229._wp !: earth radius (meter) 39 REAL(wp), PUBLIC :: grav = 9.80665_wp !: gravity (m/s2) 39 40 40 REAL(wp), PUBLIC :: rtt = 273.16_wp !: triple point of temperature [Kelvin]41 REAL(wp), PUBLIC :: rt0 = 273.15_wp !: freezing point of fresh water [Kelvin]41 REAL(wp), PUBLIC :: rtt = 273.16_wp !: triple point of temperature (Kelvin) 42 REAL(wp), PUBLIC :: rt0 = 273.15_wp !: freezing point of water (Kelvin) 42 43 #if defined key_lim3 43 REAL(wp), PUBLIC :: rt0_snow = 273.16_wp !: melting point of snow [Kelvin]44 REAL(wp), PUBLIC :: rt0_ice = 273.16_wp !: melting point of ice [Kelvin]44 REAL(wp), PUBLIC :: rt0_snow = 273.16_wp !: melting point of snow (Kelvin) 45 REAL(wp), PUBLIC :: rt0_ice = 273.16_wp !: melting point of ice (Kelvin) 45 46 #else 46 REAL(wp), PUBLIC :: rt0_snow = 273.15_wp !: melting point of snow [Kelvin]47 REAL(wp), PUBLIC :: rt0_ice = 273.05_wp !: melting point of ice [Kelvin]47 REAL(wp), PUBLIC :: rt0_snow = 273.15_wp !: melting point of snow (Kelvin) 48 REAL(wp), PUBLIC :: rt0_ice = 273.05_wp !: melting point of ice (Kelvin) 48 49 #endif 50 49 51 #if defined key_cice 50 REAL(wp), PUBLIC :: rau0 = 1026._wp !: volumic mass of reference [kg/m3]52 REAL(wp), PUBLIC :: rau0 = 1026._wp !: reference volumic mass (density) (kg/m3) 51 53 #else 52 REAL(wp), PUBLIC :: rau0 = 1035._wp !: volumic mass of reference [kg/m3]54 REAL(wp), PUBLIC :: rau0 = 1035._wp !: reference volumic mass (density) (kg/m3) 53 55 #endif 54 REAL(wp), PUBLIC :: r1_rau0 !: = 1. / rau0 [m3/kg] 55 REAL(wp), PUBLIC :: rauw = 1000._wp !: volumic mass of pure water [m3/kg] 56 REAL(wp), PUBLIC :: rcp = 4.e3_wp !: ocean specific heat [J/Kelvin] 57 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 58 REAL(wp), PUBLIC :: r1_rau0_rcp !: = 1. / ( rau0 * rcp ) 59 60 REAL(wp), PUBLIC :: rhosn = 330._wp !: volumic mass of snow [kg/m3] 61 REAL(wp), PUBLIC :: emic = 0.97_wp !: emissivity of snow or ice 62 REAL(wp), PUBLIC :: sice = 6.0_wp !: salinity of ice [psu] 63 REAL(wp), PUBLIC :: soce = 34.7_wp !: salinity of sea [psu] 64 REAL(wp), PUBLIC :: cevap = 2.5e+6_wp !: latent heat of evaporation (water) 65 REAL(wp), PUBLIC :: srgamma = 0.9_wp !: correction factor for solar radiation (Oberhuber, 1974) 66 REAL(wp), PUBLIC :: vkarmn = 0.4_wp !: von Karman constant 67 REAL(wp), PUBLIC :: stefan = 5.67e-8_wp !: Stefan-Boltzmann constant 56 REAL(wp), PUBLIC :: rau0r !: reference specific volume (m3/kg) 57 REAL(wp), PUBLIC :: rcp = 4.e+3_wp !: ocean specific heat 58 REAL(wp), PUBLIC :: ro0cpr !: = 1. / ( rau0 * rcp ) 68 59 69 60 #if defined key_lim3 || defined key_cice 70 REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice [kg/m3] 71 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice 72 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow 73 REAL(wp), PUBLIC :: cpic = 2067.0_wp !: specific heat for ice 74 REAL(wp), PUBLIC :: lsub = 2.834e+6_wp !: pure ice latent heat of sublimation [J/kg] 75 REAL(wp), PUBLIC :: lfus = 0.334e+6_wp !: latent heat of fusion of fresh ice [J/kg] 76 REAL(wp), PUBLIC :: tmut = 0.054_wp !: decrease of seawater meltpoint with salinity 77 REAL(wp), PUBLIC :: xlsn !: = lfus*rhosn (volumetric latent heat fusion of snow) [J/m3] 61 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow 62 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice 63 REAL(wp), PUBLIC :: cpic = 2067.0 !: specific heat of sea ice 64 REAL(wp), PUBLIC :: lsub = 2.834e+6 !: pure ice latent heat of sublimation (J.kg-1) 65 REAL(wp), PUBLIC :: lfus = 0.334e+6 !: latent heat of fusion of fresh ice (J.kg-1) 66 REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice (kg/m3) 67 REAL(wp), PUBLIC :: tmut = 0.054 !: decrease of seawater meltpoint with salinity 78 68 #else 79 REAL(wp), PUBLIC :: rhoic = 900._wp !: volumic mass of sea ice [kg/m3] 80 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: conductivity of the ice [W/m/K] 81 REAL(wp), PUBLIC :: rcpic = 1.8837e+6_wp !: volumetric specific heat for ice [J/m3/K] 82 REAL(wp), PUBLIC :: cpic !: = rcpic / rhoic (specific heat for ice) [J/Kg/K] 83 REAL(wp), PUBLIC :: rcdsn = 0.22_wp !: conductivity of the snow [W/m/K] 84 REAL(wp), PUBLIC :: rcpsn = 6.9069e+5_wp !: volumetric specific heat for snow [J/m3/K] 85 REAL(wp), PUBLIC :: xlsn = 110.121e+6_wp !: volumetric latent heat fusion of snow [J/m3] 86 REAL(wp), PUBLIC :: lfus !: = xlsn / rhosn (latent heat of fusion of fresh ice) [J/Kg] 87 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice [J/m3] 88 REAL(wp), PUBLIC :: xsn = 2.8e+6_wp !: volumetric latent heat of sublimation of snow [J/m3] 69 REAL(wp), PUBLIC :: rcdsn = 0.22_wp !: conductivity of the snow 70 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: conductivity of the ice 71 REAL(wp), PUBLIC :: rcpsn = 6.9069e+5_wp !: density times specific heat for snow 72 REAL(wp), PUBLIC :: rcpic = 1.8837e+6_wp !: volumetric latent heat fusion of sea ice 73 REAL(wp), PUBLIC :: lfus = 0.3337e+6 !: latent heat of fusion of fresh ice (J.kg-1) 74 REAL(wp), PUBLIC :: xlsn = 110.121e+6_wp !: volumetric latent heat fusion of snow 75 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice 76 REAL(wp), PUBLIC :: xsn = 2.8e+6 !: latent heat of sublimation of snow 77 REAL(wp), PUBLIC :: rhoic = 900._wp !: volumic mass of sea ice (kg/m3) 89 78 #endif 79 REAL(wp), PUBLIC :: rhosn = 330._wp !: volumic mass of snow (kg/m3) 80 REAL(wp), PUBLIC :: emic = 0.97_wp !: emissivity of snow or ice 81 REAL(wp), PUBLIC :: sice = 6.0_wp !: reference salinity of ice (psu) 82 REAL(wp), PUBLIC :: soce = 34.7_wp !: reference salinity of sea (psu) 83 REAL(wp), PUBLIC :: cevap = 2.5e+6_wp !: latent heat of evaporation (water) 84 REAL(wp), PUBLIC :: srgamma = 0.9_wp !: correction factor for solar radiation (Oberhuber, 1974) 85 REAL(wp), PUBLIC :: vkarmn = 0.4_wp !: von Karman constant 86 REAL(wp), PUBLIC :: stefan = 5.67e-8_wp !: Stefan-Boltzmann constant 90 87 !!---------------------------------------------------------------------- 91 88 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 105 102 !!---------------------------------------------------------------------- 106 103 107 IF(lwp) WRITE(numout,*) 108 IF(lwp) WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants' 109 IF(lwp) WRITE(numout,*) ' ~~~~~~~' 104 ! ! Define additional parameters 105 rsiyea = 365.25 * rday * 2. * rpi / 6.283076 106 rsiday = rday / ( 1. + rday / rsiyea ) 107 #if defined key_cice 108 omega = 7.292116e-05 109 #else 110 omega = 2. * rpi / rsiday 111 #endif 110 112 111 ! Ocean Parameters 112 ! ---------------- 113 IF(lwp) THEN 113 rau0r = 1. / rau0 114 ro0cpr = 1. / ( rau0 * rcp ) 115 116 117 IF(lwp) THEN ! control print 118 WRITE(numout,*) 119 WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants' 120 WRITE(numout,*) ' ~~~~~~~' 114 121 WRITE(numout,*) ' Domain info' 115 122 WRITE(numout,*) ' dimension of model' … … 124 131 WRITE(numout,*) ' jpnij : ', jpnij 125 132 WRITE(numout,*) ' lateral domain boundary condition type : jperio = ', jperio 126 ENDIF 127 128 ! Define constants 129 ! ---------------- 130 IF(lwp) WRITE(numout,*) 131 IF(lwp) WRITE(numout,*) ' Constants' 132 133 IF(lwp) WRITE(numout,*) 134 IF(lwp) WRITE(numout,*) ' mathematical constant rpi = ', rpi 135 136 rsiyea = 365.25_wp * rday * 2._wp * rpi / 6.283076_wp 137 rsiday = rday / ( 1._wp + rday / rsiyea ) 138 #if defined key_cice 139 omega = 7.292116e-05 140 #else 141 omega = 2._wp * rpi / rsiday 142 #endif 143 IF(lwp) WRITE(numout,*) 144 IF(lwp) WRITE(numout,*) ' day rday = ', rday, ' s' 145 IF(lwp) WRITE(numout,*) ' sideral year rsiyea = ', rsiyea, ' s' 146 IF(lwp) WRITE(numout,*) ' sideral day rsiday = ', rsiday, ' s' 147 IF(lwp) WRITE(numout,*) ' omega omega = ', omega, ' s^-1' 148 149 IF(lwp) WRITE(numout,*) 150 IF(lwp) WRITE(numout,*) ' nb of months per year raamo = ', raamo, ' months' 151 IF(lwp) WRITE(numout,*) ' nb of hours per day rjjhh = ', rjjhh, ' hours' 152 IF(lwp) WRITE(numout,*) ' nb of minutes per hour rhhmm = ', rhhmm, ' mn' 153 IF(lwp) WRITE(numout,*) ' nb of seconds per minute rmmss = ', rmmss, ' s' 154 155 IF(lwp) WRITE(numout,*) 156 IF(lwp) WRITE(numout,*) ' earth radius ra = ', ra, ' m' 157 IF(lwp) WRITE(numout,*) ' gravity grav = ', grav , ' m/s^2' 158 159 IF(lwp) WRITE(numout,*) 160 IF(lwp) WRITE(numout,*) ' triple point of temperature rtt = ', rtt , ' K' 161 IF(lwp) WRITE(numout,*) ' freezing point of water rt0 = ', rt0 , ' K' 162 IF(lwp) WRITE(numout,*) ' melting point of snow rt0_snow = ', rt0_snow, ' K' 163 IF(lwp) WRITE(numout,*) ' melting point of ice rt0_ice = ', rt0_ice , ' K' 164 165 r1_rau0 = 1._wp / rau0 166 r1_rcp = 1._wp / rcp 167 r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 168 IF(lwp) WRITE(numout,*) 169 IF(lwp) WRITE(numout,*) ' volumic mass of pure water rauw = ', rauw , ' kg/m^3' 170 IF(lwp) WRITE(numout,*) ' volumic mass of reference rau0 = ', rau0 , ' kg/m^3' 171 IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg' 172 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 173 IF(lwp) WRITE(numout,*) ' 1. / ( rau0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp 174 175 176 #if defined key_lim3 || defined key_cice 177 xlsn = lfus * rhosn ! volumetric latent heat fusion of snow [J/m3] 178 #else 179 cpic = rcpic / rhoic ! specific heat for ice [J/Kg/K] 180 lfus = xlsn / rhosn ! latent heat of fusion of fresh ice 181 #endif 182 183 IF(lwp) THEN 133 WRITE(numout,*) 134 WRITE(numout,*) ' Constants' 135 WRITE(numout,*) 136 WRITE(numout,*) ' mathematical constant rpi = ', rpi 137 WRITE(numout,*) ' day rday = ', rday, ' s' 138 WRITE(numout,*) ' sideral year rsiyea = ', rsiyea, ' s' 139 WRITE(numout,*) ' sideral day rsiday = ', rsiday, ' s' 140 WRITE(numout,*) ' omega omega = ', omega, ' s-1' 141 WRITE(numout,*) 142 WRITE(numout,*) ' nb of months per year raamo = ', raamo, ' months' 143 WRITE(numout,*) ' nb of hours per day rjjhh = ', rjjhh, ' hours' 144 WRITE(numout,*) ' nb of minutes per hour rhhmm = ', rhhmm, ' mn' 145 WRITE(numout,*) ' nb of seconds per minute rmmss = ', rmmss, ' s' 146 WRITE(numout,*) 147 WRITE(numout,*) ' earth radius ra = ', ra, ' m' 148 WRITE(numout,*) ' gravity grav = ', grav , ' m/s^2' 149 WRITE(numout,*) 150 WRITE(numout,*) ' triple point of temperature rtt = ', rtt , ' K' 151 WRITE(numout,*) ' freezing point of water rt0 = ', rt0 , ' K' 152 WRITE(numout,*) ' melting point of snow rt0_snow = ', rt0_snow, ' K' 153 WRITE(numout,*) ' melting point of ice rt0_ice = ', rt0_ice , ' K' 154 WRITE(numout,*) 155 WRITE(numout,*) ' ocean reference volumic mass rau0 = ', rau0 , ' kg/m^3' 156 WRITE(numout,*) ' ocean reference specific volume rau0r = ', rau0r, ' m^3/Kg' 157 WRITE(numout,*) ' ocean specific heat rcp = ', rcp 158 WRITE(numout,*) ' 1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr 184 159 WRITE(numout,*) 185 160 WRITE(numout,*) ' thermal conductivity of the snow = ', rcdsn , ' J/s/m/K' 186 161 WRITE(numout,*) ' thermal conductivity of the ice = ', rcdic , ' J/s/m/K' 162 #if defined key_lim3 187 163 WRITE(numout,*) ' fresh ice specific heat = ', cpic , ' J/kg/K' 188 164 WRITE(numout,*) ' latent heat of fusion of fresh ice / snow = ', lfus , ' J/kg' 189 #if defined key_lim3 || defined key_cice190 165 WRITE(numout,*) ' latent heat of subl. of fresh ice / snow = ', lsub , ' J/kg' 166 #elif defined key_cice 167 WRITE(numout,*) ' latent heat of fusion of fresh ice / snow = ', lfus , ' J/kg' 191 168 #else 192 169 WRITE(numout,*) ' density times specific heat for snow = ', rcpsn , ' J/m^3/K' 193 170 WRITE(numout,*) ' density times specific heat for ice = ', rcpic , ' J/m^3/K' 194 171 WRITE(numout,*) ' volumetric latent heat fusion of sea ice = ', xlic , ' J/m' 172 WRITE(numout,*) ' volumetric latent heat fusion of snow = ', xlsn , ' J/m' 195 173 WRITE(numout,*) ' latent heat of sublimation of snow = ', xsn , ' J/kg' 196 174 #endif 197 WRITE(numout,*) ' volumetric latent heat fusion of snow = ', xlsn , ' J/m^3'198 175 WRITE(numout,*) ' density of sea ice = ', rhoic , ' kg/m^3' 199 176 WRITE(numout,*) ' density of snow = ', rhosn , ' kg/m^3'
Note: See TracChangeset
for help on using the changeset viewer.
