Changeset 3600 for branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC
- Timestamp:
- 2012-11-19T14:59:22+01:00 (12 years ago)
- Location:
- branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r3421 r3600 174 174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at V--F 175 175 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing T--U points (m) 176 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio 176 177 177 178 !!---------------------------------------------------------------------- … … 294 295 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 295 296 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 296 & hift (jpi,jpj) , hifu (jpi,jpj) , STAT=ierr(8) )297 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 297 298 298 299 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & -
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r3421 r3600 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) 38 39 39 40 IMPLICIT NONE … … 84 85 CALL dom_zgr ! Vertical mesh and bathymetry 85 86 CALL dom_msk ! Masks 87 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 86 88 IF( lk_vvl ) CALL dom_vvl ! Vertical variable mesh 87 89 ! … … 322 324 END SUBROUTINE dom_ctl 323 325 326 SUBROUTINE dom_stiff 327 !!---------------------------------------------------------------------- 328 !! *** ROUTINE dom_stiff *** 329 !! 330 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency 331 !! 332 !! ** Method : Compute Haney (1991) hydrostatic condition ratio 333 !! Save the maximum in the vertical direction 334 !! (this number is only relevant in s-coordinates) 335 !! 336 !! Haney, R. L., 1991: On the pressure gradient force 337 !! over steep topography in sigma coordinate ocean models. 338 !! J. Phys. Oceanogr., 21, 610???619. 339 !!---------------------------------------------------------------------- 340 INTEGER :: ji, jj, jk 341 REAL(wp) :: zrxmax 342 REAL(wp), DIMENSION(4) :: zr1 343 !!---------------------------------------------------------------------- 344 rx1(:,:) = 0.e0 345 zrxmax = 0.e0 346 zr1(:) = 0.e0 347 348 DO ji = 2, jpim1 349 DO jj = 2, jpjm1 350 DO jk = 1, jpkm1 351 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji-1,jj ,jk ) & 352 & +gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1)) & 353 & /(gdepw(ji ,jj ,jk )+gdepw(ji-1,jj ,jk ) & 354 & -gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1) + rsmall) ) 355 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw(ji+1,jj ,jk )-gdepw(ji ,jj ,jk ) & 356 & +gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1)) & 357 & /(gdepw(ji+1,jj ,jk )+gdepw(ji ,jj ,jk ) & 358 & -gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) ) 359 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw(ji ,jj+1,jk )-gdepw(ji ,jj ,jk ) & 360 & +gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1)) & 361 & /(gdepw(ji ,jj+1,jk )+gdepw(ji ,jj ,jk ) & 362 & -gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) ) 363 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji ,jj-1,jk ) & 364 & +gdepw(ji ,jj ,jk+1)-gdepw(ji ,jj-1,jk+1)) & 365 & /(gdepw(ji ,jj ,jk )+gdepw(ji ,jj-1,jk ) & 366 & -gdepw(ji, jj ,jk+1)-gdepw(ji ,jj-1,jk+1) + rsmall) ) 367 zrxmax = MAXVAL(zr1(1:4)) 368 rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) 369 END DO 370 END DO 371 END DO 372 373 CALL lbc_lnk( rx1, 'T', 1. ) 374 375 zrxmax = MAXVAL(rx1) 376 377 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain 378 379 IF(lwp) THEN 380 WRITE(numout,*) 381 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 382 WRITE(numout,*) '~~~~~~~~~' 383 ENDIF 384 385 END SUBROUTINE dom_stiff 386 387 388 324 389 !!====================================================================== 325 390 END MODULE domain -
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r3294 r3600 172 172 173 173 IF( ln_sco ) THEN ! s-coordinate 174 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) ! ! depth175 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 174 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) 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 ) … … 187 187 CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 188 188 CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 189 ! 190 CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 ) ! ! stretched system 191 CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 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 ) 192 193 ENDIF 193 194 -
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3421 r3600 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 adn Furner stretching functio 17 18 !!---------------------------------------------------------------------- 18 19 … … 28 29 !! zgr_sco : s-coordinate 29 30 !! fssig : sigma coordinate non-dimensional function 31 !! fgamma : Siddorn and Furner stretching function 30 32 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) 31 33 !!--------------------------------------------------------------------- … … 47 49 48 50 ! !!* Namelist namzgr_sco * 51 LOGICAL :: ln_s_sh94 = .false. ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 52 LOGICAL :: ln_s_sf12 = .true. ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 53 LOGICAL :: ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 54 ! 49 55 REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) 50 56 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 coordinates 59 ! Song and Haidvogel 1994 stretching parameters 51 60 REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20) 52 61 REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1) 53 REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1) 54 LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 55 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter for song and haidvogel stretching 62 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter 56 63 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 57 REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates 64 ! Siddorn and Furner stretching parameters 65 REAL(wp) :: rn_alpha = 4.4_wp ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 66 REAL(wp) :: rn_efold = 0.0_wp ! efold length scale for transition to stretched coord 67 REAL(wp) :: rn_zs = 1.0_wp ! depth of surface grid box 68 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 69 REAL(wp) :: rn_zb_a = 0.024_wp ! bathymetry scaling factor for calculating Zb 70 REAL(wp) :: rn_zb_b = -0.2_wp ! offset for calculating Zb 58 71 59 72 !! * Substitutions … … 1034 1047 END SUBROUTINE zgr_zps 1035 1048 1036 1037 FUNCTION fssig( pk ) RESULT( pf )1038 !!----------------------------------------------------------------------1039 !! *** ROUTINE eos_init ***1040 !!1041 !! ** Purpose : provide the analytical function in s-coordinate1042 !!1043 !! ** Method : the function provide the non-dimensional position of1044 !! T and W (i.e. between 0 and 1)1045 !! T-points at integer values (between 1 and jpk)1046 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)1047 !!----------------------------------------------------------------------1048 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate1049 REAL(wp) :: pf ! sigma value1050 !!----------------------------------------------------------------------1051 !1052 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &1053 & - TANH( rn_thetb * rn_theta ) ) &1054 & * ( COSH( rn_theta ) &1055 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &1056 & / ( 2._wp * SINH( rn_theta ) )1057 !1058 END FUNCTION fssig1059 1060 1061 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )1062 !!----------------------------------------------------------------------1063 !! *** ROUTINE eos_init ***1064 !!1065 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate1066 !!1067 !! ** Method : the function provides the non-dimensional position of1068 !! T and W (i.e. between 0 and 1)1069 !! T-points at integer values (between 1 and jpk)1070 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)1071 !!----------------------------------------------------------------------1072 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate1073 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient1074 REAL(wp) :: pf1 ! sigma value1075 !!----------------------------------------------------------------------1076 !1077 IF ( rn_theta == 0 ) then ! uniform sigma1078 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )1079 ELSE ! stretched sigma1080 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &1081 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) &1082 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )1083 ENDIF1084 !1085 END FUNCTION fssig11086 1087 1088 1049 SUBROUTINE zgr_sco 1089 1050 !!---------------------------------------------------------------------- … … 1110 1071 !! esigt(k) = fsdsig(k ) 1111 1072 !! esigw(k) = fsdsig(k-0.5) 1112 !! Th is routine is given as an example, it mustbe modified1113 !! following the user s desiderata. nevertheless, the output as1073 !! Three options for stretching are give, and they can be modified 1074 !! following the users requirements. Nevertheless, the output as 1114 1075 !! well as the way to compute the model levels and scale factors 1115 !! must be respected in order to insure second order a !!uracy1076 !! must be respected in order to insure second order accuracy 1116 1077 !! schemes. 1117 1078 !! 1118 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1079 !! The three methods for stretching available are: 1080 !! 1081 !! s_sh94 (Song and Haidvogel 1994) 1082 !! a sinh/tanh function that allows sigma and stretched sigma 1083 !! 1084 !! s_sf12 (Siddorn and Furner 2012?) 1085 !! allows the maintenance of fixed surface and or 1086 !! bottom cell resolutions (cf. geopotential coordinates) 1087 !! within an analytically derived stretched S-coordinate framework. 1088 !! 1089 !! s_tanh (Madec et al 1996) 1090 !! a cosh/tanh function that gives stretched coordinates 1091 !! 1119 1092 !!---------------------------------------------------------------------- 1120 1093 ! 1121 1094 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1122 1095 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1123 REAL(wp) :: z coeft, zcoefw, zrmax, ztaper ! temporary scalars1096 REAL(wp) :: zrmax, ztaper ! temporary scalars 1124 1097 ! 1125 1098 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1126 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1127 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1128 1129 NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 1130 !!---------------------------------------------------------------------- 1099 1100 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1101 rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1102 !!---------------------------------------------------------------------- 1131 1103 ! 1132 1104 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1133 1105 ! 1134 1106 CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1135 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 )1136 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )1137 1107 ! 1138 1108 REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters … … 1144 1114 WRITE(numout,*) '~~~~~~~~~~~' 1145 1115 WRITE(numout,*) ' Namelist namzgr_sco' 1146 WRITE(numout,*) ' sigma-stretching coeffs ' 1147 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max 1148 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min 1149 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta 1150 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb 1151 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax 1152 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1153 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb 1154 WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc 1155 ENDIF 1156 1157 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1158 esigt3 = 0._wp ; esigw3 = 0._wp 1159 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1160 esigwu3 = 0._wp ; esigwv3 = 0._wp 1116 WRITE(numout,*) ' stretching coeffs ' 1117 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max 1118 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min 1119 WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc 1120 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax 1121 WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 1122 WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients' 1123 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta 1124 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb 1125 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 1126 WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 1127 WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients' 1128 WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha 1129 WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold 1130 WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs 1131 WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a 1132 WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b 1133 WRITE(numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 1134 ENDIF 1161 1135 1162 1136 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate … … 1352 1326 ! non-dimensional "sigma" for model level depth at w- and t-levels 1353 1327 1354 IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths 1355 ! ! below rn_hc, with uniform sigma in shallower waters 1356 DO ji = 1, jpi 1357 DO jj = 1, jpj 1358 1359 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1360 DO jk = 1, jpk 1361 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1362 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1363 END DO 1364 ELSE ! shallow water, uniform sigma 1365 DO jk = 1, jpk 1366 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1367 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1368 END DO 1369 ENDIF 1370 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1371 ! 1372 DO jk = 1, jpkm1 1373 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1374 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1375 END DO 1376 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1377 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1378 ! 1379 ! Coefficients for vertical depth as the sum of e3w scale factors 1380 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1381 DO jk = 2, jpk 1382 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1383 END DO 1384 ! 1385 DO jk = 1, jpk 1386 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1387 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1388 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1389 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1390 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1391 END DO 1392 ! 1393 END DO ! for all jj's 1394 END DO ! for all ji's 1395 1396 DO ji = 1, jpim1 1397 DO jj = 1, jpjm1 1398 DO jk = 1, jpk 1399 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1400 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1401 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1402 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1403 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1404 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1405 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1406 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1407 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1408 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1409 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1410 ! 1411 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1412 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1413 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1414 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1415 ! 1416 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1417 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1418 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1419 END DO 1420 END DO 1421 END DO 1422 1423 CALL lbc_lnk( e3t , 'T', 1._wp ) 1424 CALL lbc_lnk( e3u , 'U', 1._wp ) 1425 CALL lbc_lnk( e3v , 'V', 1._wp ) 1426 CALL lbc_lnk( e3f , 'F', 1._wp ) 1427 CALL lbc_lnk( e3w , 'W', 1._wp ) 1428 CALL lbc_lnk( e3uw, 'U', 1._wp ) 1429 CALL lbc_lnk( e3vw, 'V', 1._wp ) 1430 1431 ! 1432 ELSE ! not ln_s_sigma 1433 ! 1434 DO jk = 1, jpk 1435 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1436 gsigt(jk) = -fssig( REAL(jk,wp) ) 1437 END DO 1438 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1439 ! 1440 ! Coefficients for vertical scale factors at w-, t- levels 1441 !!gm bug : define it from analytical function, not like juste bellow.... 1442 !!gm or betteroffer the 2 possibilities.... 1443 DO jk = 1, jpkm1 1444 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1445 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1446 END DO 1447 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1448 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1449 1450 !!gm original form 1451 !!org DO jk = 1, jpk 1452 !!org esigt(jk)=fsdsig( FLOAT(jk) ) 1453 !!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 1454 !!org END DO 1455 !!gm 1456 ! 1457 ! Coefficients for vertical depth as the sum of e3w scale factors 1458 gsi3w(1) = 0.5_wp * esigw(1) 1459 DO jk = 2, jpk 1460 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1461 END DO 1462 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1463 DO jk = 1, jpk 1464 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1465 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1466 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 1467 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 1468 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 1469 END DO 1470 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1471 DO jj = 1, jpj 1472 DO ji = 1, jpi 1473 DO jk = 1, jpk 1474 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1475 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1476 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1477 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1478 ! 1479 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1480 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1481 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1482 END DO 1483 END DO 1484 END DO 1485 ! 1486 ENDIF ! ln_s_sigma 1487 1488 1328 1329 !======================================================================== 1330 ! Song and Haidvogel 1994 (ln_s_sh94=T) 1331 ! Siddorn and Furner 2012 (ln_sf12=T) 1332 ! or tanh function (both false) 1333 !======================================================================== 1334 IF ( ln_s_sh94 ) THEN 1335 CALL s_sh94() 1336 ELSE IF ( ln_s_sf12 ) THEN 1337 CALL s_sf12() 1338 ELSE 1339 CALL s_tanh() 1340 ENDIF 1341 1342 CALL lbc_lnk( e3t , 'T', 1._wp ) 1343 CALL lbc_lnk( e3u , 'U', 1._wp ) 1344 CALL lbc_lnk( e3v , 'V', 1._wp ) 1345 CALL lbc_lnk( e3f , 'F', 1._wp ) 1346 CALL lbc_lnk( e3w , 'W', 1._wp ) 1347 CALL lbc_lnk( e3uw, 'U', 1._wp ) 1348 CALL lbc_lnk( e3vw, 'V', 1._wp ) 1349 1350 fsdepw(:,:,:) = gdepw (:,:,:) 1351 fsde3w(:,:,:) = gdep3w(:,:,:) 1489 1352 ! 1490 1353 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0 … … 1544 1407 & ' w ', MAXVAL( fse3w (:,:,:) ) 1545 1408 ENDIF 1546 ! 1409 ! END DO 1547 1410 IF(lwp) THEN ! selected vertical profiles 1548 1411 WRITE(numout,*) … … 1574 1437 ENDIF 1575 1438 1576 !!gm bug? no more necessary? if ! defined key_helsinki 1577 DO jk = 1, jpk 1439 !================================================================================ 1440 ! check the coordinate makes sense 1441 !================================================================================ 1442 DO ji = 1, jpi 1578 1443 DO jj = 1, jpj 1579 DO ji = 1, jpi 1580 IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1581 WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1582 CALL ctl_stop( ctmp1 ) 1583 ENDIF 1584 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1585 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1586 CALL ctl_stop( ctmp1 ) 1587 ENDIF 1588 END DO 1589 END DO 1590 END DO 1591 !!gm bug #endif 1444 1445 IF( hbatt(ji,jj) > 0._wp) THEN 1446 DO jk = 1, mbathy(ji,jj) 1447 ! check coordinate is monotonically increasing 1448 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1449 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1450 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1451 WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 1452 WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 1453 CALL ctl_stop( ctmp1 ) 1454 ENDIF 1455 ! and check it has never gone negative 1456 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1457 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1458 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1459 WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1460 WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1461 CALL ctl_stop( ctmp1 ) 1462 ENDIF 1463 ! and check it never exceeds the total depth 1464 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 1465 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1466 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1467 WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1468 CALL ctl_stop( ctmp1 ) 1469 ENDIF 1470 END DO 1471 1472 DO jk = 1, mbathy(ji,jj)-1 1473 ! and check it never exceeds the total depth 1474 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 1475 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1476 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1477 WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1478 CALL ctl_stop( ctmp1 ) 1479 ENDIF 1480 END DO 1481 1482 ENDIF 1483 1484 END DO 1485 END DO 1592 1486 ! 1593 1487 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1488 ! 1489 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1490 ! 1491 END SUBROUTINE zgr_sco 1492 1493 !!====================================================================== 1494 SUBROUTINE s_sh94() 1495 1496 !!---------------------------------------------------------------------- 1497 !! *** ROUTINE s_sh94 *** 1498 !! 1499 !! ** Purpose : stretch the s-coordinate system 1500 !! 1501 !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 1502 !! mixed S/sigma coordinate 1503 !! 1504 !! Reference : Song and Haidvogel 1994. 1505 !!---------------------------------------------------------------------- 1506 ! 1507 INTEGER :: ji, jj, jk ! dummy loop argument 1508 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1509 ! 1510 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1511 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1512 1513 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1514 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1515 1516 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1517 esigt3 = 0._wp ; esigw3 = 0._wp 1518 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1519 esigwu3 = 0._wp ; esigwv3 = 0._wp 1520 1521 DO ji = 1, jpi 1522 DO jj = 1, jpj 1523 1524 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1525 DO jk = 1, jpk 1526 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1527 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1528 END DO 1529 ELSE ! shallow water, uniform sigma 1530 DO jk = 1, jpk 1531 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1532 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1533 END DO 1534 ENDIF 1535 ! 1536 DO jk = 1, jpkm1 1537 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1538 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1539 END DO 1540 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1541 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1542 ! 1543 ! Coefficients for vertical depth as the sum of e3w scale factors 1544 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1545 DO jk = 2, jpk 1546 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1547 END DO 1548 ! 1549 DO jk = 1, jpk 1550 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1551 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1552 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1553 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1554 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1555 END DO 1556 ! 1557 END DO ! for all jj's 1558 END DO ! for all ji's 1559 1560 DO ji = 1, jpim1 1561 DO jj = 1, jpjm1 1562 DO jk = 1, jpk 1563 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1564 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1565 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1566 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1567 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1568 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1569 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1570 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1571 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1572 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1573 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1574 ! 1575 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1576 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1577 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1578 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1579 ! 1580 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1581 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1582 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1583 END DO 1584 END DO 1585 END DO 1586 1594 1587 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1595 1588 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1596 ! 1597 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 1598 ! 1599 END SUBROUTINE zgr_sco 1589 1590 END SUBROUTINE s_sh94 1591 1592 SUBROUTINE s_sf12 1593 1594 !!---------------------------------------------------------------------- 1595 !! *** ROUTINE s_sf12 *** 1596 !! 1597 !! ** Purpose : stretch the s-coordinate system 1598 !! 1599 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 1600 !! mixed S/sigma/Z coordinate 1601 !! 1602 !! This method allows the maintenance of fixed surface and or 1603 !! bottom cell resolutions (cf. geopotential coordinates) 1604 !! within an analytically derived stretched S-coordinate framework. 1605 !! 1606 !! 1607 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 1608 !!---------------------------------------------------------------------- 1609 ! 1610 INTEGER :: ji, jj, jk ! dummy loop argument 1611 REAL(wp) :: fsmth ! smoothing around critical depth 1612 REAL(wp) :: zss, zbb ! Surface and bottom cell thickness in sigma space 1613 ! 1614 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1615 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1616 1617 ! 1618 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1619 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1620 1621 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1622 esigt3 = 0._wp ; esigw3 = 0._wp 1623 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1624 esigwu3 = 0._wp ; esigwv3 = 0._wp 1625 1626 DO ji = 1, jpi 1627 DO jj = 1, jpj 1628 1629 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 1630 1631 zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 1632 ! could be changed by users but care must be taken to do so carefully 1633 zbb = 1.0_wp-(zbb/hbatt(ji,jj)) 1634 1635 zss = rn_zs / hbatt(ji,jj) 1636 1637 IF (rn_efold /= 0.0_wp) THEN 1638 fsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 1639 ELSE 1640 fsmth = 1.0_wp 1641 ENDIF 1642 1643 DO jk = 1, jpk 1644 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1645 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 1646 ENDDO 1647 gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth ) 1648 gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth ) 1649 1650 ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 1651 1652 DO jk = 1, jpk 1653 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1654 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 1655 END DO 1656 1657 ELSE ! shallow water, z coordinates 1658 1659 DO jk = 1, jpk 1660 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 1661 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 1662 END DO 1663 1664 ENDIF 1665 1666 DO jk = 1, jpkm1 1667 esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1668 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1669 END DO 1670 esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 )) 1671 esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk)) 1672 1673 ! Coefficients for vertical depth as the sum of e3w scale factors 1674 gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 1675 DO jk = 2, jpk 1676 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1677 END DO 1678 1679 DO jk = 1, jpk 1680 gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigt3(ji,jj,jk) 1681 gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigw3(ji,jj,jk) 1682 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsi3w3(ji,jj,jk) 1683 END DO 1684 1685 ENDDO ! for all jj's 1686 ENDDO ! for all ji's 1687 1688 DO ji=1,jpi 1689 DO jj=1,jpj 1690 1691 DO jk = 1, jpk 1692 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 1693 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1694 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 1695 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1696 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) + & 1697 hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 1698 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1699 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 1700 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1701 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 1702 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1703 1704 e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*esigt3(ji,jj,jk) 1705 e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*esigtu3(ji,jj,jk) 1706 e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*esigtv3(ji,jj,jk) 1707 e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*esigtf3(ji,jj,jk) 1708 ! 1709 e3w(ji,jj,jk)=hbatt(ji,jj)*esigw3(ji,jj,jk) 1710 e3uw(ji,jj,jk)=hbatu(ji,jj)*esigwu3(ji,jj,jk) 1711 e3vw(ji,jj,jk)=hbatv(ji,jj)*esigwv3(ji,jj,jk) 1712 END DO 1713 1714 ENDDO 1715 ENDDO 1716 1717 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1718 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1719 1720 END SUBROUTINE s_sf12 1721 1722 SUBROUTINE s_tanh() 1723 1724 !!---------------------------------------------------------------------- 1725 !! *** ROUTINE s_tanh*** 1726 !! 1727 !! ** Purpose : stretch the s-coordinate system 1728 !! 1729 !! ** Method : s-coordinate stretch 1730 !! 1731 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1732 !!---------------------------------------------------------------------- 1733 1734 INTEGER :: ji, jj, jk ! dummy loop argument 1735 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1736 1737 DO jk = 1, jpk 1738 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1739 gsigt(jk) = -fssig( REAL(jk,wp) ) 1740 END DO 1741 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1742 ! 1743 ! Coefficients for vertical scale factors at w-, t- levels 1744 !!gm bug : define it from analytical function, not like juste bellow.... 1745 !!gm or betteroffer the 2 possibilities.... 1746 DO jk = 1, jpkm1 1747 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1748 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1749 END DO 1750 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1751 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1752 ! 1753 ! Coefficients for vertical depth as the sum of e3w scale factors 1754 gsi3w(1) = 0.5_wp * esigw(1) 1755 DO jk = 2, jpk 1756 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1757 END DO 1758 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1759 DO jk = 1, jpk 1760 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1761 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1762 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 1763 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 1764 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 1765 END DO 1766 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1767 DO jj = 1, jpj 1768 DO ji = 1, jpi 1769 DO jk = 1, jpk 1770 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1771 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1772 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1773 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1774 ! 1775 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1776 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1777 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1778 END DO 1779 END DO 1780 END DO 1781 END SUBROUTINE s_tanh 1782 1783 FUNCTION fssig( pk ) RESULT( pf ) 1784 !!---------------------------------------------------------------------- 1785 !! *** ROUTINE fssig *** 1786 !! 1787 !! ** Purpose : provide the analytical function in s-coordinate 1788 !! 1789 !! ** Method : the function provide the non-dimensional position of 1790 !! T and W (i.e. between 0 and 1) 1791 !! T-points at integer values (between 1 and jpk) 1792 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1793 !!---------------------------------------------------------------------- 1794 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 1795 REAL(wp) :: pf ! sigma value 1796 !!---------------------------------------------------------------------- 1797 ! 1798 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 1799 & - TANH( rn_thetb * rn_theta ) ) & 1800 & * ( COSH( rn_theta ) & 1801 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1802 & / ( 2._wp * SINH( rn_theta ) ) 1803 ! 1804 END FUNCTION fssig 1805 1806 1807 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 1808 !!---------------------------------------------------------------------- 1809 !! *** ROUTINE fssig1 *** 1810 !! 1811 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 1812 !! 1813 !! ** Method : the function provides the non-dimensional position of 1814 !! T and W (i.e. between 0 and 1) 1815 !! T-points at integer values (between 1 and jpk) 1816 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1817 !!---------------------------------------------------------------------- 1818 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 1819 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 1820 REAL(wp) :: pf1 ! sigma value 1821 !!---------------------------------------------------------------------- 1822 ! 1823 IF ( rn_theta == 0 ) then ! uniform sigma 1824 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 1825 ELSE ! stretched sigma 1826 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 1827 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1828 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1829 ENDIF 1830 ! 1831 END FUNCTION fssig1 1832 1833 1834 FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam ) 1835 !!---------------------------------------------------------------------- 1836 !! *** ROUTINE fgamma *** 1837 !! 1838 !! ** Purpose : provide analytical function for the s-coordinate 1839 !! 1840 !! ** Method : the function provides the non-dimensional position of 1841 !! T and W (i.e. between 0 and 1) 1842 !! T-points at integer values (between 1 and jpk) 1843 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1844 !! 1845 !! This method allows the maintenance of fixed surface and or 1846 !! bottom cell resolutions (cf. geopotential coordinates) 1847 !! within an analytically derived stretched S-coordinate framework. 1848 !! 1849 !! Reference : Siddorn and Furner, in prep 1850 !!---------------------------------------------------------------------- 1851 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 1852 REAL(wp) :: gam(jpk) ! stretched coordinate 1853 REAL(wp), INTENT(in ) :: Zbb ! Bottom box depth 1854 REAL(wp), INTENT(in ) :: Zss ! surface box depth 1855 REAL(wp), INTENT(in ) :: F ! Smoothing parameter 1856 REAL(wp) :: a1,a2,a3 ! local variables 1857 REAL(wp) :: n1,n2 ! local variables 1858 REAL(wp) :: A,B,X ! local variables 1859 integer :: jk 1860 !!---------------------------------------------------------------------- 1861 ! 1862 1863 n1 = 1./(jpk-1.) 1864 n2 = 1. - n1 1865 1866 a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp) 1867 a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp) 1868 a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1) 1869 1870 A = Zbb - a3*(Zss-a1)-a2 1871 A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) ) 1872 B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1) 1873 X = 1.0_wp-A/2.0_wp-B 1874 1875 DO jk = 1, jpk 1876 gam(jk) = A*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+B*pk1(jk)**3.0_wp + X*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 1877 gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F) 1878 ENDDO 1879 1880 ! 1881 END FUNCTION fgamma 1600 1882 1601 1883 !!====================================================================== -
branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90
r3294 r3600 19 19 jpidta = 198, & !: first horizontal dimension > or = to jpi 20 20 jpjdta = 224, & !: second > or = to jpj 21 jpkdta = 33, & !: number of levels > or = to jpk21 jpkdta = 51, & !: number of levels > or = to jpk 22 22 ! total domain matrix size 23 23 jpiglo = jpidta, & !: first dimension of global domain --> i
Note: See TracChangeset
for help on using the changeset viewer.