Skip to content
Snippets Groups Projects
zdfosm.F90 212 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
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 3076 3077 3078 3079 3080 3081 3082 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 3209 3210 3211 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
      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
      !
      ! 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 jk = 1, jpk
         avt(:,:,jk) = avtb(jk) * tmask(:,:,jk)
      END DO
      !
      ! 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(A2D(0)) = posmdia2d(:,:)
            CALL iom_put( cdname, osmdia2d(A2D(nn_hls)) )
         ELSE   ! Halo present
            CALL iom_put( cdname, posmdia2d )
Guillaume Samson's avatar
Guillaume Samson committed
         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(A2D(0),:) = posmdia3d(:,:,:)
            CALL iom_put( cdname, osmdia3d(A2D(nn_hls),:) )
         ELSE   ! Halo present
            CALL iom_put( cdname, posmdia3d )
Guillaume Samson's avatar
Guillaume Samson committed
         END IF
      END IF
      !
   END SUBROUTINE zdf_osm_iomput_3d

   !!======================================================================

END MODULE zdfosm