Ignore:
Timestamp:
2016-11-02T15:24:08+01:00 (4 years ago)
Author:
jcastill
Message:

First trial to changes needed for wave coupling

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7167 r7168  
    2727   PRIVATE 
    2828 
     29   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
    2930   PUBLIC   sbc_wave    ! routine called in sbcmod 
    3031    
    31    INTEGER , PARAMETER ::   jpfld  = 4           ! number of files to read for stokes drift 
    32    INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point 
    33    INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point 
    34    INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave hight      (m)      at T-point 
    35    INTEGER , PARAMETER ::   jp_wmp = 4           ! index of mean wave period            (s)      at T-point 
     32   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
     33   LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE. 
     34   LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE. 
     35   LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE. 
     36   LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE. 
     37   LOGICAL, PUBLIC     ::   cpl_wper=.FALSE. 
     38   LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE. 
     39   LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE. 
     40   LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE. 
     41 
     42   INTEGER             ::   jpfld                ! number of files to read for stokes drift 
     43   INTEGER             ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point 
     44   INTEGER             ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point 
     45   INTEGER             ::   jp_swh               ! index of significant wave hight      (m)      at T-point 
     46   INTEGER             ::   jp_wmp               ! index of mean wave period            (s)      at T-point 
    3647 
    3748   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     
    5768CONTAINS 
    5869 
     70   SUBROUTINE sbc_stokes( ) 
     71      !!--------------------------------------------------------------------- 
     72      !!                     ***  ROUTINE sbc_stokes  *** 
     73      !! 
     74      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al., 
     75      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1) 
     76      !! 
     77      !! ** Method  : - Calculate Stokes transport speed  
     78      !!              - Calculate horizontal divergence  
     79      !!              - Integrate the horizontal divergenze from the bottom  
     80      !! ** action   
     81      !!--------------------------------------------------------------------- 
     82      INTEGER                ::   jj,ji,jk  
     83      REAL(wp)                       ::  ztransp, zsp0, zk, zus, zvs 
     84      REAL(wp), DIMENSION(jpi,jpj)   ::  zfac  
     85      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
     86 
     87 
     88      CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
     89      DO jk = 1, jpk 
     90         DO jj = 1, jpj 
     91            DO ji = 1, jpi 
     92            ! On T grid 
     93            ! Stokes transport speed estimated from Hs and Tmean 
     94            ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
     95            ! Stokes surface speed 
     96            zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 + sf_sd(jp_vsd)%fnow(ji,jj,1)**2) 
     97            ! Wavenumber scale 
     98            zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
     99            ! Depth attenuation 
     100            zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 
     101            END DO 
     102         END DO 
     103 
     104         DO jj = 1, jpjm1 
     105            DO ji = 1, jpim1 
     106               ! Into the U and V Grid  
     107               zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 
     108                  &                                + zfac(ji+1,jj) * tmask(ji+1,jj,1) ) 
     109               zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 
     110                  &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) ) 
     111               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) & 
     112                  &                                         + zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) ) 
     113               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) & 
     114                  &                                         + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) ) 
     115               usd3d(ji,jj,jk) = usd2d(ji,jj) * zus 
     116               vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs 
     117            END DO 
     118         END DO 
     119      END DO 
     120      ! 
     121      CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
     122      CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
     123      ! 
     124      DO jk = 1, jpkm1                       ! e3t * Horizontal divergence 
     125         DO jj = 2, jpjm1 
     126            DO ji = fs_2, fs_jpim1   ! vector opt. 
     127               ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * fse3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     & 
     128                  &                 - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     & 
     129                  &                 + e1v(ji,jj  ) * fse3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     & 
     130                  &                 - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj) 
     131            END DO 
     132         END DO 
     133         IF( .NOT. AGRIF_Root() ) THEN 
     134            IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   : ,jk) = 0._wp      ! east 
     135            IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   : ,jk) = 0._wp      ! west 
     136            IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :  ,nlcj-1,jk) = 0._wp      ! north 
     137            IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2  ,jk) = 0._wp      ! south 
     138         ENDIF 
     139      END DO 
     140      CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
     141      ! 
     142      DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence 
     143         wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 
     144      END DO 
     145#if defined key_bdy 
     146      IF( lk_bdy ) THEN 
     147         DO jk = 1, jpkm1 
     148            wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
     149         END DO 
     150      ENDIF 
     151#endif 
     152      CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 
     153      ! 
     154   END SUBROUTINE sbc_stokes 
     155 
     156   SUBROUTINE sbc_qiao( ) 
     157      !!--------------------------------------------------------------------- 
     158      !!                     ***  ROUTINE sbc_qiao  *** 
     159      !! 
     160      !! ** Purpose :   Qiao formulation for wave enhanced turbulence 
     161      !!                2010 (DOI: 10.1007/s10236-010-0326)  
     162      !! 
     163      !! ** Method  : -  
     164      !! ** action   
     165      !!--------------------------------------------------------------------- 
     166      INTEGER                ::   jj,ji 
     167 
     168      ! Calculate the module of the stokes drift on T grid 
     169      !------------------------------------------------- 
     170      DO jj = 1, jpj 
     171         DO ji = 1, jpi 
     172            tsd2d(ji,jj) = ((zusd2dt(ji,jj) * tmask(ji,jj,1))**2.0  + & 
     173            &               (zvsd2dt(ji,jj) * tmask(ji,jj,1))**2.0)**0.5 
     174         END DO 
     175      END DO 
     176      ! 
     177   END SUBROUTINE sbc_qiao 
     178 
    59179   SUBROUTINE sbc_wave( kt ) 
    60180      !!--------------------------------------------------------------------- 
     
    71191      !! ** action   
    72192      !!--------------------------------------------------------------------- 
    73       USE zdf_oce,  ONLY : ln_zdfqiao 
     193      USE zdf_oce    ,  ONLY : ln_zdfqiao 
     194 
     195      IMPLICIT NONE 
    74196 
    75197      INTEGER, INTENT( in  ) ::   kt       ! ocean time step 
    76198      ! 
    77199      INTEGER                ::   ierror   ! return error code 
    78       INTEGER                ::   ifpr, jj,ji,jk  
     200      INTEGER                ::   ifpr 
    79201      INTEGER                ::   ios      ! Local integer output status for namelist read 
    80202      ! 
    81203      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    82       REAL(wp)                       ::  ztransp, zsp0, zk, zus, zvs 
    83       REAL(wp), DIMENSION(jpi,jpj)   ::  zfac  
    84       REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
    85       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     204      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    86205      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    87206                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     
    103222         ! 
    104223         IF ( ln_cdgw ) THEN 
    105             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    106             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    107             ! 
    108                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    109             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    110             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     224            IF ( .NOT. cpl_wdrag ) THEN 
     225               ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     226               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     227               ! 
     228                                      ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
     229               IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
     230               CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     231            ENDIF 
    111232            ALLOCATE( cdn_wave(jpi,jpj) ) 
    112233            cdn_wave(:,:) = 0.0 
     
    114235 
    115236         IF ( ln_tauoc ) THEN 
    116             ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
    117             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    118             ! 
    119                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
    120             IF( sn_tauoc%ln_tint )   ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
    121             CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     237            IF ( .NOT. cpl_wstrf ) THEN 
     238               ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     239               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     240               ! 
     241                                       ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     242               IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     243               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     244            ENDIF 
    122245            ALLOCATE( tauoc_wave(jpi,jpj) ) 
    123246            tauoc_wave(:,:) = 0.0 
    124         ENDIF 
     247         ENDIF 
    125248 
    126249         IF ( ln_sdw ) THEN 
    127             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; 
    128             slf_i(jp_swh) = sn_swh ; slf_i(jp_wmp) = sn_wmp; 
    129             ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
    130             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    131             ! 
    132             DO ifpr= 1, jpfld 
    133                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    134                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    135             END DO 
    136  
    137             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     250            ! Find out how many fields have to be read from file if not coupled 
     251            jpfld=0 
     252            jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 
     253            IF( .NOT. cpl_sdrftx ) THEN 
     254               jpfld=jpfld+1 
     255               jp_usd=jpfld 
     256            ENDIF 
     257            IF( .NOT. cpl_sdrfty ) THEN 
     258               jpfld=jpfld+1 
     259               jp_vsd=jpfld 
     260            ENDIF 
     261            IF( .NOT. cpl_hsig ) THEN 
     262               jpfld=jpfld+1 
     263               jp_swh=jpfld 
     264            ENDIF 
     265            IF( .NOT. cpl_wper ) THEN 
     266               jpfld=jpfld+1 
     267               jp_wmp=jpfld 
     268            ENDIF 
     269 
     270            ! Read from file only the non-coupled fields  
     271            IF( jpfld > 0 ) THEN 
     272               ALLOCATE( slf_i(jpfld) ) 
     273               IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 
     274               IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 
     275               IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 
     276               IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 
     277               ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
     278               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     279               ! 
     280               DO ifpr= 1, jpfld 
     281                  ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     282                  IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     283               END DO 
     284 
     285               CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     286            ENDIF 
    138287            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) ) 
    139288            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
     
    145294            swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ; 
    146295            IF ( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
    147                ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
    148                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
    149                                       ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
    150                IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
    151                CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     296               IF ( .NOT. cpl_wnum ) THEN 
     297                  ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     298                  IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
     299                                         ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     300                  IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     301                  CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     302               ENDIF 
    152303               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
    153304               wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp 
     
    156307      ENDIF 
    157308      ! 
    158       IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==! 
     309      IF ( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
    159310         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    160311         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    161312      ENDIF 
    162313 
    163       IF ( ln_tauoc ) THEN             !==  Wave induced stress  ==! 
     314      IF ( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
    164315         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
    165316         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     
    168319      IF ( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    169320         ! 
    170          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
    171          swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
    172          wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
    173          zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    174          zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     321         ! Read from file only if the field is not coupled 
     322         IF( jpfld > 0 ) THEN 
     323            CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
     324            IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
     325            IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     326            IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     327            IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     328         ENDIF 
    175329         ! 
     330         ! Read also wave number if needed, so that it is available in coupling routines 
     331         IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
     332            CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing 
     333            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     334         ENDIF 
     335            
    176336         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
    177337         !(DOI: 10.1175/JPO-D-14-0020.1)==!  
    178338         ! 
    179          CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
    180          DO jk = 1, jpk 
    181             DO jj = 1, jpj 
    182                DO ji = 1, jpi 
    183                ! On T grid 
    184                ! Stokes transport speed estimated from Hs and Tmean 
    185                ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
    186                ! Stokes surface speed 
    187                zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 +  sf_sd(jp_vsd)%fnow(ji,jj,1)**2) 
    188                ! Wavenumber scale 
    189                zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
    190                ! Depth attenuation 
    191                zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 
    192                END DO 
    193             END DO 
    194  
    195             DO jj = 1, jpjm1 
    196                DO ji = 1, jpim1 
    197                   ! Into the U and V Grid  
    198                   zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 
    199                  &                                + zfac(ji+1,jj) * tmask(ji+1,jj,1) ) 
    200  
    201                   zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 
    202                  &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) ) 
    203  
    204                   usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) & 
    205                  &                                             +  zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) ) 
    206  
    207                   vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) & 
    208                  &                                              + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) ) 
    209  
    210                   usd3d(ji,jj,jk) = usd2d(ji,jj) * zus 
    211                   vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs 
    212                END DO 
    213             END DO 
    214          END DO 
    215          ! 
    216          CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
    217          CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
    218          ! 
    219          DO jk = 1, jpkm1                       ! e3t * Horizontal divergence 
    220             DO jj = 2, jpjm1 
    221                DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                   ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * fse3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     & 
    223                      &                 - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     & 
    224                      &                 + e1v(ji,jj  ) * fse3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     & 
    225                      &                 - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj) 
    226                END DO   
    227             END DO   
    228             IF( .NOT. AGRIF_Root() ) THEN 
    229                IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,jk) = 0._wp      ! east 
    230                IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,jk) = 0._wp      ! west 
    231                IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,jk) = 0._wp      ! north 
    232                IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,jk) = 0._wp      ! south 
    233             ENDIF 
    234          END DO 
    235          CALL lbc_lnk( ze3hdiv, 'T', 1. )  
    236          ! 
    237          DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence 
    238             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 
    239          END DO 
    240 #if defined key_bdy 
    241          IF( lk_bdy ) THEN 
    242             DO jk = 1, jpkm1 
    243                wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    244             END DO 
    245          ENDIF 
    246 #endif 
    247          CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 
    248  
    249          IF ( ln_zdfqiao ) THEN 
    250             wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
    251             
    252             ! Calculate the module of the stokes drift on T grid 
    253             !------------------------------------------------- 
    254             DO jj = 1, jpj 
    255                DO ji = 1, jpi 
    256                    tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0  +     & 
    257                 &               (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5 
    258                END DO 
    259             END DO 
    260          ENDIF 
    261          ! 
     339         ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 
     340         IF( jpfld == 4 ) THEN 
     341            CALL sbc_stokes() 
     342            IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
     343               CALL sbc_qiao() 
     344            ENDIF 
     345         ENDIF 
    262346      ENDIF 
    263347      ! 
Note: See TracChangeset for help on using the changeset viewer.