Newer
Older
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
rpsi3p = 1._wp
rpsi2 = 0.5_wp
!
SELECT CASE ( nn_stab_func )
CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions
CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions
CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified)
END SELECT
!
CASE( 1 ) ! k-eps
!
IF(lwp) WRITE(numout,*) ' ==>> k-eps closure chosen'
IF(lwp) WRITE(numout,*)
rpp = 3._wp
rmm = 1.5_wp
rnn = -1._wp
rsc_tke = 1._wp
rsc_psi = 1.2_wp ! Schmidt number for psi
rpsi1 = 1.44_wp
rpsi3p = 1._wp
rpsi2 = 1.92_wp
!
SELECT CASE ( nn_stab_func )
CASE( 0, 1 ) ; rpsi3m = -0.52_wp ! G88 or KC stability functions
CASE( 2 ) ; rpsi3m = -0.629_wp ! Canuto A stability functions
CASE( 3 ) ; rpsi3m = -0.566 ! Canuto B stability functions
END SELECT
!
CASE( 2 ) ! k-omega
!
IF(lwp) WRITE(numout,*) ' ==>> k-omega closure chosen'
IF(lwp) WRITE(numout,*)
rpp = -1._wp
rmm = 0.5_wp
rnn = -1._wp
rsc_tke = 2._wp
rsc_psi = 2._wp
rpsi1 = 0.555_wp
rpsi3p = 1._wp
rpsi2 = 0.833_wp
!
SELECT CASE ( nn_stab_func )
CASE( 0, 1 ) ; rpsi3m = -0.58_wp ! G88 or KC stability functions
CASE( 2 ) ; rpsi3m = -0.64_wp ! Canuto A stability functions
CASE( 3 ) ; rpsi3m = -0.64_wp ! Canuto B stability functions caution : constant not identified)
END SELECT
!
CASE( 3 ) ! generic
!
IF(lwp) WRITE(numout,*) ' ==>> generic closure chosen'
IF(lwp) WRITE(numout,*)
rpp = 2._wp
rmm = 1._wp
rnn = -0.67_wp
rsc_tke = 0.8_wp
rsc_psi = 1.07_wp
rpsi1 = 1._wp
rpsi3p = 1._wp
rpsi2 = 1.22_wp
!
SELECT CASE ( nn_stab_func )
CASE( 0, 1 ) ; rpsi3m = 0.1_wp ! G88 or KC stability functions
CASE( 2 ) ; rpsi3m = 0.05_wp ! Canuto A stability functions
CASE( 3 ) ; rpsi3m = 0.05_wp ! Canuto B stability functions caution : constant not identified)
END SELECT
!
END SELECT
!
SELECT CASE ( nn_stab_func ) !* set the parameters of the stability functions
!
CASE ( 0 ) ! Galperin stability functions
!
IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Galperin'
rc2 = 0._wp
rc3 = 0._wp
rc_diff = 1._wp
rc0 = 0.5544_wp
rcm_sf = 0.9884_wp
rghmin = -0.28_wp
rgh0 = 0.0233_wp
rghcri = 0.02_wp
!
CASE ( 1 ) ! Kantha-Clayson stability functions
!
IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Kantha-Clayson'
rc2 = 0.7_wp
rc3 = 0.2_wp
rc_diff = 1._wp
rc0 = 0.5544_wp
rcm_sf = 0.9884_wp
rghmin = -0.28_wp
rgh0 = 0.0233_wp
rghcri = 0.02_wp
!
CASE ( 2 ) ! Canuto A stability functions
!
IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto A'
rs0 = 1.5_wp * rl1 * rl5*rl5
rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
rs4 = 2._wp * rl5
rs5 = 2._wp * rl4
rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2) &
& + 0.75_wp * rl1 * ( rl6 - rl7 )
rd0 = 3._wp * rl5*rl5
rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
rc0 = 0.5268_wp
rf6 = 8._wp / (rc0**6._wp)
rc_diff = SQRT(2._wp) / (rc0**3._wp)
rcm_sf = 0.7310_wp
rghmin = -0.28_wp
rgh0 = 0.0329_wp
rghcri = 0.03_wp
!
CASE ( 3 ) ! Canuto B stability functions
!
IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto B'
rs0 = 1.5_wp * rm1 * rm5*rm5
rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
rs4 = 2._wp * rm5
rs5 = 2._wp * rm4
rs6 = (2._wp/3._wp) * rm5 * (3._wp*rm3*rm3-rm2*rm2) - 0.5_wp * rm5*rm1*(3._wp*rm3-rm2) + 0.75_wp * rm1*(rm6-rm7)
rd0 = 3._wp * rm5*rm5
rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
rc0 = 0.5268_wp !! rc0 = 0.5540_wp (Warner ...) to verify !
rf6 = 8._wp / ( rc0**6._wp )
rc_diff = SQRT(2._wp)/(rc0**3.)
rcm_sf = 0.7470_wp
rghmin = -0.28_wp
rgh0 = 0.0444_wp
rghcri = 0.0414_wp
!
END SELECT
! !* Set Schmidt number for psi diffusion in the wave breaking case
! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
IF( ln_sigpsi ) THEN
ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
! Verification: retrieve Burchard (2001) results by uncomenting the line below:
! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn
rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.))
ELSE
rsc_psi0 = rsc_psi
ENDIF
! !* Shear free turbulence parameters
!
ra_sf = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) &
& - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
IF ( rn_crban==0._wp ) THEN
rl_sf = vkarmn
ELSE
rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke &
& + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) &
& *SQRT(rsc_tke*(rsc_tke &
& + 24._wp*rsc_psi0*rpsi2)) ) &
& /(12._wp*rnn**2.) )
ENDIF
!
IF(lwp) THEN !* Control print
WRITE(numout,*)
WRITE(numout,*) ' Limit values :'
WRITE(numout,*) ' Parameter m = ', rmm
WRITE(numout,*) ' Parameter n = ', rnn
WRITE(numout,*) ' Parameter p = ', rpp
WRITE(numout,*) ' rpsi1 = ', rpsi1
WRITE(numout,*) ' rpsi2 = ', rpsi2
WRITE(numout,*) ' rpsi3m = ', rpsi3m
WRITE(numout,*) ' rpsi3p = ', rpsi3p
WRITE(numout,*) ' rsc_tke = ', rsc_tke
WRITE(numout,*) ' rsc_psi = ', rsc_psi
WRITE(numout,*) ' rsc_psi0 = ', rsc_psi0
WRITE(numout,*) ' rc0 = ', rc0
WRITE(numout,*)
WRITE(numout,*) ' Shear free turbulence parameters:'
WRITE(numout,*) ' rcm_sf = ', rcm_sf
WRITE(numout,*) ' ra_sf = ', ra_sf
WRITE(numout,*) ' rl_sf = ', rl_sf
ENDIF
! !* Constants initialization
rc02 = rc0 * rc0 ; rc02r = 1. / rc02
rc03 = rc02 * rc0
rc04 = rc03 * rc0
rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking
rsbc_tke2 = rn_Dt * rn_crban / rl_sf ! Neumann + Wave breaking
zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer
rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness
rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness
rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi
rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
!
rfact_tke = -0.5_wp / rsc_tke * rn_Dt ! Cst used for the Diffusion term of tke
rfact_psi = -0.5_wp / rsc_psi * rn_Dt ! Cst used for the Diffusion term of tke
!
! !* Wall proximity function
!!gm tmask or wmask ????
zwall(:,:,:) = 1._wp * tmask(:,:,:)
! !* read or initialize all required files
CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n)
!
END SUBROUTINE zdf_gls_init
SUBROUTINE gls_rst( kt, cdrw )
!!---------------------------------------------------------------------
!! *** ROUTINE gls_rst ***
!!
!! ** Purpose : Read or write TKE file (en) in restart file
!!
!! ** Method : use of IOM library
!! if the restart does not contain TKE, en is either
!! set to rn_emin or recomputed (nn_igls/=0)
!!----------------------------------------------------------------------
USE zdf_oce , ONLY : en, avt_k, avm_k ! ocean vertical physics
!!
INTEGER , INTENT(in) :: kt ! ocean time-step
CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
!
INTEGER :: jit, jk ! dummy loop indices
INTEGER :: id1, id2, id3, id4
INTEGER :: ji, jj, ikbu, ikbv
REAL(wp):: cbx, cby
!!----------------------------------------------------------------------
!
IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
! ! ---------------
IF( ln_rstart ) THEN !* Read the restart file
id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
!
IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist
CALL iom_get( numror, jpdom_auto, 'en' , en , kfill = jpfillcopy ) ! we devide by en -> must be != 0.
CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k )
CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k )
CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n, kfill = jpfillcopy ) ! we devide by hmxl_n -> must be != 0.
ELSE
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values'
en (:,:,:) = rn_emin
hmxl_n(:,:,:) = 0.05_wp
! avt_k, avm_k already set to the background value in zdf_phy_init
ENDIF
ELSE !* Start from rest
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>> start from rest, set en and hmxl_n by background values'
en (:,:,:) = rn_emin
hmxl_n(:,:,:) = 0.05_wp
! avt_k, avm_k already set to the background value in zdf_phy_init
ENDIF
!
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
! ! -------------------
IF(lwp) WRITE(numout,*) '---- gls-rst ----'
CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k )
CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k )
CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
!
ENDIF
!
END SUBROUTINE gls_rst
!!======================================================================
END MODULE zdfgls