Changeset 3453
- Timestamp:
- 2012-08-14T15:38:51+02:00 (12 years ago)
- Location:
- branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM
- Files:
-
- 4 edited
-
CONFIG/AMM12/EXP00/namelist (modified) (3 diffs)
-
NEMO/OPA_SRC/DOM/domwri.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/par_AMM_12km.h90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12/EXP00/namelist
r3309 r3453 65 65 &namzgr_sco ! s-coordinate or hybrid z-s-coordinate 66 66 !----------------------------------------------------------------------- 67 rn_sbot_min = 10. ! minimum depth of s-bottom surface (>0) (m) 68 rn_sbot_max = 7000. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 67 ln_s_sh94 = .false. ! Song & Haidvogel 1994 hybrid S-sigma (T)| 68 ln_s_sf12 = .true. ! Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 69 ln_sigcrit = .true. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 70 rn_sbot_min = 10.0 ! minimum depth of s-bottom surface (>0) (m) 71 rn_sbot_max = 7000.0 ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 72 rn_hc = 50.0 ! critical depth for transition to stretched coordinates 73 rn_rmax = 0.3 ! maximum cut-off r-value allowed (0<r_max<1) 74 ! SH94 stretching coefficients 69 75 rn_theta = 6.0 ! surface control parameter (0<=theta<=20) 70 rn_thetb = 1.00 ! bottom control parameter (0<=thetb<= 1) 71 rn_rmax = 0.30 ! maximum cut-off r-value allowed (0<r_max<1) 72 ln_s_sigma = .true. ! hybrid s-sigma coordinates 73 rn_bb = 0.8 ! stretching with s-sigma 74 rn_hc = 150.0 ! critical depth with s-sigma 76 rn_thetb = 1.0 ! bottom control parameter (0<=thetb<= 1) 77 rn_bb = 0.8 ! stretching with SH94 s-sigma 78 ! SF12 stretching coefficient 79 rn_alpha = 4.4 ! stretching with SH94 s-sigma 80 rn_efold = 0.0 ! efold length scale for transition to stretched coord 81 rn_zs = 1.0 ! depth of surface grid box 82 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 83 rn_zb_a = 0.024 ! bathymetry scaling factor for calculating Zb 84 rn_zb_b = -0.2 ! offset for calculating Zb 75 85 / 76 86 !----------------------------------------------------------------------- … … 79 89 nn_bathy = 1 ! compute (=0) or read (=1) the bathymetry file 80 90 nn_closea = 0 ! remove (=0) or keep (=1) closed seas and lakes (ORCA) 81 nn_msh = 0! create (=1) a mesh file or not (=0)91 nn_msh = 1 ! create (=1) a mesh file or not (=0) 82 92 rn_hmin = -3. ! min depth of the ocean (>0) or min number of ocean level (<0) 83 93 rn_e3zps_min= 20. ! partial step thickness is set larger than the minimum of … … 419 429 bn_v3d = 'amm12_bdyV_u3d' , 24 , 'vomecrty' , .true. , .false. , 'daily' , '' , '' 420 430 bn_tem = 'amm12_bdyT_tra' , 24 , 'votemper' , .true. , .false. , 'daily' , '' , '' 421 bn_ sal= 'amm12_bdyT_tra' , 24 , 'vosaline' , .true. , .false. , 'daily' , '' , ''431 bn_tem = 'amm12_bdyT_tra' , 24 , 'vosaline' , .true. , .false. , 'daily' , '' , '' 422 432 cn_dir = 'bdydta/' 423 433 ln_full_vel = .false. -
branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r3294 r3453 188 188 CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 189 189 ! 190 CALL iom_rstput( 0, 0, inum4, 'gdept _0' , gdept_0) ! ! stretched system191 CALL iom_rstput( 0, 0, inum4, 'gdepw _0' , gdepw_0 )190 CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept ) ! ! stretched system 191 CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw ) 192 192 ENDIF 193 193 -
branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3421 r3453 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 using the Siddorn and Furner 2012? 1730 !! mixed S/sigma/Z coordinate 1731 !! 1732 !! 1733 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1734 !!---------------------------------------------------------------------- 1735 1736 INTEGER :: ji, jj, jk ! dummy loop argument 1737 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1738 1739 DO jk = 1, jpk 1740 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1741 gsigt(jk) = -fssig( REAL(jk,wp) ) 1742 END DO 1743 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1744 ! 1745 ! Coefficients for vertical scale factors at w-, t- levels 1746 !!gm bug : define it from analytical function, not like juste bellow.... 1747 !!gm or betteroffer the 2 possibilities.... 1748 DO jk = 1, jpkm1 1749 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1750 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1751 END DO 1752 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1753 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1754 ! 1755 ! Coefficients for vertical depth as the sum of e3w scale factors 1756 gsi3w(1) = 0.5_wp * esigw(1) 1757 DO jk = 2, jpk 1758 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1759 END DO 1760 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1761 DO jk = 1, jpk 1762 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1763 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1764 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 1765 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 1766 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 1767 END DO 1768 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1769 DO jj = 1, jpj 1770 DO ji = 1, jpi 1771 DO jk = 1, jpk 1772 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1773 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1774 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1775 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1776 ! 1777 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1778 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1779 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1780 END DO 1781 END DO 1782 END DO 1783 END SUBROUTINE s_tanh 1784 1785 FUNCTION fssig( pk ) RESULT( pf ) 1786 !!---------------------------------------------------------------------- 1787 !! *** ROUTINE fssig *** 1788 !! 1789 !! ** Purpose : provide the analytical function in s-coordinate 1790 !! 1791 !! ** Method : the function provide the non-dimensional position of 1792 !! T and W (i.e. between 0 and 1) 1793 !! T-points at integer values (between 1 and jpk) 1794 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1795 !!---------------------------------------------------------------------- 1796 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 1797 REAL(wp) :: pf ! sigma value 1798 !!---------------------------------------------------------------------- 1799 ! 1800 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 1801 & - TANH( rn_thetb * rn_theta ) ) & 1802 & * ( COSH( rn_theta ) & 1803 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1804 & / ( 2._wp * SINH( rn_theta ) ) 1805 ! 1806 END FUNCTION fssig 1807 1808 1809 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 1810 !!---------------------------------------------------------------------- 1811 !! *** ROUTINE fssig1 *** 1812 !! 1813 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 1814 !! 1815 !! ** Method : the function provides the non-dimensional position of 1816 !! T and W (i.e. between 0 and 1) 1817 !! T-points at integer values (between 1 and jpk) 1818 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1819 !!---------------------------------------------------------------------- 1820 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 1821 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 1822 REAL(wp) :: pf1 ! sigma value 1823 !!---------------------------------------------------------------------- 1824 ! 1825 IF ( rn_theta == 0 ) then ! uniform sigma 1826 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 1827 ELSE ! stretched sigma 1828 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 1829 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1830 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1831 ENDIF 1832 ! 1833 END FUNCTION fssig1 1834 1835 1836 FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam ) 1837 !!---------------------------------------------------------------------- 1838 !! *** ROUTINE fgamma *** 1839 !! 1840 !! ** Purpose : provide analytical function for the s-coordinate 1841 !! 1842 !! ** Method : the function provides the non-dimensional position of 1843 !! T and W (i.e. between 0 and 1) 1844 !! T-points at integer values (between 1 and jpk) 1845 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1846 !! 1847 !! This method allows the maintenance of fixed surface and or 1848 !! bottom cell resolutions (cf. geopotential coordinates) 1849 !! within an analytically derived stretched S-coordinate framework. 1850 !! 1851 !! Reference : Siddorn and Furner, in prep 1852 !!---------------------------------------------------------------------- 1853 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 1854 REAL(wp) :: gam(jpk) ! stretched coordinate 1855 REAL(wp), INTENT(in ) :: Zbb ! Bottom box depth 1856 REAL(wp), INTENT(in ) :: Zss ! surface box depth 1857 REAL(wp), INTENT(in ) :: F ! Smoothing parameter 1858 REAL(wp) :: a1,a2,a3 ! local variables 1859 REAL(wp) :: n1,n2 ! local variables 1860 REAL(wp) :: A,B,X ! local variables 1861 integer :: jk 1862 !!---------------------------------------------------------------------- 1863 ! 1864 1865 n1 = 1./(jpk-1.) 1866 n2 = 1. - n1 1867 1868 a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp) 1869 a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp) 1870 a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1) 1871 1872 A = Zbb - a3*(Zss-a1)-a2 1873 A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) ) 1874 B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1) 1875 X = 1.0_wp-A/2.0_wp-B 1876 1877 DO jk = 1, jpk 1878 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) ) 1879 gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F) 1880 ENDDO 1881 1882 ! 1883 END FUNCTION fgamma 1600 1884 1601 1885 !!====================================================================== -
branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90
r3294 r3453 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.
