Newer
Older
3001
3002
3003
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: Kmm ! Time level
!!
INTEGER :: ios ! Local integer
INTEGER :: ji, jj, jk ! Dummy loop indices
REAL(wp) :: z1_t2
!!
REAL(wp), PARAMETER :: pp_large = -1e10_wp
!!
NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave, nn_osm_wave, &
& ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd, ln_kpprimix, rn_riinfty, &
& rn_difri, ln_convmix, rn_difconv, nn_osm_wave, nn_osm_SD_reduce, &
& ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter
!! Namelist for Fox-Kemper parametrization
NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat, &
& rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
IF(lwm) WRITE ( numond, namzdf_osm )
IF(lwp) THEN ! Control print
WRITE(numout,*)
WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters'
WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la
WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle
WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la
WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd
WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0
WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes
WRITE(numout,*) ' Horizontal average flag nn_ave = ', nn_ave
WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave
SELECT CASE (nn_osm_wave)
CASE(0)
WRITE(numout,*) ' Calculated assuming constant La#=0.3'
CASE(1)
WRITE(numout,*) ' Calculated from Pierson Moskowitz wind-waves'
CASE(2)
WRITE(numout,*) ' Calculated from ECMWF wave fields'
END SELECT
WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce = ', nn_osm_SD_reduce
WRITE(numout,*) ' Fraction of hbl to average SD over/fit'
WRITE(numout,*) ' Exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac
SELECT CASE (nn_osm_SD_reduce)
CASE(0)
WRITE(numout,*) ' No reduction'
CASE(1)
WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL'
CASE(2)
WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL'
END SELECT
WRITE(numout,*) ' Reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter = ', ln_zdfosm_ice_shelter
WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm
WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, &
& 'm^2/s'
WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix
WRITE(numout,*) ' Local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
WRITE(numout,*) ' Maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri
WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix
WRITE(numout,*) ' Diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv
ENDIF
!
! ! Check wave coupling settings !
! ! Further work needed - see ticket #2447 !
IF ( nn_osm_wave == 2 ) THEN
IF (.NOT. ( ln_wave .AND. ln_sdw )) &
& CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' )
END IF
!
! TEMP: Specifically, the issue is that rCdU_bot is accessed on halo points but is no longer defined there. We must remove halo calculations from zdfosm where possible.
IF( ln_tile ) THEN
CALL ctl_warn( 'zdf_osm_init: tiling is not currently working with OSMOSIS- it is disabled' )
ln_tile = .FALSE.
CALL dom_tile_init
ENDIF
!
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132
3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
3169
3170
3171
3172
3173
3174
3175
3176
3177
3178
3179
3180
3181
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
! Flags associated with diagnostic output
IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) ) ln_dia_pyc_shr = .TRUE.
IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE.
!
! Allocate zdfosm arrays
IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
!
IF( ln_osm_mle ) THEN ! Initialise Fox-Kemper parametrization
READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
903 IF( ios /= 0 ) CALL ctl_nam( ios, 'namosm_mle in reference namelist' )
READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
904 IF( ios > 0 ) CALL ctl_nam( ios, 'namosm_mle in configuration namelist' )
IF(lwm) WRITE ( numond, namosm_mle )
!
IF(lwp) THEN ! Namelist print
WRITE(numout,*)
WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
WRITE(numout,*) '~~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namosm_mle : '
WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle
WRITE(numout,*) ' Magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce
WRITE(numout,*) ' Scale of ML front (ML radius of deform.) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, &
& 'm'
WRITE(numout,*) ' Maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', &
& rn_osm_mle_time, 's'
WRITE(numout,*) ' Reference latitude (deg) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, &
& 'deg'
WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c
WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', &
& rn_osm_mle_thresh, 'm^2/s'
WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's'
WRITE(numout,*) ' Switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit
WRITE(numout,*) ' hmle limit (fraction of zmld) (ln_osm_hmle_limit = .T.) rn_osm_hmle_limit = ', rn_osm_hmle_limit
END IF
END IF
!
IF(lwp) THEN
WRITE(numout,*)
IF ( ln_osm_mle ) THEN
WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation'
ELSE
WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
END IF
END IF
!
IF( ln_osm_mle ) THEN ! MLE initialisation
!
rb_c = grav * rn_osm_mle_rho_c / rho0 ! Mixed Layer buoyancy criteria
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 '
IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
!
IF( nn_osm_mle == 1 ) THEN
rc_f = rn_osm_mle_ce / ( 5e3_wp * 2.0_wp * omega * SIN( rad * rn_osm_mle_lat ) )
END IF
! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
z1_t2 = 2e-5_wp
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
r1_ft(ji,jj) = MIN( 1.0_wp / ( ABS( ff_t(ji,jj)) + epsln ), ABS( ff_t(ji,jj) ) / z1_t2**2 )
END_2D
! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
!
END IF
!
CALL osm_rst( nit000, Kmm, 'READ' ) ! Read or initialize hbl, dh, hmle
!
IF ( ln_zdfddm ) THEN
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm '
END IF
END IF
!
! Set constants not in namelist
! -----------------------------
IF(lwp) THEN
WRITE(numout,*)
END IF
!
dstokes(:,:) = pp_large
IF (nn_osm_wave == 0) THEN
dstokes(:,:) = rn_osm_dstokes
END IF
!
! Horizontal average : initialization of weighting arrays
! -------------------
SELECT CASE ( nn_ave )
CASE ( 0 ) ! no horizontal average
IF(lwp) WRITE(numout,*) ' no horizontal average on avt'
IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
! Weighting mean arrays etmean
! ( 1 1 )
! avt = 1/4 ( 1 1 )
!
etmean(:,:,:) = 0.0_wp
!
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.0_wp, umask(ji-1,jj, jk) + umask(ji,jj,jk) + &
& vmask(ji, jj-1,jk) + vmask(ji,jj,jk) )
END_3D
CASE ( 1 ) ! horizontal average
IF(lwp) WRITE(numout,*) ' horizontal average on avt'
! Weighting mean arrays etmean
! ( 1/2 1 1/2 )
! avt = 1/8 ( 1 2 1 )
! ( 1/2 1 1/2 )
etmean(:,:,:) = 0.0_wp
!
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
etmean(ji,jj,jk) = tmask(ji, jj,jk) / MAX( 1.0_wp, 2.0_wp * tmask(ji,jj,jk) + &
& 0.5_wp * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) + &
& tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) + &
& 1.0_wp * ( tmask(ji-1,jj, jk) + tmask(ji, jj+1,jk) + &
& tmask(ji, jj-1,jk) + tmask(ji+1,jj, jk) ) )
END_3D
CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
CALL ctl_stop( ctmp1 )
END SELECT
!
! Initialization of vertical eddy coef. to the background value
! -------------------------------------------------------------
DO_3D( 0, 0, 0, 0, 1, jpk )
avt(ji,jj,jk) = avtb(jk) * tmask(ji,jj,jk)
END_3D
3212
3213
3214
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227
3228
3229
3230
3231
3232
3233
3234
3235
3236
3237
3238
3239
3240
3241
3242
3243
3244
3245
3246
3247
3248
3249
3250
3251
3252
3253
3254
3255
3256
3257
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269
3270
3271
3272
3273
3274
3275
3276
3277
3278
3279
3280
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290
3291
3292
3293
3294
3295
3296
3297
3298
3299
3300
3301
3302
3303
3304
3305
3306
3307
3308
3309
3310
3311
3312
3313
3314
3315
3316
3317
3318
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
3331
3332
3333
3334
3335
3336
3337
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
3348
3349
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360
3361
3362
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
3373
3374
3375
3376
3377
3378
3379
3380
3381
3382
3383
3384
3385
3386
3387
3388
3389
3390
3391
3392
3393
3394
3395
3396
3397
3398
3399
3400
3401
3402
3403
3404
3405
3406
3407
3408
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419
3420
3421
3422
3423
3424
3425
3426
3427
3428
3429
3430
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
3443
3444
3445
3446
3447
3448
3449
3450
3451
3452
3453
3454
3455
3456
3457
3458
3459
!
! Zero the surface flux for non local term and osm mixed layer depth
! ------------------------------------------------------------------
ghamt(:,:,:) = 0.0_wp
ghams(:,:,:) = 0.0_wp
ghamu(:,:,:) = 0.0_wp
ghamv(:,:,:) = 0.0_wp
!
IF ( ln_dia_osm ) THEN ! Initialise auxiliary arrays for diagnostic output
osmdia2d(:,:) = 0.0_wp
osmdia3d(:,:,:) = 0.0_wp
END IF
!
END SUBROUTINE zdf_osm_init
SUBROUTINE osm_rst( kt, Kmm, cdrw )
!!---------------------------------------------------------------------
!! *** ROUTINE osm_rst ***
!!
!! ** Purpose : Read or write BL fields in restart file
!!
!! ** Method : use of IOM library. If the restart does not contain
!! required fields, they are recomputed from stratification
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! Ocean time step index
INTEGER , INTENT(in ) :: Kmm ! Ocean time level index (middle)
CHARACTER(len=*), INTENT(in ) :: cdrw ! "READ"/"WRITE" flag
!!
INTEGER :: id1, id2, id3 ! iom enquiry index
INTEGER :: ji, jj, jk ! Dummy loop indices
INTEGER :: iiki, ikt ! Local integer
REAL(wp) :: zhbf ! Tempory scalars
REAL(wp) :: zN2_c ! Local scalar
REAL(wp) :: rho_c = 0.01_wp ! Density criterion for mixed layer depth
INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! Level of mixed-layer depth (pycnocline top)
!!----------------------------------------------------------------------
!
!!-----------------------------------------------------------------------------
! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
!!-----------------------------------------------------------------------------
IF( TRIM(cdrw) == 'READ' .AND. ln_rstart) THEN
id1 = iom_varid( numror, 'wn', ldstop = .FALSE. )
IF( id1 > 0 ) THEN ! 'wn' exists; read
CALL iom_get( numror, jpdom_auto, 'wn', ww )
WRITE(numout,*) ' ===>>>> : wn read from restart file'
ELSE
ww(:,:,:) = 0.0_wp
WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
END IF
!
id1 = iom_varid( numror, 'hbl', ldstop = .FALSE. )
id2 = iom_varid( numror, 'dh', ldstop = .FALSE. )
IF( id1 > 0 .AND. id2 > 0 ) THEN ! 'hbl' exists; read and return
CALL iom_get( numror, jpdom_auto, 'hbl', hbl )
CALL iom_get( numror, jpdom_auto, 'dh', dh )
hml(:,:) = hbl(:,:) - dh(:,:) ! Initialise ML depth
WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file'
IF( ln_osm_mle ) THEN
id3 = iom_varid( numror, 'hmle', ldstop = .FALSE. )
IF( id3 > 0 ) THEN
CALL iom_get( numror, jpdom_auto, 'hmle', hmle )
WRITE(numout,*) ' ===>>>> : hmle read from restart file'
ELSE
WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl'
hmle(:,:) = hbl(:,:) ! Initialise MLE depth
END IF
END IF
RETURN
ELSE ! 'hbl' & 'dh' not in restart file, recalculate
WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
END IF
END IF
!
!!-----------------------------------------------------------------------------
! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
!!-----------------------------------------------------------------------------
IF ( TRIM(cdrw) == 'WRITE' ) THEN
IF(lwp) WRITE(numout,*) '---- osm-rst ----'
CALL iom_rstput( kt, nitrst, numrow, 'wn', ww )
CALL iom_rstput( kt, nitrst, numrow, 'hbl', hbl )
CALL iom_rstput( kt, nitrst, numrow, 'dh', dh )
IF ( ln_osm_mle ) THEN
CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle )
END IF
RETURN
END IF
!
!!-----------------------------------------------------------------------------
! Getting hbl, no restart file with hbl, so calculate from surface stratification
!!-----------------------------------------------------------------------------
IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
! w-level of the mixing and mixed layers
CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )
CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )
imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point
hbl(:,:) = 0.0_wp ! Here hbl used as a dummy variable, integrating vertically N^2
zN2_c = grav * rho_c * r1_rho0 ! Convert density criteria into N^2 criteria
!
hbl(:,:) = 0.0_wp ! Here hbl used as a dummy variable, integrating vertically N^2
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
ikt = mbkt(ji,jj)
hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0.0_wp ) * e3w(ji,jj,jk,Kmm)
IF ( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
END_3D
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
iiki = MAX( 4, imld_rst(ji,jj) )
hbl(ji,jj) = gdepw(ji,jj,iiki,Kmm ) ! Turbocline depth
dh(ji,jj) = e3t(ji,jj,iiki-1,Kmm ) ! Turbocline depth
hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
END_2D
!
WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
!
IF( ln_osm_mle ) THEN
hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
END IF
!
ww(:,:,:) = 0.0_wp
WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
!
END SUBROUTINE osm_rst
SUBROUTINE tra_osm( kt, Kmm, pts, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_osm ***
!!
!! ** Purpose : compute and add to the tracer trend the non-local tracer flux
!!
!! ** Method : ???
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! Time step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! Time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! Active tracers and RHS of tracer equation
!!
INTEGER :: ji, jj, jk
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
!!----------------------------------------------------------------------
!
IF ( kt == nit000 ) THEN
IF ( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
IF(lwp) WRITE(numout,*) '~~~~~~~ '
END IF
END IF
!
IF ( l_trdtra ) THEN ! Save ta and sa trends
ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
END IF
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
& - ( ghamt(ji,jj,jk ) &
& - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm)
pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) &
& - ( ghams(ji,jj,jk ) &
& - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
END_3D
!
IF ( l_trdtra ) THEN ! Save the non-local tracer flux trends for diagnostics
ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt )
CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds )
DEALLOCATE( ztrdt, ztrds )
END IF
!
IF ( sn_cfctl%l_prtctl ) THEN
CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, &
& tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
END IF
!
END SUBROUTINE tra_osm
SUBROUTINE trc_osm( kt ) ! Dummy routine
!!----------------------------------------------------------------------
!! *** ROUTINE trc_osm ***
!!
!! ** Purpose : compute and add to the passive tracer trend the non-local
!! passive tracer flux
!!
!!
!! ** Method : ???
!!
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt
!!----------------------------------------------------------------------
!
WRITE(*,*) 'trc_osm: Not written yet', kt
!
END SUBROUTINE trc_osm
SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_osm ***
!!
!! ** Purpose : compute and add to the velocity trend the non-local flux
!! copied/modified from tra_osm
!!
!! ** Method : ???
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! Ocean time step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! Ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! Ocean velocities and RHS of momentum equation
!!
INTEGER :: ji, jj, jk ! dummy loop indices
!!----------------------------------------------------------------------
!
IF ( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
IF(lwp) WRITE(numout,*) '~~~~~~~ '
END IF
!
! Code saving tracer trends removed, replace with trdmxl_oce
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Add non-local u and v fluxes
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( ghamu(ji,jj,jk ) - &
& ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( ghamv(ji,jj,jk ) - &
& ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm)
END_3D
!
! Code for saving tracer trends removed
!
END SUBROUTINE dyn_osm
SUBROUTINE zdf_osm_iomput_2d( cdname, posmdia2d )
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_osm_iomput_2d ***
!!
!! ** Purpose : Wrapper for subroutine iom_put that accepts 2D arrays
!! with and without halo
!!
!!----------------------------------------------------------------------
CHARACTER(LEN=*), INTENT(in ) :: cdname
REAL(wp), DIMENSION(:,:), INTENT(in ) :: posmdia2d
!!----------------------------------------------------------------------
!
IF ( ln_dia_osm .AND. iom_use( cdname ) ) THEN
IF ( SIZE( posmdia2d, 1 ) == ntei-ntsi+1 .AND. SIZE( posmdia2d, 2 ) == ntej-ntsj+1 ) THEN ! Halo absent
osmdia2d(T2D(0)) = posmdia2d(:,:)
CALL iom_put( cdname, osmdia2d(T2D(nn_hls)) )
CALL iom_put( cdname, posmdia2d )
END IF
END IF
!
END SUBROUTINE zdf_osm_iomput_2d
SUBROUTINE zdf_osm_iomput_3d( cdname, posmdia3d )
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_osm_iomput_3d ***
!!
!! ** Purpose : Wrapper for subroutine iom_put that accepts 3D arrays
!! with and without halo
!!
!!----------------------------------------------------------------------
CHARACTER(LEN=*), INTENT(in ) :: cdname
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: posmdia3d
!!----------------------------------------------------------------------
!
IF ( ln_dia_osm .AND. iom_use( cdname ) ) THEN
IF ( SIZE( posmdia3d, 1 ) == ntei-ntsi+1 .AND. SIZE( posmdia3d, 2 ) == ntej-ntsj+1 ) THEN ! Halo absent
osmdia3d(T2D(0),:) = posmdia3d(:,:,:)
CALL iom_put( cdname, osmdia3d(T2D(nn_hls),:) )
CALL iom_put( cdname, posmdia3d )
END IF
END IF
!
END SUBROUTINE zdf_osm_iomput_3d
!!======================================================================
END MODULE zdfosm