New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7403 for branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2016-11-30T17:56:53+01:00 (8 years ago)
Author:
timgraham
Message:

Merge dev_INGV_METO_merge_2016 into branch

Location:
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC
Files:
33 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r6665 r7403  
    66   !! History :  3.2  !  2009-11  (S. Masson)  Original code 
    77   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 
    8    !!---------------------------------------------------------------------- 
    9 #if defined key_diaar5 
    10    !!---------------------------------------------------------------------- 
    11    !!   'key_diaar5'  :                           activate ar5 diagnotics 
    128   !!---------------------------------------------------------------------- 
    139   !!   dia_ar5       : AR5 diagnostics 
    … …  
    2420   USE phycst         ! physical constant 
    2521   USE in_out_manager  ! I/O manager 
     22   USE zdfddm 
     23   USE zdf_oce 
    2624 
    2725   IMPLICIT NONE 
    … …  
    2927 
    3028   PUBLIC   dia_ar5        ! routine called in step.F90 module 
    31    PUBLIC   dia_ar5_init   ! routine called in opa.F90 module 
    3229   PUBLIC   dia_ar5_alloc  ! routine called in nemogcm.F90 module 
    33  
    34    LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE.   ! coupled flag 
     30   PUBLIC   dia_ar5_hst    ! heat/salt transport 
    3531 
    3632   REAL(wp)                         ::   vol0         ! ocean volume (interior domain) 
    … …  
    3935   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain) 
    4036   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity 
     37 
     38   LOGICAL  :: l_ar5 
    4139       
     40   !! * Substitutions 
     41#  include "zdfddm_substitute.h90" 
     42#  include "vectopt_loop_substitute.h90" 
    4243   !!---------------------------------------------------------------------- 
    4344   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    7374      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
    7475      REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass 
     76      REAL(wp) ::   zaw, zbw, zrw 
    7577      ! 
    7678      REAL(wp), POINTER, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace  
     79      REAL(wp), POINTER, DIMENSION(:,:)     :: zpe                         ! 2D workspace  
    7780      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zrhd , zrhop               ! 3D workspace 
    7881      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace 
    … …  
    8083      IF( nn_timing == 1 )   CALL timing_start('dia_ar5') 
    8184  
    82       CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres ) 
    83       CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    ) 
    84       CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 ) 
    85  
    86       zarea_ssh(:,:) = area(:,:) * sshn(:,:) 
    87  
    88       !                                         ! total volume of liquid seawater 
    89       zvolssh = SUM( zarea_ssh(:,:) )  
    90       IF( lk_mpp )   CALL mpp_sum( zvolssh ) 
    91       zvol = vol0 + zvolssh 
     85      IF( kt == nit000 )     CALL dia_ar5_init 
     86 
     87      IF( l_ar5 ) THEN 
     88         CALL wrk_alloc( jpi , jpj              , zarea_ssh , zbotpres ) 
     89         CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    ) 
     90         CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 ) 
     91         zarea_ssh(:,:) = area(:,:) * sshn(:,:) 
     92      ENDIF 
     93      ! 
     94      IF( iom_use( 'voltot' ) .OR. iom_use( 'sshtot' )  .OR. iom_use( 'sshdyn' )  ) THEN     
     95         !                                         ! total volume of liquid seawater 
     96         zvolssh = SUM( zarea_ssh(:,:) )  
     97         IF( lk_mpp )   CALL mpp_sum( zvolssh ) 
     98         zvol = vol0 + zvolssh 
    9299       
    93       CALL iom_put( 'voltot', zvol               ) 
    94       CALL iom_put( 'sshtot', zvolssh / area_tot ) 
    95  
    96       !                      
    97       ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh 
    98       ztsn(:,:,:,jp_sal) = sn0(:,:,:) 
    99       CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity 
    100       ! 
    101       zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    102       DO jk = 1, jpkm1 
    103          zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 
    104       END DO 
    105       IF( ln_linssh ) THEN 
    106          IF( ln_isfcav ) THEN 
    107             DO ji=1,jpi 
    108                DO jj=1,jpj 
    109                   zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    110                END DO 
    111             END DO 
    112          ELSE 
    113             zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
    114          END IF 
     100         CALL iom_put( 'voltot', zvol               ) 
     101         CALL iom_put( 'sshtot', zvolssh / area_tot ) 
     102         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) ) 
     103         ! 
     104      ENDIF 
     105 
     106      IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN     
     107         !                      
     108         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh 
     109         ztsn(:,:,:,jp_sal) = sn0(:,:,:) 
     110         CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity 
     111         ! 
     112         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     113         DO jk = 1, jpkm1 
     114            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 
     115         END DO 
     116         IF( ln_linssh ) THEN 
     117            IF( ln_isfcav ) THEN 
     118               DO ji = 1, jpi 
     119                  DO jj = 1, jpj 
     120                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     121                  END DO 
     122               END DO 
     123            ELSE 
     124               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     125            END IF 
    115126!!gm 
    116127!!gm   riceload should be added in both ln_linssh=T or F, no? 
    117128!!gm 
    118       END IF 
    119       !                                          
    120       zarho = SUM( area(:,:) * zbotpres(:,:) )  
    121       IF( lk_mpp )   CALL mpp_sum( zarho ) 
    122       zssh_steric = - zarho / area_tot 
    123       CALL iom_put( 'sshthster', zssh_steric ) 
     129         END IF 
     130         !                                          
     131         zarho = SUM( area(:,:) * zbotpres(:,:) )  
     132         IF( lk_mpp )   CALL mpp_sum( zarho ) 
     133         zssh_steric = - zarho / area_tot 
     134         CALL iom_put( 'sshthster', zssh_steric ) 
    124135       
    125       !                                         ! steric sea surface height 
    126       CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density 
    127       zrhop(:,:,jpk) = 0._wp 
    128       CALL iom_put( 'rhop', zrhop ) 
    129       ! 
    130       zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     136         !                                         ! steric sea surface height 
     137         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density 
     138         zrhop(:,:,jpk) = 0._wp 
     139         CALL iom_put( 'rhop', zrhop ) 
     140         ! 
     141         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     142         DO jk = 1, jpkm1 
     143            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 
     144         END DO 
     145         IF( ln_linssh ) THEN 
     146            IF ( ln_isfcav ) THEN 
     147               DO ji = 1,jpi 
     148                  DO jj = 1,jpj 
     149                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     150                  END DO 
     151               END DO 
     152            ELSE 
     153               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     154            END IF 
     155         END IF 
     156         !     
     157         zarho = SUM( area(:,:) * zbotpres(:,:) )  
     158         IF( lk_mpp )   CALL mpp_sum( zarho ) 
     159         zssh_steric = - zarho / area_tot 
     160         CALL iom_put( 'sshsteric', zssh_steric ) 
     161       
     162         !                                         ! ocean bottom pressure 
     163         zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
     164         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 
     165         CALL iom_put( 'botpres', zbotpres ) 
     166         ! 
     167      ENDIF 
     168 
     169      IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' )  .OR. iom_use( 'saltot' )  ) THEN     
     170         !                                         ! Mean density anomalie, temperature and salinity 
     171         ztemp = 0._wp 
     172         zsal  = 0._wp 
     173         DO jk = 1, jpkm1 
     174            DO jj = 1, jpj 
     175               DO ji = 1, jpi 
     176                  zztmp = area(ji,jj) * e3t_n(ji,jj,jk) 
     177                  ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem) 
     178                  zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal) 
     179               END DO 
     180            END DO 
     181         END DO 
     182         IF( ln_linssh ) THEN 
     183            IF( ln_isfcav ) THEN 
     184               DO ji = 1, jpi 
     185                  DO jj = 1, jpj 
     186                     ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
     187                     zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
     188                  END DO 
     189               END DO 
     190            ELSE 
     191               ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) ) 
     192               zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) ) 
     193            END IF 
     194         ENDIF 
     195         IF( lk_mpp ) THEN   
     196            CALL mpp_sum( ztemp ) 
     197            CALL mpp_sum( zsal  ) 
     198         END IF 
     199         ! 
     200         zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater 
     201         ztemp = ztemp / zvol                            ! potential temperature in liquid seawater 
     202         zsal  = zsal  / zvol                            ! Salinity of liquid seawater 
     203         ! 
     204         CALL iom_put( 'masstot', zmass ) 
     205         CALL iom_put( 'temptot', ztemp ) 
     206         CALL iom_put( 'saltot' , zsal  ) 
     207         ! 
     208      ENDIF 
     209 
     210      IF( iom_use( 'tnpeo' )) THEN     
     211      ! Work done against stratification by vertical mixing 
     212      ! Exclude points where rn2 is negative as convection kicks in here and 
     213      ! work is not being done against stratification 
     214          CALL wrk_alloc( jpi, jpj, zpe ) 
     215          zpe(:,:) = 0._wp 
     216          IF( lk_zdfddm ) THEN 
     217             DO ji=1,jpi 
     218                DO jj=1,jpj 
     219                   DO jk=1,jpk 
     220                      zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
     221                         &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
     222                      ! 
     223                      zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw 
     224                      zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw 
     225                      ! 
     226                      zpe(ji, jj) = zpe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * & 
     227                           &       grav * (avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  & 
     228                           &       - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) ) 
     229 
     230                   ENDDO 
     231                ENDDO 
     232             ENDDO 
     233          ELSE 
     234             DO ji = 1, jpi 
     235                DO jj = 1, jpj 
     236                   DO jk = 1, jpk 
     237                       zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk) 
     238                   ENDDO 
     239                ENDDO 
     240             ENDDO 
     241          ENDIF 
     242          CALL lbc_lnk( zpe, 'T', 1._wp)          
     243          CALL iom_put( 'tnpeo', zpe ) 
     244          CALL wrk_dealloc( jpi, jpj, zpe ) 
     245      ENDIF 
     246      ! 
     247      IF( l_ar5 ) THEN 
     248        CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres ) 
     249        CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    ) 
     250        CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 ) 
     251      ENDIF 
     252      ! 
     253      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5') 
     254      ! 
     255   END SUBROUTINE dia_ar5 
     256 
     257   SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva )  
     258      !!---------------------------------------------------------------------- 
     259      !!                    ***  ROUTINE dia_ar5_htr *** 
     260      !!---------------------------------------------------------------------- 
     261      !! Wrapper for heat transport calculations 
     262      !! Called from all advection and/or diffusion routines 
     263      !!---------------------------------------------------------------------- 
     264      INTEGER                         , INTENT(in )  :: ktra  ! tracer index 
     265      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf' 
     266      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion 
     267      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion 
     268      ! 
     269      INTEGER    ::  ji, jj, jk 
     270      REAL(wp), POINTER, DIMENSION(:,:)  :: z2d 
     271 
     272     
     273 
     274      CALL wrk_alloc( jpi, jpj, z2d ) 
     275      z2d(:,:) = pua(:,:,1)  
    131276      DO jk = 1, jpkm1 
    132          zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 
    133       END DO 
    134       IF( ln_linssh ) THEN 
    135          IF ( ln_isfcav ) THEN 
    136             DO ji=1,jpi 
    137                DO jj=1,jpj 
    138                   zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    139                END DO 
     277         DO jj = 2, jpjm1 
     278            DO ji = fs_2, fs_jpim1   ! vector opt. 
     279               z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk)  
    140280            END DO 
    141          ELSE 
    142             zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
    143          END IF 
    144       END IF 
    145       !     
    146       zarho = SUM( area(:,:) * zbotpres(:,:) )  
    147       IF( lk_mpp )   CALL mpp_sum( zarho ) 
    148       zssh_steric = - zarho / area_tot 
    149       CALL iom_put( 'sshsteric', zssh_steric ) 
    150        
    151       !                                         ! ocean bottom pressure 
    152       zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
    153       zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 
    154       CALL iom_put( 'botpres', zbotpres ) 
    155  
    156       !                                         ! Mean density anomalie, temperature and salinity 
    157       ztemp = 0._wp 
    158       zsal  = 0._wp 
    159       DO jk = 1, jpkm1 
    160          DO jj = 1, jpj 
    161             DO ji = 1, jpi 
    162                zztmp = area(ji,jj) * e3t_n(ji,jj,jk) 
    163                ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem) 
    164                zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal) 
    165             END DO 
    166          END DO 
    167       END DO 
    168       IF( ln_linssh ) THEN 
    169          IF( ln_isfcav ) THEN 
    170             DO ji=1,jpi 
    171                DO jj=1,jpj 
    172                   ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
    173                   zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
    174                END DO 
    175             END DO 
    176          ELSE 
    177             ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) ) 
    178             zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) ) 
    179          END IF 
    180       ENDIF 
    181       IF( lk_mpp ) THEN   
    182          CALL mpp_sum( ztemp ) 
    183          CALL mpp_sum( zsal  ) 
    184       END IF 
    185       ! 
    186       zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater 
    187       ztemp = ztemp / zvol                            ! potential temperature in liquid seawater 
    188       zsal  = zsal  / zvol                            ! Salinity of liquid seawater 
    189       ! 
    190       CALL iom_put( 'masstot', zmass ) 
    191       CALL iom_put( 'temptot', ztemp ) 
    192       CALL iom_put( 'saltot' , zsal  ) 
    193       ! 
    194       CALL wrk_dealloc( jpi , jpj              , zarea_ssh , zbotpres ) 
    195       CALL wrk_dealloc( jpi , jpj , jpk        , zrhd      , zrhop    ) 
    196       CALL wrk_dealloc( jpi , jpj , jpk , jpts , ztsn                 ) 
    197       ! 
    198       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5') 
    199       ! 
    200    END SUBROUTINE dia_ar5 
     281         END DO 
     282       END DO 
     283       CALL lbc_lnk( z2d, 'U', -1. ) 
     284       IF( cptr == 'adv' ) THEN 
     285          IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in i-direction 
     286          IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rau0     * z2d )  ! advective salt transport in i-direction 
     287       ENDIF 
     288       IF( cptr == 'ldf' ) THEN 
     289          IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in i-direction 
     290          IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rau0     * z2d ) ! diffusive salt transport in i-direction 
     291       ENDIF 
     292       ! 
     293       z2d(:,:) = pva(:,:,1)  
     294       DO jk = 1, jpkm1 
     295          DO jj = 2, jpjm1 
     296             DO ji = fs_2, fs_jpim1   ! vector opt. 
     297                z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk)  
     298             END DO 
     299          END DO 
     300       END DO 
     301       CALL lbc_lnk( z2d, 'V', -1. ) 
     302       IF( cptr == 'adv' ) THEN 
     303          IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in j-direction 
     304          IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rau0     * z2d )  ! advective salt transport in j-direction 
     305       ENDIF 
     306       IF( cptr == 'ldf' ) THEN 
     307          IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in j-direction 
     308          IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rau0     * z2d ) ! diffusive salt transport in j-direction 
     309       ENDIF 
     310           
     311       CALL wrk_dealloc( jpi, jpj, z2d ) 
     312 
     313   END SUBROUTINE dia_ar5_hst 
    201314 
    202315 
    … …  
    217330      IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init') 
    218331      ! 
    219       CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta ) 
    220       !                                      ! allocate dia_ar5 arrays 
    221       IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' ) 
    222  
    223       area(:,:) = e1e2t(:,:) * tmask_i(:,:) 
    224  
    225       area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot ) 
    226  
    227       vol0        = 0._wp 
    228       thick0(:,:) = 0._wp 
    229       DO jk = 1, jpkm1 
    230          vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) ) 
    231          thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) 
    232       END DO 
    233       IF( lk_mpp )   CALL mpp_sum( vol0 ) 
    234  
    235  
    236       CALL iom_open ( 'sali_ref_clim_monthly', inum ) 
    237       CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  ) 
    238       CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 ) 
    239       CALL iom_close( inum ) 
    240  
    241       sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
    242       sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
    243       IF( ln_zps ) THEN               ! z-coord. partial steps 
    244          DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    245             DO ji = 1, jpi 
    246                ik = mbkt(ji,jj) 
    247                IF( ik > 1 ) THEN 
    248                   zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    249                   sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
    250                ENDIF 
     332      l_ar5 = .FALSE. 
     333      IF(   iom_use( 'voltot'  ) .OR. iom_use( 'sshtot'    )  .OR. iom_use( 'sshdyn' )  .OR.  &  
     334         &  iom_use( 'masstot' ) .OR. iom_use( 'temptot'   )  .OR. iom_use( 'saltot' ) .OR.  &     
     335         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) L_ar5 = .TRUE. 
     336   
     337      IF( l_ar5 ) THEN 
     338         ! 
     339         CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta ) 
     340         !                                      ! allocate dia_ar5 arrays 
     341         IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' ) 
     342 
     343         area(:,:) = e1e2t(:,:) * tmask_i(:,:) 
     344 
     345         area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot ) 
     346 
     347         vol0        = 0._wp 
     348         thick0(:,:) = 0._wp 
     349         DO jk = 1, jpkm1 
     350            vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) ) 
     351            thick0(:,:) = thick0(:,:) +    tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) 
     352         END DO 
     353         IF( lk_mpp )   CALL mpp_sum( vol0 ) 
     354 
     355 
     356         CALL iom_open ( 'sali_ref_clim_monthly', inum ) 
     357         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  ) 
     358         CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 ) 
     359         CALL iom_close( inum ) 
     360 
     361         sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
     362         sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
     363         IF( ln_zps ) THEN               ! z-coord. partial steps 
     364            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
     365               DO ji = 1, jpi 
     366                  ik = mbkt(ji,jj) 
     367                  IF( ik > 1 ) THEN 
     368                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     369                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
     370                  ENDIF 
     371               END DO 
    251372            END DO 
    252          END DO 
    253       ENDIF 
    254       ! 
    255       CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta ) 
     373         ENDIF 
     374         ! 
     375         CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta ) 
     376         ! 
     377      ENDIF 
    256378      ! 
    257379      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init') 
    258380      ! 
    259381   END SUBROUTINE dia_ar5_init 
    260  
    261 #else 
    262    !!---------------------------------------------------------------------- 
    263    !!   Default option :                                         NO diaar5 
    264    !!---------------------------------------------------------------------- 
    265    LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .FALSE.   ! coupled flag 
    266 CONTAINS 
    267    SUBROUTINE dia_ar5_init    ! Dummy routine 
    268    END SUBROUTINE dia_ar5_init 
    269    SUBROUTINE dia_ar5( kt )   ! Empty routine 
    270       INTEGER ::   kt 
    271       WRITE(*,*) 'dia_ar5: You should not have seen this print! error?', kt 
    272    END SUBROUTINE dia_ar5 
    273 #endif 
    274382 
    275383   !!====================================================================== 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r6140 r7403  
    99   !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation 
    1010   !!            3.6  ! 2014-12  (C. Ethe) use of IOM 
     11   !!            3.6  ! 2016-06  (T. Graham) Addition of diagnostics for CMIP6 
    1112   !!---------------------------------------------------------------------- 
    1213 
    … …  
    3839   PUBLIC   dia_ptr_init   ! call in step module 
    3940   PUBLIC   dia_ptr        ! call in step module 
     41   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines 
    4042 
    4143   !                                  !!** namelist  namptr  ** 
    42    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf   !: Heat TRansports (adv, diff, overturn.) 
    43    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf   !: Salt TRansports (adv, diff, overturn.) 
    44     
     44   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_adv, htr_ldf, htr_eiv   !: Heat TRansports (adv, diff, Bolus.) 
     45   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   str_adv, str_ldf, str_eiv   !: Salt TRansports (adv, diff, Bolus.) 
     46   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_ove, str_ove   !: heat Salt TRansports ( overturn.) 
     47   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   htr_btr, str_btr   !: heat Salt TRansports ( barotropic ) 
    4548 
    4649   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F) 
    4750   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation 
    48    INTEGER        ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)  
     51   INTEGER, PUBLIC ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)  
    4952 
    5053   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup 
    … …  
    7578      ! 
    7679      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    77       REAL(wp) ::   zv, zsfc               ! local scalar 
     80      REAL(wp) ::   zsfc,zvfc               ! local scalar 
    7881      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace 
    7982      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace 
    8083      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace 
    8184      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace 
    82       CHARACTER( len = 10 )  :: cl1 
     85      REAL(wp), DIMENSION(jpj)     ::  vsum   ! 1D workspace 
     86      REAL(wp), DIMENSION(jpj,jpts)     ::  tssum   ! 1D workspace 
     87  
     88      ! 
     89      !overturning calculation 
     90      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse 
     91      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace 
     93 
     94 
     95      CHARACTER( len = 12 )  :: cl1 
    8396      !!---------------------------------------------------------------------- 
    8497      ! 
    … …  
    109122            END DO 
    110123         ENDIF 
     124         IF( iom_use("sopstove") .OR. iom_use("sophtove") .OR. iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN 
     125            ! define fields multiplied by scalar 
     126            zmask(:,:,:) = 0._wp 
     127            zts(:,:,:,:) = 0._wp 
     128            zvn(:,:,:) = 0._wp 
     129            DO jk = 1, jpkm1 
     130               DO jj = 1, jpjm1 
     131                  DO ji = 1, jpi 
     132                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk) 
     133                     zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc 
     134                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid 
     135                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc 
     136                     zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc 
     137                  ENDDO 
     138               ENDDO 
     139             ENDDO 
     140         ENDIF 
     141         IF( iom_use("sopstove") .OR. iom_use("sophtove") ) THEN 
     142             sjk(:,:,1) = ptr_sjk( zmask(:,:,:), btmsk(:,:,1) ) 
     143             r1_sjk(:,:,1) = 0._wp 
     144             WHERE( sjk(:,:,1) /= 0._wp )   r1_sjk(:,:,1) = 1._wp / sjk(:,:,1) 
     145 
     146             ! i-mean T and S, j-Stream-Function, global 
     147             tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1) 
     148             sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1) 
     149             v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) ) 
     150 
     151             htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 ) 
     152             str_ove(:,1) = SUM( v_msf(:,:,1)*sn_jk(:,:,1) ,2 ) 
     153 
     154             z2d(1,:) = htr_ove(:,1) * rc_pwatt        !  (conversion in PW) 
     155             DO ji = 1, jpi 
     156               z2d(ji,:) = z2d(1,:) 
     157             ENDDO 
     158             cl1 = 'sophtove' 
     159             CALL iom_put( TRIM(cl1), z2d ) 
     160             z2d(1,:) = str_ove(:,1) * rc_ggram        !  (conversion in Gg) 
     161             DO ji = 1, jpi 
     162               z2d(ji,:) = z2d(1,:) 
     163             ENDDO 
     164             cl1 = 'sopstove' 
     165             CALL iom_put( TRIM(cl1), z2d ) 
     166             IF( ln_subbas ) THEN 
     167                DO jn = 2, nptr 
     168                    sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
     169                    r1_sjk(:,:,jn) = 0._wp 
     170                    WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn) 
     171 
     172                    ! i-mean T and S, j-Stream-Function, basin 
     173                    tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
     174                    sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
     175                    v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) )  
     176                    htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 ) 
     177                    str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 ) 
     178 
     179                    z2d(1,:) = htr_ove(:,jn) * rc_pwatt !  (conversion in PW) 
     180                    DO ji = 1, jpi 
     181                        z2d(ji,:) = z2d(1,:) 
     182                    ENDDO 
     183                    cl1 = TRIM('sophtove_'//clsubb(jn)) 
     184                    CALL iom_put( cl1, z2d ) 
     185                    z2d(1,:) = str_ove(:,jn) * rc_ggram        ! (conversion in Gg) 
     186                    DO ji = 1, jpi 
     187                        z2d(ji,:) = z2d(1,:) 
     188                    ENDDO 
     189                    cl1 = TRIM('sopstove_'//clsubb(jn)) 
     190                    CALL iom_put( cl1, z2d ) 
     191                END DO 
     192             ENDIF 
     193         ENDIF 
     194         IF( iom_use("sopstbtr") .OR. iom_use("sophtbtr") ) THEN 
     195         ! Calculate barotropic heat and salt transport here  
     196             sjk(:,1,1) = ptr_sj( zmask(:,:,:), btmsk(:,:,1) ) 
     197             r1_sjk(:,1,1) = 0._wp 
     198             WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1) 
     199             
     200            vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1)) 
     201            tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) ) 
     202            tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) ) 
     203            htr_btr(:,1) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,1) 
     204            str_btr(:,1) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,1) 
     205            z2d(1,:) = htr_btr(:,1) * rc_pwatt        !  (conversion in PW) 
     206            DO ji = 2, jpi 
     207               z2d(ji,:) = z2d(1,:) 
     208            ENDDO 
     209            cl1 = 'sophtbtr' 
     210            CALL iom_put( TRIM(cl1), z2d ) 
     211            z2d(1,:) = str_btr(:,1) * rc_ggram        !  (conversion in Gg) 
     212            DO ji = 2, jpi 
     213              z2d(ji,:) = z2d(1,:) 
     214            ENDDO 
     215            cl1 = 'sopstbtr' 
     216            CALL iom_put( TRIM(cl1), z2d ) 
     217            IF( ln_subbas ) THEN 
     218                DO jn = 2, nptr 
     219                    sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) ) 
     220                    r1_sjk(:,1,jn) = 0._wp 
     221                    WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn) 
     222                    vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn)) 
     223                    tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) 
     224                    tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) 
     225                    htr_btr(:,jn) = vsum * tssum(:,jp_tem) * r1_sjk(:,1,jn) 
     226                    str_btr(:,jn) = vsum * tssum(:,jp_sal) * r1_sjk(:,1,jn) 
     227                    z2d(1,:) = htr_btr(:,jn) * rc_pwatt !  (conversion in PW) 
     228                    DO ji = 1, jpi 
     229                        z2d(ji,:) = z2d(1,:) 
     230                    ENDDO 
     231                    cl1 = TRIM('sophtbtr_'//clsubb(jn)) 
     232                    CALL iom_put( cl1, z2d ) 
     233                    z2d(1,:) = str_btr(:,jn) * rc_ggram        ! (conversion in Gg) 
     234                    DO ji = 1, jpi 
     235                        z2d(ji,:) = z2d(1,:) 
     236                    ENDDO 
     237                    cl1 = TRIM('sopstbtr_'//clsubb(jn)) 
     238                    CALL iom_put( cl1, z2d ) 
     239               ENDDO 
     240            ENDIF !ln_subbas 
     241         ENDIF !iom_use("sopstbtr....) 
    111242         ! 
    112243      ELSE 
    … …  
    148279         !                                ! Advective and diffusive heat and salt transport 
    149280         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN    
    150             z2d(1,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW) 
     281            z2d(1,:) = htr_adv(:,1) * rc_pwatt        !  (conversion in PW) 
    151282            DO ji = 1, jpi 
    152283               z2d(ji,:) = z2d(1,:) 
    … …  
    154285            cl1 = 'sophtadv'                  
    155286            CALL iom_put( TRIM(cl1), z2d ) 
    156             z2d(1,:) = str_adv(:) * rc_ggram        ! (conversion in Gg) 
     287            z2d(1,:) = str_adv(:,1) * rc_ggram        ! (conversion in Gg) 
    157288            DO ji = 1, jpi 
    158289               z2d(ji,:) = z2d(1,:) 
    … …  
    160291            cl1 = 'sopstadv' 
    161292            CALL iom_put( TRIM(cl1), z2d ) 
     293            IF( ln_subbas ) THEN 
     294              DO jn=2,nptr 
     295               z2d(1,:) = htr_adv(:,jn) * rc_pwatt        !  (conversion in PW) 
     296               DO ji = 1, jpi 
     297                 z2d(ji,:) = z2d(1,:) 
     298               ENDDO 
     299               cl1 = TRIM('sophtadv_'//clsubb(jn))                  
     300               CALL iom_put( cl1, z2d ) 
     301               z2d(1,:) = str_adv(:,jn) * rc_ggram        ! (conversion in Gg) 
     302               DO ji = 1, jpi 
     303                  z2d(ji,:) = z2d(1,:) 
     304               ENDDO 
     305               cl1 = TRIM('sopstadv_'//clsubb(jn))                  
     306               CALL iom_put( cl1, z2d )               
     307              ENDDO 
     308            ENDIF 
    162309         ENDIF 
    163310         ! 
    164311         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN    
    165             z2d(1,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW)  
     312            z2d(1,:) = htr_ldf(:,1) * rc_pwatt        !  (conversion in PW)  
    166313            DO ji = 1, jpi 
    167314               z2d(ji,:) = z2d(1,:) 
    … …  
    169316            cl1 = 'sophtldf' 
    170317            CALL iom_put( TRIM(cl1), z2d ) 
    171             z2d(1,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg) 
     318            z2d(1,:) = str_ldf(:,1) * rc_ggram        !  (conversion in Gg) 
    172319            DO ji = 1, jpi 
    173320               z2d(ji,:) = z2d(1,:) 
    … …  
    175322            cl1 = 'sopstldf' 
    176323            CALL iom_put( TRIM(cl1), z2d ) 
     324            IF( ln_subbas ) THEN 
     325              DO jn=2,nptr 
     326               z2d(1,:) = htr_ldf(:,jn) * rc_pwatt        !  (conversion in PW) 
     327               DO ji = 1, jpi 
     328                 z2d(ji,:) = z2d(1,:) 
     329               ENDDO 
     330               cl1 = TRIM('sophtldf_'//clsubb(jn))                  
     331               CALL iom_put( cl1, z2d ) 
     332               z2d(1,:) = str_ldf(:,jn) * rc_ggram        ! (conversion in Gg) 
     333               DO ji = 1, jpi 
     334                  z2d(ji,:) = z2d(1,:) 
     335               ENDDO 
     336               cl1 = TRIM('sopstldf_'//clsubb(jn))                  
     337               CALL iom_put( cl1, z2d )               
     338              ENDDO 
     339            ENDIF 
     340         ENDIF 
     341 
     342         IF( iom_use("sophteiv") .OR. iom_use("sopsteiv") ) THEN  
     343            z2d(1,:) = htr_eiv(:,1) * rc_pwatt        !  (conversion in PW)  
     344            DO ji = 1, jpi 
     345               z2d(ji,:) = z2d(1,:) 
     346            ENDDO 
     347            cl1 = 'sophteiv' 
     348            CALL iom_put( TRIM(cl1), z2d ) 
     349            z2d(1,:) = str_eiv(:,1) * rc_ggram        !  (conversion in Gg) 
     350            DO ji = 1, jpi 
     351               z2d(ji,:) = z2d(1,:) 
     352            ENDDO 
     353            cl1 = 'sopsteiv' 
     354            CALL iom_put( TRIM(cl1), z2d ) 
     355            IF( ln_subbas ) THEN 
     356               DO jn=2,nptr 
     357                  z2d(1,:) = htr_eiv(:,jn) * rc_pwatt        !  (conversion in PW) 
     358                  DO ji = 1, jpi 
     359                     z2d(ji,:) = z2d(1,:) 
     360                  ENDDO 
     361                  cl1 = TRIM('sophteiv_'//clsubb(jn))                  
     362                  CALL iom_put( cl1, z2d ) 
     363                  z2d(1,:) = str_eiv(:,jn) * rc_ggram        ! (conversion in Gg) 
     364                  DO ji = 1, jpi 
     365                     z2d(ji,:) = z2d(1,:) 
     366                  ENDDO 
     367                  cl1 = TRIM('sopsteiv_'//clsubb(jn))  
     368                  CALL iom_put( cl1, z2d )               
     369               ENDDO 
     370            ENDIF 
    177371         ENDIF 
    178372         ! 
    … …  
    254448         ! Initialise arrays to zero because diatpr is called before they are first calculated 
    255449         ! Note that this means diagnostics will not be exactly correct when model run is restarted. 
    256          htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp   
    257          htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp  
     450         htr_adv(:,:) = 0._wp  ;  str_adv(:,:) =  0._wp  
     451         htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp  
     452         htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp  
     453         htr_ove(:,:) = 0._wp  ;   str_ove(:,:) =  0._wp 
     454         htr_btr(:,:) = 0._wp  ;   str_btr(:,:) =  0._wp 
    258455         ! 
    259456      ENDIF  
    … …  
    261458   END SUBROUTINE dia_ptr_init 
    262459 
     460   SUBROUTINE dia_ptr_hst( ktra, cptr, pva )  
     461      !!---------------------------------------------------------------------- 
     462      !!                    ***  ROUTINE dia_ptr_hst *** 
     463      !!---------------------------------------------------------------------- 
     464      !! Wrapper for heat and salt transport calculations to calculate them for each basin 
     465      !! Called from all advection and/or diffusion routines 
     466      !!---------------------------------------------------------------------- 
     467      INTEGER                         , INTENT(in )  :: ktra  ! tracer index 
     468      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv' 
     469      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion 
     470      INTEGER                                        :: jn    ! 
     471 
     472      IF( cptr == 'adv' ) THEN 
     473         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pva(:,:,:) ) 
     474         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pva(:,:,:) ) 
     475      ENDIF 
     476      IF( cptr == 'ldf' ) THEN 
     477         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pva(:,:,:) ) 
     478         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pva(:,:,:) ) 
     479      ENDIF 
     480      IF( cptr == 'eiv' ) THEN 
     481         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pva(:,:,:) ) 
     482         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) ) 
     483      ENDIF 
     484      ! 
     485      IF( ln_subbas ) THEN 
     486         ! 
     487         IF( cptr == 'adv' ) THEN 
     488             IF( ktra == jp_tem ) THEN  
     489                DO jn = 2, nptr 
     490                   htr_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     491                END DO 
     492             ENDIF 
     493             IF( ktra == jp_sal ) THEN  
     494                DO jn = 2, nptr 
     495                   str_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     496                END DO 
     497             ENDIF 
     498         ENDIF 
     499         IF( cptr == 'ldf' ) THEN 
     500             IF( ktra == jp_tem ) THEN  
     501                DO jn = 2, nptr 
     502                    htr_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     503                 END DO 
     504             ENDIF 
     505             IF( ktra == jp_sal ) THEN  
     506                DO jn = 2, nptr 
     507                   str_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     508                END DO 
     509             ENDIF 
     510         ENDIF 
     511         IF( cptr == 'eiv' ) THEN 
     512             IF( ktra == jp_tem ) THEN  
     513                DO jn = 2, nptr 
     514                    htr_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     515                 END DO 
     516             ENDIF 
     517             IF( ktra == jp_sal ) THEN  
     518                DO jn = 2, nptr 
     519                   str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     520                END DO 
     521             ENDIF 
     522         ENDIF 
     523         ! 
     524      ENDIF 
     525   END SUBROUTINE dia_ptr_hst 
     526 
    263527 
    264528   FUNCTION dia_ptr_alloc() 
    … …  
    271535      ierr(:) = 0 
    272536      ! 
    273       ALLOCATE( btmsk(jpi,jpj,nptr) ,           & 
    274          &      htr_adv(jpj) , str_adv(jpj) ,   & 
    275          &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  ) 
     537      ALLOCATE( btmsk(jpi,jpj,nptr) ,              & 
     538         &      htr_adv(jpj,nptr) , str_adv(jpj,nptr) ,   & 
     539         &      htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) ,   & 
     540         &      htr_ove(jpj,nptr) , str_ove(jpj,nptr) ,   & 
     541         &      htr_btr(jpj,nptr) , str_btr(jpj,nptr) ,   & 
     542         &      htr_ldf(jpj,nptr) , str_ldf(jpj,nptr) , STAT=ierr(1)  ) 
    276543         ! 
    277544      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2)) 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6387 r7403  
    302302      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence 
    303303      ! 
    304       IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
     304      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
    305305         z3d(:,:,jpk) = 0.e0 
     306         z2d(:,:) = 0.e0 
    306307         DO jk = 1, jpkm1 
    307308            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 
     309            z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    308310         END DO 
    309311         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction 
     312         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum 
    310313      ENDIF 
    311314       
    … …  
    370373         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction 
    371374      ENDIF 
     375 
     376      ! Vertical integral of temperature 
     377      IF( iom_use("tosmint") ) THEN 
     378         z2d(:,:)=0._wp 
     379         DO jk = 1, jpkm1 
     380            DO jj = 2, jpjm1 
     381               DO ji = fs_2, fs_jpim1   ! vector opt. 
     382                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
     383               END DO 
     384            END DO 
     385         END DO 
     386         CALL lbc_lnk( z2d, 'T', -1. ) 
     387         CALL iom_put( "tosmint", z2d )  
     388      ENDIF 
     389 
     390      ! Vertical integral of salinity 
     391      IF( iom_use("somint") ) THEN 
     392         z2d(:,:)=0._wp 
     393         DO jk = 1, jpkm1 
     394            DO jj = 2, jpjm1 
     395               DO ji = fs_2, fs_jpim1   ! vector opt. 
     396                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     397               END DO 
     398            END DO 
     399         END DO 
     400         CALL lbc_lnk( z2d, 'T', -1. ) 
     401         CALL iom_put( "somint", z2d )  
     402      ENDIF 
     403 
     404      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2) 
    372405      ! 
    373406      CALL wrk_dealloc( jpi , jpj      , z2d ) 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90

    r6140 r7403  
    2424   USE ldfslp          ! lateral diffusion: slope of iso-neutral surfaces 
    2525   USE ldfc1d_c2d      ! lateral diffusion: 1D & 2D cases  
    26    USE diaar5, ONLY:   lk_diaar5 
     26   USE diaptr 
    2727   ! 
    2828   USE trc_oce, ONLY: lk_offline ! offline flag 
    … …  
    730730      CALL iom_put( "woce_eiv", zw3d ) 
    731731      ! 
     732      !       
     733      ! 
     734      CALL wrk_alloc( jpi,jpj,   zw2d ) 
     735      ! 
     736      zztmp = 0.5_wp * rau0 * rcp  
     737      IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 
     738        zw2d(:,:)   = 0._wp  
     739        zw3d(:,:,:) = 0._wp  
     740        DO jk = 1, jpkm1 
     741           DO jj = 2, jpjm1 
     742              DO ji = fs_2, fs_jpim1   ! vector opt. 
     743                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   & 
     744                    &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) )  
     745                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     746              END DO 
     747           END DO 
     748        END DO 
     749        CALL lbc_lnk( zw2d, 'U', -1. ) 
     750        CALL lbc_lnk( zw3d, 'U', -1. ) 
     751        CALL iom_put( "ueiv_heattr"  , zztmp * zw2d )                  ! heat transport in i-direction 
     752        CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction 
     753      ENDIF 
     754      zw2d(:,:)   = 0._wp  
     755      zw3d(:,:,:) = 0._wp  
     756      DO jk = 1, jpkm1 
     757         DO jj = 2, jpjm1 
     758            DO ji = fs_2, fs_jpim1   ! vector opt. 
     759               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   & 
     760                  &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) )  
     761               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     762            END DO 
     763         END DO 
     764      END DO 
     765      CALL lbc_lnk( zw2d, 'V', -1. ) 
     766      CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in j-direction 
     767      CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in j-direction 
     768      ! 
     769      IF( ln_diaptr )  CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) 
     770      ! 
     771      zztmp = 0.5_wp * 0.5 
     772      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN 
     773        zw2d(:,:) = 0._wp  
     774        zw3d(:,:,:) = 0._wp  
     775        DO jk = 1, jpkm1 
     776           DO jj = 2, jpjm1 
     777              DO ji = fs_2, fs_jpim1   ! vector opt. 
     778                 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   & 
     779                    &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) )  
     780                 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     781              END DO 
     782           END DO 
     783        END DO 
     784        CALL lbc_lnk( zw2d, 'U', -1. ) 
     785        CALL lbc_lnk( zw3d, 'U', -1. ) 
     786        CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction 
     787        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction 
     788      ENDIF 
     789      zw2d(:,:) = 0._wp  
     790      zw3d(:,:,:) = 0._wp  
     791      DO jk = 1, jpkm1 
     792         DO jj = 2, jpjm1 
     793            DO ji = fs_2, fs_jpim1   ! vector opt. 
     794               zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   & 
     795                  &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) )  
     796               zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     797            END DO 
     798         END DO 
     799      END DO 
     800      CALL lbc_lnk( zw2d, 'V', -1. ) 
     801      CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction 
     802      CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction 
     803      ! 
     804      IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d ) 
     805      ! 
     806      CALL wrk_dealloc( jpi,jpj,   zw2d ) 
    732807      CALL wrk_dealloc( jpi,jpj,jpk,   zw3d ) 
    733       !       
    734       ! 
    735       IF( lk_diaar5 ) THEN                               !==  eiv heat transport: calculate and output  ==! 
    736          CALL wrk_alloc( jpi,jpj,   zw2d ) 
    737          ! 
    738          zztmp = 0.5_wp * rau0 * rcp  
    739          zw2d(:,:) = 0._wp  
    740          DO jk = 1, jpkm1 
    741             DO jj = 2, jpjm1 
    742                DO ji = fs_2, fs_jpim1   ! vector opt. 
    743                   zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   & 
    744                      &                              * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) )  
    745                END DO 
    746             END DO 
    747          END DO 
    748          CALL lbc_lnk( zw2d, 'U', -1. ) 
    749          CALL iom_put( "ueiv_heattr", zw2d )                  ! heat transport in i-direction 
    750          zw2d(:,:) = 0._wp  
    751          DO jk = 1, jpkm1 
    752             DO jj = 2, jpjm1 
    753                DO ji = fs_2, fs_jpim1   ! vector opt. 
    754                   zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   & 
    755                      &                              * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) )  
    756                END DO 
    757             END DO 
    758          END DO 
    759          CALL lbc_lnk( zw2d, 'V', -1. ) 
    760          CALL iom_put( "veiv_heattr", zw2d )                  !  heat transport in i-direction 
    761          ! 
    762          CALL wrk_dealloc( jpi,jpj,   zw2d ) 
    763       ENDIF 
    764808      ! 
    765809      IF( nn_timing == 1 )  CALL timing_stop( 'ldf_eiv_dia')       

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90

    r5836 r7403  
    6666   INTEGER                    ::   nsnd         ! total number of fields sent  
    6767   INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
    68    INTEGER, PUBLIC, PARAMETER ::   nmaxfld=50   ! Maximum number of coupling fields 
     68   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=55   ! Maximum number of coupling fields 
    6969   INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields 
    7070   INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r6140 r7403  
    6565   LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model 
    6666   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
     67   LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used 
     68   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
    6769   ! 
    6870   LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs 
    … …  
    120122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s] 
    121123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1) 
    122 #if defined key_cpl_carbon_cycle 
    123124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    124 #endif 
    125125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
    126126 
    … …  
    166166         ! 
    167167      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
    168 #if defined key_cpl_carbon_cycle 
    169168         &      atm_co2(jpi,jpj) ,                                        & 
    170 #endif 
    171169         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    172170         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r6813 r7403  
    745745 
    746746      !! Neutral coefficients at 10m: 
    747       IF( ln_cdgw ) THEN      ! wave drag case 
     747      IF( ln_wave .AND. ln_cdgw ) THEN      ! wave drag case 
    748748         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
    749749         ztmp0   (:,:) = cdn_wave(:,:) 
    … …  
    791791         END IF 
    792792        
    793          IF( ln_cdgw ) THEN      ! surface wave case 
     793         IF( ln_wave .AND. ln_cdgw ) THEN      ! surface wave case 
    794794            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
    795795            Cd      = sqrt_Cd * sqrt_Cd 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90

    r6140 r7403  
    1717   USE fldread        ! read input fields 
    1818   USE sbc_oce        ! Surface boundary condition: ocean fields 
    19    USE sbcwave  ,ONLY :   cdn_wave !wave module 
    2019   ! 
    2120   USE iom            ! I/O manager library 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6722 r7403  
    1818   !!   sbc_cpl_snd     : send     fields to the atmosphere 
    1919   !!---------------------------------------------------------------------- 
    20    USE dom_oce        ! ocean space and time domain 
    21    USE sbc_oce        ! Surface boundary condition: ocean fields 
    22    USE sbc_ice        ! Surface boundary condition: ice fields 
    23    USE sbcapr         ! Stochastic param. : ??? 
    24    USE sbcdcy         ! surface boundary condition: diurnal cycle 
    25    USE phycst         ! physical constants 
     20   USE dom_oce         ! ocean space and time domain 
     21   USE sbc_oce         ! Surface boundary condition: ocean fields 
     22   USE trc_oce         ! share SMS/Ocean variables 
     23   USE sbc_ice         ! Surface boundary condition: ice fields 
     24   USE sbcapr          ! Stochastic param. : ??? 
     25   USE sbcdcy          ! surface boundary condition: diurnal cycle 
     26   USE sbcwave         ! surface boundary condition: waves 
     27   USE phycst          ! physical constants 
    2628#if defined key_lim3 
    2729   USE ice            ! ice variables 
    … …  
    3638   USE albedo         !  
    3739   USE eosbn2         !  
    38    USE sbcrnf  , ONLY : l_rnfcpl 
    39 #if defined key_cpl_carbon_cycle 
    40    USE p4zflx, ONLY : oce_co2 
    41 #endif 
     40   USE sbcrnf, ONLY : l_rnfcpl 
    4241#if defined key_cice 
    4342   USE ice_domain_size, only: ncat 
    … …  
    106105   INTEGER, PARAMETER ::   jpr_e3t1st = 41   ! first T level thickness  
    107106   INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level 
    108    INTEGER, PARAMETER ::   jprcv      = 42   ! total number of fields received 
     107   INTEGER, PARAMETER ::   jpr_mslp   = 43   ! mean sea level pressure  
     108   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig  
     109   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux  
     110   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1  
     111   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2  
     112   INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period 
     113   INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber 
     114   INTEGER, PARAMETER ::   jpr_wstrf  = 50   ! Stress fraction adsorbed by waves 
     115   INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient 
     116   INTEGER, PARAMETER ::   jprcv      = 51   ! total number of fields received   
    109117 
    110118   INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere 
    … …  
    136144   INTEGER, PARAMETER ::   jps_e3t1st = 27   ! first level depth (vvl) 
    137145   INTEGER, PARAMETER ::   jps_fraqsr = 28   ! fraction of solar net radiation absorbed in the first ocean level 
    138    INTEGER, PARAMETER ::   jpsnd      = 28   ! total number of fields sended 
     146   INTEGER, PARAMETER ::   jps_ficet  = 29   ! total ice fraction   
     147   INTEGER, PARAMETER ::   jps_ocxw   = 30   ! currents on grid 1   
     148   INTEGER, PARAMETER ::   jps_ocyw   = 31   ! currents on grid 2 
     149   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level  
     150   INTEGER, PARAMETER ::   jpsnd      = 32   ! total number of fields sent  
    139151 
    140152   !                                  !!** namelist namsbc_cpl ** 
    … …  
    150162   !                                   ! Received from the atmosphere 
    151163   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 
    152    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                         
     164   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp                            
     165   ! Send to waves  
     166   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev  
     167   ! Received from waves  
     168   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag 
    153169   !                                   ! Other namelist parameters 
    154170   INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
    … …  
    163179   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
    164180 
    165    INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo          ! OASIS info argument 
     181   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]  
     182   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)  
     183 
     184   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument 
    166185 
    167186   !! Substitution 
    … …  
    178197      !!             ***  FUNCTION sbc_cpl_alloc  *** 
    179198      !!---------------------------------------------------------------------- 
    180       INTEGER :: ierr(3) 
     199      INTEGER :: ierr(4) 
    181200      !!---------------------------------------------------------------------- 
    182201      ierr(:) = 0 
    … …  
    189208      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
    190209      ! 
     210      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) )  
     211 
    191212      sbc_cpl_alloc = MAXVAL( ierr ) 
    192213      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc ) 
    … …  
    214235      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    215236      !! 
    216       NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      & 
    217          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
    218          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   & 
    219          &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask 
     237      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      &  
     238         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &  
     239         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &  
     240         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   & 
     241         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   & 
     242         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp  
    220243      !!--------------------------------------------------------------------- 
    221244      ! 
    … …  
    258281         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    259282         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     283         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'  
     284         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'  
     285         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'  
     286         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'  
     287         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'  
     288         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'  
     289         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'  
     290         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'  
    260291         WRITE(numout,*)'  sent fields (multiple ice categories)' 
    261292         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')' 
    262293         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')' 
    263294         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')' 
     295         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')'  
    264296         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')' 
    265297         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref  
    … …  
    267299         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd 
    268300         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')' 
     301         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')'  
     302         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')'  
     303         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'  
     304         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref  
     305         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor  
     306         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd  
    269307         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    270308         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
    … …  
    305343      !  
    306344      ! Vectors: change of sign at north fold ONLY if on the local grid 
     345      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled 
    307346      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 
    308347       
    … …  
    372411         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp. 
    373412      ENDIF 
    374       ! 
     413      ENDIF 
     414 
    375415      !                                                      ! ------------------------- ! 
    376416      !                                                      !    freshwater budget      !   E-P 
    … …  
    467507      !                                                      !      Atmospheric CO2      ! 
    468508      !                                                      ! ------------------------- ! 
    469       srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE. 
     509      srcv(jpr_co2 )%clname = 'O_AtmCO2'    
     510      IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )  THEN 
     511         srcv(jpr_co2 )%laction = .TRUE. 
     512         l_co2cpl = .TRUE. 
     513         IF(lwp) WRITE(numout,*) 
     514         IF(lwp) WRITE(numout,*) '   Atmospheric pco2 received from oasis ' 
     515         IF(lwp) WRITE(numout,*) 
     516      ENDIF 
     517 
     518      !                                                      ! ------------------------- !  
     519      !                                                      ! Mean Sea Level Pressure   !  
     520      !                                                      ! ------------------------- !  
     521      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE.  
     522 
    470523      !                                                      ! ------------------------- ! 
    471524      !                                                      !   topmelt and botmelt     !    
    … …  
    481534         srcv(jpr_topm:jpr_botm)%laction = .TRUE. 
    482535      ENDIF 
     536      !                                                      ! ------------------------- ! 
     537      !                                                      !      Wave breaking        !     
     538      !                                                      ! ------------------------- !  
     539      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height 
     540      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN 
     541         srcv(jpr_hsig)%laction = .TRUE. 
     542         cpl_hsig = .TRUE. 
     543      ENDIF 
     544      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy 
     545      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN 
     546         srcv(jpr_phioc)%laction = .TRUE. 
     547         cpl_phioc = .TRUE. 
     548      ENDIF 
     549      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction 
     550      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN 
     551         srcv(jpr_sdrftx)%laction = .TRUE. 
     552         cpl_sdrftx = .TRUE. 
     553      ENDIF 
     554      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction 
     555      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN 
     556         srcv(jpr_sdrfty)%laction = .TRUE. 
     557         cpl_sdrfty = .TRUE. 
     558      ENDIF 
     559      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period 
     560      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN 
     561         srcv(jpr_wper)%laction = .TRUE. 
     562         cpl_wper = .TRUE. 
     563      ENDIF 
     564      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number 
     565      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN 
     566         srcv(jpr_wnum)%laction = .TRUE. 
     567         cpl_wnum = .TRUE. 
     568      ENDIF 
     569      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave 
     570      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN 
     571         srcv(jpr_wstrf)%laction = .TRUE. 
     572         cpl_wstrf = .TRUE. 
     573      ENDIF 
     574      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient 
     575      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN 
     576         srcv(jpr_wdrag)%laction = .TRUE. 
     577         cpl_wdrag = .TRUE. 
     578      ENDIF 
     579      !  
    483580      !                                                      ! ------------------------------- ! 
    484581      !                                                      !   OPA-SAS coupling - rcv by opa !    
    … …  
    635732      !                                                      ! ------------------------- ! 
    636733      ssnd(jps_fice)%clname = 'OIceFrc' 
     734      ssnd(jps_ficet)%clname = 'OIceFrcT'  
    637735      ssnd(jps_hice)%clname = 'OIceTck' 
    638736      ssnd(jps_hsnw)%clname = 'OSnwTck' 
    … …  
    643741      ENDIF 
    644742       
     743      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE.  
     744 
    645745      SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) 
    646746      CASE( 'none'         )       ! nothing to do 
    … …  
    663763      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1' 
    664764      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1' 
     765      ssnd(jps_ocxw)%clname = 'O_OCurxw'  
     766      ssnd(jps_ocyw)%clname = 'O_OCuryw'  
    665767      ! 
    666768      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold 
    … …  
    683785      END SELECT 
    684786 
     787      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold  
     788         
     789      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN  
     790         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V'  
     791      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN  
     792         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' )  
     793      ENDIF  
     794      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1.  
     795      SELECT CASE( TRIM( sn_snd_crtw%cldes ) )  
     796         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE.  
     797         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE.  
     798         CASE( 'weighted oce and ice' )   !   nothing to do  
     799         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.  
     800         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' )  
     801      END SELECT  
     802 
    685803      !                                                      ! ------------------------- ! 
    686804      !                                                      !          CO2 flux         ! 
    687805      !                                                      ! ------------------------- ! 
    688806      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE. 
     807 
     808      !                                                      ! ------------------------- !  
     809      !                                                      !     Sea surface height    !  
     810      !                                                      ! ------------------------- !  
     811      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE.  
    689812 
    690813      !                                                      ! ------------------------------- ! 
    … …  
    781904      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   & 
    782905         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 
    783       ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
     906      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    784907 
    785908      CALL wrk_dealloc( jpi,jpj,   zacs, zaos ) 
    … …  
    835958      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    836959      !!---------------------------------------------------------------------- 
    837       INTEGER, INTENT(in) ::   kt       ! ocean model time step index 
    838       INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation  
    839       INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3) 
    840  
     960      USE zdf_oce,  ONLY : ln_zdfqiao 
     961 
     962      IMPLICIT NONE 
     963 
     964      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index 
     965      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation  
     966      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3) 
    841967      !! 
    842968      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module?? 
    … …  
    9841110      ENDIF 
    9851111 
    986 #if defined key_cpl_carbon_cycle 
    9871112      !                                                      ! ================== ! 
    9881113      !                                                      ! atmosph. CO2 (ppm) ! 
    9891114      !                                                      ! ================== ! 
    9901115      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 
    991 #endif 
     1116      !  
     1117      !                                                      ! ========================= !  
     1118      !                                                      ! Mean Sea Level Pressure   !   (taum)  
     1119      !                                                      ! ========================= !  
     1120      !  
     1121      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH  
     1122          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields  
     1123 
     1124          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization  
     1125          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)  
     1126          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                         !atmospheric pressure  
     1127     
     1128          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)  
     1129      END IF  
     1130      ! 
     1131      IF( ln_sdw ) THEN  ! Stokes Drift correction activated 
     1132      !                                                      ! ========================= !  
     1133      !                                                      !       Stokes drift u      ! 
     1134      !                                                      ! ========================= !  
     1135         IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
     1136      ! 
     1137      !                                                      ! ========================= !  
     1138      !                                                      !       Stokes drift v      ! 
     1139      !                                                      ! ========================= !  
     1140         IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
     1141      ! 
     1142      !                                                      ! ========================= !  
     1143      !                                                      !      Wave mean period     ! 
     1144      !                                                      ! ========================= !  
     1145         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1) 
     1146      ! 
     1147      !                                                      ! ========================= !  
     1148      !                                                      !  Significant wave height  ! 
     1149      !                                                      ! ========================= !  
     1150         IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
     1151      ! 
     1152      !                                                      ! ========================= !  
     1153      !                                                      !    Vertical mixing Qiao   ! 
     1154      !                                                      ! ========================= !  
     1155         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1) 
     1156 
     1157         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode 
     1158         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction & 
     1159                                                                    .OR. srcv(jpr_hsig)%laction ) THEN 
     1160            CALL sbc_stokes() 
     1161            IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 
     1162         ENDIF 
     1163         IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 
     1164      ENDIF 
     1165      !                                                      ! ========================= !  
     1166      !                                                      ! Stress adsorbed by waves  ! 
     1167      !                                                      ! ========================= !  
     1168      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1) 
     1169 
     1170      !                                                      ! ========================= !  
     1171      !                                                      !   Wave drag coefficient   ! 
     1172      !                                                      ! ========================= !  
     1173      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1) 
    9921174 
    9931175      !  Fields received by SAS when OASIS coupling 
    … …  
    19192101         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info ) 
    19202102      ENDIF 
    1921       ! 
    1922 #if defined key_cpl_carbon_cycle 
    19232103      !                                                      ! ------------------------- ! 
    19242104      !                                                      !  CO2 flux from PISCES     !  
    19252105      !                                                      ! ------------------------- ! 
    1926       IF( ssnd(jps_co2)%laction )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info ) 
    1927       ! 
    1928 #endif 
     2106      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info ) 
     2107      ! 
    19292108      !                                                      ! ------------------------- ! 
    19302109      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      ! 
    … …  
    20632242      ENDIF 
    20642243      ! 
     2244      !                                                      ! ------------------------- !  
     2245      !                                                      !  Surface current to waves !  
     2246      !                                                      ! ------------------------- !  
     2247      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN  
     2248          !      
     2249          !                                                  j+1  j     -----V---F  
     2250          ! surface velocity always sent from T point                    !       |  
     2251          !                                                       j      |   T   U  
     2252          !                                                              |       |  
     2253          !                                                   j   j-1   -I-------|  
     2254          !                                               (for I)        |       |  
     2255          !                                                             i-1  i   i  
     2256          !                                                              i      i+1 (for I)  
     2257          SELECT CASE( TRIM( sn_snd_crtw%cldes ) )  
     2258          CASE( 'oce only'             )      ! C-grid ==> T  
     2259             DO jj = 2, jpjm1  
     2260                DO ji = fs_2, fs_jpim1   ! vector opt.  
     2261                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )  
     2262                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) )   
     2263                END DO  
     2264             END DO  
     2265          CASE( 'weighted oce and ice' )     
     2266             SELECT CASE ( cp_ice_msh )  
     2267             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T  
     2268                DO jj = 2, jpjm1  
     2269                   DO ji = fs_2, fs_jpim1   ! vector opt.  
     2270                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)    
     2271                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)  
     2272                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)  
     2273                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)  
     2274                   END DO  
     2275                END DO  
     2276             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T  
     2277                DO jj = 2, jpjm1  
     2278                   DO ji = 2, jpim1   ! NO vector opt.  
     2279                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)    
     2280                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)    
     2281                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &  
     2282                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2283                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &  
     2284                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2285                   END DO  
     2286                END DO  
     2287             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T  
     2288                DO jj = 2, jpjm1  
     2289                   DO ji = 2, jpim1   ! NO vector opt.  
     2290                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)    
     2291                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)    
     2292                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &  
     2293                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2294                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &  
     2295                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2296                   END DO  
     2297                END DO  
     2298             END SELECT  
     2299             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )  
     2300          CASE( 'mixed oce-ice'        )  
     2301             SELECT CASE ( cp_ice_msh )  
     2302             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T  
     2303                DO jj = 2, jpjm1  
     2304                   DO ji = fs_2, fs_jpim1   ! vector opt.  
     2305                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &  
     2306                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)  
     2307                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &  
     2308                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)  
     2309                   END DO  
     2310                END DO  
     2311             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T  
     2312                DO jj = 2, jpjm1  
     2313                   DO ji = 2, jpim1   ! NO vector opt.  
     2314                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &     
     2315                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &  
     2316                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2317                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &   
     2318                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &  
     2319                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2320                   END DO  
     2321                END DO  
     2322             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T  
     2323                DO jj = 2, jpjm1  
     2324                   DO ji = 2, jpim1   ! NO vector opt.  
     2325                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &     
     2326                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &  
     2327                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2328                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &   
     2329                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &  
     2330                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)  
     2331                   END DO  
     2332                END DO  
     2333             END SELECT  
     2334          END SELECT  
     2335         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. )  
     2336         !  
     2337         !  
     2338         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components  
     2339         !                                                                        ! Ocean component  
     2340            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component   
     2341            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component   
     2342            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components   
     2343            zoty1(:,:) = ztmp2(:,:)   
     2344            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component  
     2345               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component   
     2346               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component   
     2347               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components   
     2348               zity1(:,:) = ztmp2(:,:)  
     2349            ENDIF  
     2350         ENDIF  
     2351         !  
     2352!         ! spherical coordinates to cartesian -> 2 components to 3 components  
     2353!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN  
     2354!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents  
     2355!            ztmp2(:,:) = zoty1(:,:)  
     2356!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )  
     2357!            !  
     2358!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities  
     2359!               ztmp1(:,:) = zitx1(:,:)  
     2360!               ztmp1(:,:) = zity1(:,:)  
     2361!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )  
     2362!            ENDIF  
     2363!         ENDIF  
     2364         !  
     2365         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid  
     2366         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid  
     2367         !   
     2368      ENDIF  
     2369      !  
     2370      IF( ssnd(jps_ficet)%laction ) THEN  
     2371         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info )  
     2372      END IF  
     2373      !                                                      ! ------------------------- !  
     2374      !                                                      !   Water levels to waves   !  
     2375      !                                                      ! ------------------------- !  
     2376      IF( ssnd(jps_wlev)%laction ) THEN  
     2377         IF( ln_apr_dyn ) THEN   
     2378            IF( kt /= nit000 ) THEN   
     2379               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )   
     2380            ELSE   
     2381               ztmp1(:,:) = sshb(:,:)   
     2382            ENDIF   
     2383         ELSE   
     2384            ztmp1(:,:) = sshn(:,:)   
     2385         ENDIF   
     2386         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )  
     2387      END IF  
    20652388      ! 
    20662389      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r6140 r7403  
    9090      INTEGER, INTENT( in ) :: kt                   ! ocean time step 
    9191      ! 
    92       INTEGER               :: ji, jj               ! loop index 
     92      INTEGER               :: ji, jj, jk           ! loop index 
     93      INTEGER               :: ikt, ikb             ! loop index 
    9394      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)  
     95      REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 
     96      REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d 
    9497      !!--------------------------------------------------------------------- 
    9598      ! 
    … …  
    161164         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 
    162165         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 
    163          CALL lbc_lnk(fwfisf(:,:)         ,'T',1.) 
    164          CALL lbc_lnk(qisf(:,:)           ,'T',1.) 
     166         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.) 
     167         CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
     168 
     169!============================================================================================================================================= 
     170         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 
     171            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
     172            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        ) 
     173 
     174            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s) 
     175            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2) 
     176            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2) 
     177            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2) 
     178 
     179            DO jj = 1,jpj 
     180               DO ji = 1,jpi 
     181                  ikt = misfkt(ji,jj) 
     182                  ikb = misfkb(ji,jj) 
     183                  DO jk = ikt, ikb - 1 
     184                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
     185                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
     186                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
     187                  END DO 
     188                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
     189                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
     190                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
     191               END DO 
     192            END DO 
     193 
     194            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 
     195            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 
     196            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 
     197            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  )) 
     198 
     199            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
     200            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        ) 
     201         END IF 
     202 
     203         ! output 
     204         CALL iom_put('qlatisf'  , qisf) 
     205         CALL iom_put('fwfisf', fwfisf) 
     206!============================================================================================================================================= 
    165207 
    166208         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
    … …  
    177219         END IF 
    178220         !  
    179          ! output 
    180          CALL iom_put('qisf'  , qisf) 
    181          CALL iom_put('fwfisf', fwfisf) 
    182  
    183221         ! deallocation 
    184222         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  ) 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r6460 r7403  
    8989      NAMELIST/namsbc/ nn_fsbc  , ln_ana   , ln_flx, ln_blk_clio, ln_blk_core, ln_blk_mfs,   & 
    9090         &             ln_cpl   , ln_mixcpl, nn_components      , nn_limflx  ,               & 
    91          &             ln_traqsr, ln_dm2dc ,                                                 &   
     91         &             ln_traqsr, ln_dm2dc ,                                                 & 
    9292         &             nn_ice   , nn_ice_embd,                                               & 
    9393         &             ln_rnf   , ln_ssr   , ln_isf   , nn_fwb    , ln_apr_dyn,              & 
    94          &             ln_wave  ,                                                            & 
    95          &             nn_lsm    
     94         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauoc  , ln_stcor  ,              & 
     95         &             nn_lsm 
    9696      INTEGER  ::   ios 
    9797      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm 
    … …  
    153153         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea 
    154154         WRITE(numout,*) '              nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm 
    155          WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave   
     155         WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave 
     156         WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw        = ', ln_sdw 
     157         WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc      = ', ln_tauoc 
     158         WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor      = ', ln_stcor 
     159         WRITE(numout,*) '                 neutral drag coefficient (CORE, MFS)    ln_cdgw       = ', ln_cdgw 
    156160      ENDIF 
    157161      ! 
    … …  
    220224         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 
    221225       
     226      IF ( ln_wave ) THEN 
     227      !Activated wave module but neither drag nor stokes drift activated 
     228         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN 
     229            CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F') 
     230      !drag coefficient read from wave model definable only with mfs bulk formulae and core  
     231         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
     232             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
     233         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN 
     234             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     235         ENDIF 
     236      ELSE 
     237      IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor )                &  
     238         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
     239         &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
     240         &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
     241         &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      &   
     242         &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
     243      ENDIF  
    222244      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
    223245      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl 
    … …  
    357379            &                CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS 
    358380      END SELECT 
    359  
     381      IF ( ln_wave .AND. ln_tauoc) THEN                                 ! Wave stress subctracted 
     382            utau(:,:) = utau(:,:)*tauoc_wave(:,:) 
     383            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 
     384            taum(:,:) = taum(:,:)*tauoc_wave(:,:) 
     385      ! 
     386            SELECT CASE( nsbc ) 
     387            CASE(  0,1,2,3,5,-1 )  ; 
     388                IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. & 
     389                        & If not requested select ln_tauoc=.false' 
     390            END SELECT 
     391      ! 
     392      END IF 
    360393      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    361394 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r6140 r7403  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient  
    7    !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift  
    8    !!---------------------------------------------------------------------- 
    9  
    10    !!---------------------------------------------------------------------- 
    11    !!   sbc_wave      : read drag coefficient from wave model in netcdf files  
     6   !! History :  3.3  !   2011-09  (M. Adani)  Original code: Drag Coefficient  
     7   !!         :  3.4  !   2012-10  (M. Adani)  Stokes Drift  
     8   !!            3.6  !   2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation 
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   sbc_wave      : wave data from wave model in netcdf files  
    1213   !!---------------------------------------------------------------------- 
    1314   USE oce            !  
    14    USE sbc_oce        ! Surface boundary condition: ocean fields 
     15   USE sbc_oce       ! Surface boundary condition: ocean fields 
    1516   USE bdy_oce        ! 
    1617   USE domvvl         ! 
    17    ! 
    1818   USE iom            ! I/O manager library 
    1919   USE in_out_manager ! I/O manager 
    2020   USE lib_mpp        ! distribued memory computing library 
    21    USE fldread        ! read input fields 
     21   USE fldread       ! read input fields 
    2222   USE wrk_nemo       ! 
     23   USE phycst         ! physical constants  
    2324 
    2425   IMPLICIT NONE 
    2526   PRIVATE 
    2627 
    27    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
     28   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
     29   PUBLIC   sbc_wave    ! routine called in sbcmod 
    2830    
    29    INTEGER , PARAMETER ::   jpfld  = 3   ! maximum number of files to read for srokes drift 
    30    INTEGER , PARAMETER ::   jp_usd = 1   ! index of stokes drift  (i-component) (m/s)    at T-point 
    31    INTEGER , PARAMETER ::   jp_vsd = 2   ! index of stokes drift  (j-component) (m/s)    at T-point 
    32    INTEGER , PARAMETER ::   jp_wn  = 3   ! index of wave number                 (1/m)    at T-point 
    33  
    34    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    35    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    36  
    37    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:)   :: cdn_wave  
    38    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:,:) :: usd3d, vsd3d, wsd3d  
    39    REAL(wp),         ALLOCATABLE, DIMENSION (:,:)   :: usd2d, vsd2d, uwavenum, vwavenum  
     31   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
     32   LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE. 
     33   LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE. 
     34   LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE. 
     35   LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE. 
     36   LOGICAL, PUBLIC     ::   cpl_wper=.FALSE. 
     37   LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE. 
     38   LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE. 
     39   LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE. 
     40 
     41   INTEGER ::   jpfld                ! number of files to read for stokes drift 
     42   INTEGER ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point 
     43   INTEGER ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point 
     44   INTEGER ::   jp_swh               ! index of significant wave hight      (m)      at T-point 
     45   INTEGER ::   jp_wmp               ! index of mean wave period            (s)      at T-point 
     46 
     47   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     48   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
     49   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
     50   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     51   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave  
     52   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum 
     53   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave 
     54   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d 
     55   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt 
     56   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d  
     57   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt 
    4058 
    4159   !! * Substitutions 
    … …  
    4866CONTAINS 
    4967 
     68   SUBROUTINE sbc_stokes( ) 
     69      !!--------------------------------------------------------------------- 
     70      !!                     ***  ROUTINE sbc_stokes  *** 
     71      !! 
     72      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al., 
     73      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1) 
     74      !! 
     75      !! ** Method  : - Calculate Stokes transport speed  
     76      !!              - Calculate horizontal divergence  
     77      !!              - Integrate the horizontal divergenze from the bottom  
     78      !! ** action   
     79      !!--------------------------------------------------------------------- 
     80      INTEGER                ::   jj,ji,jk  
     81      REAL(wp)                       ::  ztransp, zfac, zsp0, zk, zus, zvs 
     82      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
     83      !!--------------------------------------------------------------------- 
     84      ! 
     85 
     86      CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
     87      DO jk = 1, jpk 
     88         DO jj = 1, jpj 
     89            DO ji = 1, jpi 
     90               ! On T grid 
     91               ! Stokes transport speed estimated from Hs and Tmean 
     92               ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
     93               ! Stokes surface speed 
     94               zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2) 
     95               ! Wavenumber scale 
     96               zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
     97               ! Depth attenuation 
     98               zfac = EXP(-2.0_wp*zk*gdept_n(ji,jj,jk))/(1.0_wp+8.0_wp*zk*gdept_n(ji,jj,jk)) 
     99               ! 
     100               usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
     101               vsd3dt(ji,jj,jk) = zfac * zvsd2dt(ji,jj) * tmask(ji,jj,jk) 
     102            END DO 
     103         END DO 
     104      END DO  
     105      ! Into the U and V Grid 
     106      DO jk = 1, jpkm1 
     107         DO jj = 1, jpjm1 
     108            DO ji = 1, fs_jpim1 
     109               usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   & 
     110                               &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) ) 
     111               vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   & 
     112                               &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) ) 
     113            END DO 
     114         END DO 
     115      END DO 
     116      ! 
     117      CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
     118      CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
     119      ! 
     120      DO jk = 1, jpkm1               ! Horizontal divergence 
     121         DO jj = 2, jpj 
     122            DO ji = fs_2, jpi 
     123               ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     & 
     124                  &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     & 
     125                  &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     & 
     126                  &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj) 
     127            END DO 
     128         END DO 
     129      END DO 
     130      ! 
     131      IF( .NOT. AGRIF_Root() ) THEN 
     132         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east 
     133         IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west 
     134         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north 
     135         IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south 
     136      ENDIF 
     137      ! 
     138      CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
     139      ! 
     140      DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence 
     141         wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - e3t_n(:,:,jk) * ze3hdiv(:,:,jk) 
     142      END DO 
     143#if defined key_bdy 
     144      IF( lk_bdy ) THEN 
     145         DO jk = 1, jpkm1 
     146            wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
     147         END DO 
     148      ENDIF 
     149#endif 
     150      CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 
     151      ! 
     152   END SUBROUTINE sbc_stokes 
     153 
     154   SUBROUTINE sbc_qiao 
     155      !!--------------------------------------------------------------------- 
     156      !!                     ***  ROUTINE sbc_qiao  *** 
     157      !! 
     158      !! ** Purpose :   Qiao formulation for wave enhanced turbulence 
     159      !!                2010 (DOI: 10.1007/s10236-010-0326)  
     160      !! 
     161      !! ** Method  : -  
     162      !! ** action   
     163      !!--------------------------------------------------------------------- 
     164      INTEGER :: jj, ji 
     165 
     166      ! Calculate the module of the stokes drift on T grid 
     167      !------------------------------------------------- 
     168      DO jj = 1, jpj 
     169         DO ji = 1, jpi 
     170            tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj) * zusd2dt(ji,jj) + zvsd2dt(ji,jj) * zvsd2dt(ji,jj) ) 
     171         END DO 
     172      END DO 
     173      ! 
     174   END SUBROUTINE sbc_qiao 
     175 
    50176   SUBROUTINE sbc_wave( kt ) 
    51177      !!--------------------------------------------------------------------- 
    52       !!                     ***  ROUTINE sbc_apr  *** 
    53       !! 
    54       !! ** Purpose :   read drag coefficient from wave model  in netcdf files. 
     178      !!                     ***  ROUTINE sbc_wave  *** 
     179      !! 
     180      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
    55181      !! 
    56182      !! ** Method  : - Read namelist namsbc_wave 
    57183      !!              - Read Cd_n10 fields in netcdf files  
    58184      !!              - Read stokes drift 2d in netcdf files  
    59       !!              - Read wave number      in netcdf files  
    60       !!              - Compute 3d stokes drift using monochromatic 
    61       !! ** action  :    
    62       !!--------------------------------------------------------------------- 
    63       INTEGER, INTENT( in  ) ::   kt       ! ocean time step 
     185      !!              - Read wave number in netcdf files  
     186      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     187      !!                formulation 
     188      !! ** action   
     189      !!--------------------------------------------------------------------- 
     190      USE zdf_oce,  ONLY : ln_zdfqiao 
     191 
     192      INTEGER, INTENT( in  ) :: kt       ! ocean time step 
    64193      ! 
    65194      INTEGER                ::   ierror   ! return error code 
    66       INTEGER                ::   ifpr, jj,ji,jk  
    67       INTEGER                ::   ios     ! Local integer output status for namelist read 
    68       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     195      INTEGER                ::   ifpr 
     196      INTEGER                ::   ios      ! Local integer output status for namelist read 
     197      ! 
    69198      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    70       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read 
    71       REAL(wp), DIMENSION(:,:,:), POINTER ::   zusd_t, zvsd_t, ze3hdiv   ! 3D workspace 
    72       !! 
    73       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn, ln_cdgw , ln_sdw 
     199      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
     200      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
     201                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     202      !! 
     203      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc 
    74204      !!--------------------------------------------------------------------- 
    75205      ! 
    … …  
    80210         READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    81211901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    82          ! 
     212          
    83213         REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    84214         READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    … …  
    86216         IF(lwm) WRITE ( numond, namsbc_wave ) 
    87217         ! 
    88          IF(lwp) THEN               ! Control print 
    89             WRITE(numout,*) '        Namelist namsbc_wave : surface wave setting'  
    90             WRITE(numout,*) '           wave drag coefficient                      ln_cdgw  = ', ln_cdgw   
    91             WRITE(numout,*) '           wave stokes drift                          ln_sdw   = ', ln_sdw 
    92          ENDIF 
    93          ! 
    94          IF( .NOT.( ln_cdgw .OR. ln_sdw ) )    & 
    95             &  CALL ctl_warn( 'ln_sbcwave=T but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
    96          IF( ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )   &        
    97             &  CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
    98          ! 
    99218         IF( ln_cdgw ) THEN 
    100             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    101             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    102             ! 
    103                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    104             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    105             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     219            IF( .NOT. cpl_wdrag ) THEN 
     220               ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     221               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     222               ! 
     223                                      ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
     224               IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
     225               CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     226            ENDIF 
    106227            ALLOCATE( cdn_wave(jpi,jpj) ) 
    107             cdn_wave(:,:) = 0.0 
    108          ENDIF 
     228         ENDIF 
     229 
     230         IF( ln_tauoc ) THEN 
     231            IF( .NOT. cpl_wstrf ) THEN 
     232               ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     233               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     234               ! 
     235                                       ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     236               IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     237               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     238            ENDIF 
     239            ALLOCATE( tauoc_wave(jpi,jpj) ) 
     240         ENDIF 
     241 
    109242         IF( ln_sdw ) THEN 
    110             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 
    111             ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    112             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    113             ! 
    114             DO ifpr= 1, jpfld 
    115                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    116                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    117             END DO 
    118             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    119             ALLOCATE( usd2d(jpi,jpj) , vsd2d(jpi,jpj) , uwavenum(jpi,jpj) , vwavenum(jpi,jpj) ) 
     243            ! Find out how many fields have to be read from file if not coupled 
     244            jpfld=0 
     245            jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 
     246            IF( .NOT. cpl_sdrftx ) THEN 
     247               jpfld=jpfld+1 
     248               jp_usd=jpfld 
     249            ENDIF 
     250            IF( .NOT. cpl_sdrfty ) THEN 
     251               jpfld=jpfld+1 
     252               jp_vsd=jpfld 
     253            ENDIF 
     254            IF( .NOT. cpl_hsig ) THEN 
     255               jpfld=jpfld+1 
     256               jp_swh=jpfld 
     257            ENDIF 
     258            IF( .NOT. cpl_wper ) THEN 
     259               jpfld=jpfld+1 
     260               jp_wmp=jpfld 
     261            ENDIF 
     262 
     263            ! Read from file only the non-coupled fields  
     264            IF( jpfld > 0 ) THEN 
     265               ALLOCATE( slf_i(jpfld) ) 
     266               IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 
     267               IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 
     268               IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 
     269               IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 
     270               ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
     271               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     272               ! 
     273               DO ifpr= 1, jpfld 
     274                  ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     275                  IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     276               END DO 
     277 
     278               CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     279            ENDIF 
    120280            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    121             usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;    uwavenum(:,:) = 0._wp 
    122             vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ;    vwavenum(:,:) = 0._wp 
     281            ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) ) 
     282            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 
     283            ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) ) 
     284            usd3d(:,:,:) = 0._wp 
     285            vsd3d(:,:,:) = 0._wp 
    123286            wsd3d(:,:,:) = 0._wp 
    124          ENDIF 
    125       ENDIF 
    126       ! 
    127       IF( ln_cdgw ) THEN               !==  Neutral drag coefficient  ==! 
     287            IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
     288               IF( .NOT. cpl_wnum ) THEN 
     289                  ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     290                  IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
     291                                         ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     292                  IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     293                  CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     294               ENDIF 
     295               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
     296            ENDIF 
     297         ENDIF 
     298      ENDIF 
     299      ! 
     300      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
    128301         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    129302         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    130303      ENDIF 
    131       ! 
    132       IF( ln_sdw )  THEN               !==  Computation of the 3d Stokes Drift  ==! 
     304 
     305      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
     306         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
     307         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     308      ENDIF 
     309 
     310      IF( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    133311         ! 
    134          CALL wrk_alloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv ) 
     312         ! Read from file only if the field is not coupled 
     313         IF( jpfld > 0 ) THEN 
     314            CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
     315            IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
     316            IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     317            IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     318            IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     319         ENDIF 
    135320         ! 
    136          CALL fld_read( kt, nn_fsbc, sf_sd )    !* read drag coefficient from external forcing 
     321         ! Read also wave number if needed, so that it is available in coupling routines 
     322         IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
     323            CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing 
     324            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     325         ENDIF 
     326            
     327         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
     328         !(DOI: 10.1175/JPO-D-14-0020.1)==!  
    137329         ! 
    138          DO jk = 1, jpkm1                       !* distribute it on the vertical 
    139             zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) ) 
    140             zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) ) 
    141          END DO 
    142          DO jk = 1, jpkm1                       !* interpolate the stokes drift from t-point to u- and v-points 
    143             DO jj = 1, jpjm1 
    144                DO ji = 1, jpim1 
    145                    usd3d(ji,jj,jk) = 0.5_wp * ( zusd_t(ji  ,jj,jk) + zusd_t(ji+1,jj,jk) ) * umask(ji,jj,jk) 
    146                    vsd3d(ji,jj,jk) = 0.5_wp * ( zvsd_t(ji  ,jj,jk) + zvsd_t(ji,jj+1,jk) ) * vmask(ji,jj,jk) 
    147                END DO 
    148             END DO 
    149          END DO 
    150          CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
    151          CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
    152          ! 
    153          DO jk = 1, jpkm1                       !* e3t * Horizontal divergence  ==! 
    154             DO jj = 2, jpjm1 
    155                DO ji = fs_2, fs_jpim1   ! vector opt. 
    156                   ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     & 
    157                      &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     & 
    158                      &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     & 
    159                      &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj) 
    160                END DO   
    161             END DO   
    162             IF( .NOT. AGRIF_Root() ) THEN 
    163                IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,jk) = 0._wp      ! east 
    164                IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,jk) = 0._wp      ! west 
    165                IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,jk) = 0._wp      ! north 
    166                IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,jk) = 0._wp      ! south 
    167             ENDIF 
    168          END DO 
    169          CALL lbc_lnk( ze3hdiv, 'T', 1. )  
    170          ! 
    171          DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence 
    172             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 
    173          END DO 
    174 #if defined key_bdy 
    175          IF( lk_bdy ) THEN 
    176             DO jk = 1, jpkm1 
    177                wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    178             END DO 
    179          ENDIF 
    180 #endif 
    181          CALL wrk_dealloc( jpi,jpj,jpk,   zusd_t, zvsd_t, ze3hdiv ) 
    182          !  
     330         ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 
     331         IF( jpfld == 4 ) THEN 
     332            CALL sbc_stokes() 
     333            IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
     334               CALL sbc_qiao() 
     335            ENDIF 
     336         ENDIF 
    183337      ENDIF 
    184338      ! 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r6140 r7403  
    99   !!            3.7  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes  
    1010   !!             -   !  2014-12  (G. Madec) suppression of cross land advection option 
     11   !!            3.6  !  2015-06  (E. Clementi) Addition of Stokes drift in case of wave coupling 
    1112   !!---------------------------------------------------------------------- 
    1213 
    … …  
    2627   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff. 
    2728   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces 
     29   USE trd_oce         ! trends: ocean variables 
     30   USE trdtra          ! trends manager: tracers  
    2831   ! 
    2932   USE in_out_manager ! I/O manager 
    … …  
    3336   USE wrk_nemo       ! Memory Allocation 
    3437   USE timing         ! Timing 
    35  
    36    USE diaptr          ! Poleward heat transport  
     38   USE sbcwave        ! wave module 
     39   USE sbc_oce        ! surface boundary condition: ocean 
     40   USE diaptr         ! Poleward heat transport  
    3741 
    3842   IMPLICIT NONE 
    … …  
    8690      INTEGER ::   jk   ! dummy loop index 
    8791      REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn 
     92      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace 
    8893      !!---------------------------------------------------------------------- 
    8994      ! 
    … …  
    9398      ! 
    9499      !                                          ! set time step 
     100      zun(:,:,:) = 0.0 
     101      zvn(:,:,:) = 0.0 
     102      zwn(:,:,:) = 0.0 
     103      !     
    95104      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    96105         r2dt = rdt                                 ! = rdt (restarting with Euler time stepping) 
    … …  
    100109      ! 
    101110      !                                         !==  effective transport  ==! 
    102       DO jk = 1, jpkm1 
    103          zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)                  ! eulerian transport only 
    104          zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    105          zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk) 
    106       END DO 
     111      IF( ln_wave .AND. ln_sdw )  THEN 
     112         DO jk = 1, jpkm1 
     113            zun(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) *      & 
     114                        &  ( un(:,:,jk) + usd3d(:,:,jk) )                       ! eulerian transport + Stokes Drift 
     115            zvn(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) *      & 
     116                        &  ( vn(:,:,jk) + vsd3d(:,:,jk) ) 
     117            zwn(:,:,jk) = e1e2t(:,:) *                    & 
     118                        &  ( wn(:,:,jk) + wsd3d(:,:,jk) ) 
     119         END DO 
     120      ELSE 
     121         DO jk = 1, jpkm1 
     122            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * un(:,:,jk)               ! eulerian transport only 
     123            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     124            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk) 
     125         END DO 
     126      ENDIF 
    107127      ! 
    108128      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections 
    … …  
    127147      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF  
    128148!!gm ??? 
     149      ! 
     150      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     151         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     152         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     153         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     154      ENDIF 
    129155      ! 
    130156      SELECT CASE ( nadv )                      !==  compute advection trend and add it to general trend  ==! 
    … …  
    145171      END SELECT 
    146172      ! 
    147       !                                         ! print mean trends (used for debugging) 
     173      IF( l_trdtra )   THEN                      ! save the advective trends for further diagnostics 
     174         DO jk = 1, jpkm1 
     175            ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk) 
     176            ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk) 
     177         END DO 
     178         CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt ) 
     179         CALL trd_tra( kt, 'TRA', jp_sal, jptra_totad, ztrds ) 
     180         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     181      ENDIF 
     182      !                                              ! print mean trends (used for debugging) 
    148183      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv  - Ta: ', mask1=tmask,               & 
    149184         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90

    r6140 r7403  
    1818   USE trdtra         ! trends manager: tracers  
    1919   USE diaptr         ! poleward transport diagnostics 
     20   USE diaar5         ! AR5 diagnostics 
    2021   ! 
    2122   USE in_out_manager ! I/O manager 
    … …  
    3233    
    3334   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6 
     35 
     36   LOGICAL :: l_trd   ! flag to compute trends 
     37   LOGICAL :: l_ptr   ! flag to compute poleward transport 
     38   LOGICAL :: l_hst   ! flag to compute heat/salt transport 
    3439 
    3540   !! * Substitutions 
    … …  
    8893      ENDIF 
    8994      ! 
     95      l_trd = .FALSE. 
     96      l_hst = .FALSE. 
     97      l_ptr = .FALSE. 
     98      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )        l_trd = .TRUE. 
     99      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE.  
     100      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     101         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
     102      ! 
    90103      !                     
    91104      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers 
    … …  
    184197         END DO 
    185198         !                             ! trend diagnostics 
    186          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN 
     199         IF( l_trd ) THEN 
    187200            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
    188201            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    189202            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 
    190203         END IF 
    191          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    192          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    193            IF( jn == jp_tem )   htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    194            IF( jn == jp_sal )   str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    195          ENDIF 
     204         !                                 ! "Poleward" heat and salt transports  
     205         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
     206         !                                 !  heat and salt transport 
     207         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) ) 
    196208         ! 
    197209      END DO 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r6771 r7403  
    2020   USE trdtra         ! tracers trends 
    2121   USE diaptr         ! poleward transport diagnostics 
     22   USE diaar5         ! AR5 diagnostics 
     23   USE phycst, ONLY: rau0_rcp 
    2224   ! 
    2325   USE in_out_manager ! I/O manager 
     26   USE iom 
    2427   USE lib_mpp        ! MPP library 
    2528   USE lbclnk         ! ocean lateral boundary condition (or mpp link)  
    … …  
    3639 
    3740   LOGICAL  ::   l_trd   ! flag to compute trends 
     41   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
     42   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport 
    3843   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6 
    3944 
    … …  
    8085      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      - 
    8186      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
    82       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz 
     87      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz, zptry 
     88      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d 
    8389      !!---------------------------------------------------------------------- 
    8490      ! 
    … …  
    94100      ! 
    95101      l_trd = .FALSE. 
    96       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    97       ! 
    98       IF( l_trd )  THEN 
     102      l_hst = .FALSE. 
     103      l_ptr = .FALSE. 
     104      IF( ( cdtype == 'TRA'   .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )     l_trd = .TRUE. 
     105      IF(   cdtype == 'TRA'   .AND. ln_diaptr )                                              l_ptr = .TRUE.  
     106      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     107         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
     108      ! 
     109      IF( l_trd .OR. l_hst )  THEN 
    99110         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
    100111         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp 
    101112      ENDIF 
    102113      ! 
     114      IF( l_ptr ) THEN   
     115         CALL wrk_alloc( jpi, jpj, jpk, zptry ) 
     116         zptry(:,:,:) = 0._wp 
     117      ENDIF 
    103118      !                          ! surface & bottom value : flux set to zero one for all 
    104119      zwz(:,:, 1 ) = 0._wp             
    … …  
    161176         CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign) 
    162177         !                 
    163          IF( l_trd )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
     178         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    164179            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
    165180         END IF 
    166181         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    167          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    168            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    169            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    170          ENDIF 
     182         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:)  
    171183         ! 
    172184         !        !==  anti-diffusive flux : high order minus low order  ==! 
    … …  
    292304         END DO 
    293305         ! 
    294          IF( l_trd ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
     306         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
    295307            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    296308            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    297309            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    298             ! 
     310         ENDIF 
     311            ! 
     312         IF( l_trd ) THEN  
    299313            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
    300314            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
    301315            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
    302316            ! 
    303             CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    304317         END IF 
    305          !                    ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    306          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    307            IF( jn == jp_tem )   htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) ) 
    308            IF( jn == jp_sal )   str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) ) 
     318         !                                !  heat/salt transport 
     319         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 
     320 
     321         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     322         IF( l_ptr ) THEN   
     323            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
     324            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 
    309325         ENDIF 
    310326         ! 
    311327      END DO                     ! end of tracer loop 
    312328      ! 
    313       CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
     329                              CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 
     330      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
     331      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry ) 
    314332      ! 
    315333      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct') 
    … …  
    357375      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav 
    358376      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz 
     377      REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry 
    359378      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs 
    360379      !!---------------------------------------------------------------------- 
    … …  
    373392      ! 
    374393      l_trd = .FALSE. 
    375       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
    376       ! 
    377       IF( l_trd )  THEN 
     394      l_hst = .FALSE. 
     395      l_ptr = .FALSE. 
     396      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
     397      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
     398      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     399         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
     400      ! 
     401      IF( l_trd .OR. l_hst )  THEN 
    378402         CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
    379403         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp 
    380404      ENDIF 
    381405      ! 
     406      IF( l_ptr ) THEN   
     407         CALL wrk_alloc( jpi, jpj,jpk, zptry ) 
     408         zptry(:,:,:) = 0._wp 
     409      ENDIF 
    382410      zwi(:,:,:) = 0._wp 
    383411      z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 
    … …  
    445473         CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign) 
    446474         !                 
    447          IF( l_trd )  THEN                ! trend diagnostics (contribution of upstream fluxes) 
     475         IF( l_trd .OR. l_hst )  THEN                ! trend diagnostics (contribution of upstream fluxes) 
    448476            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:) 
    449477         END IF 
    450478         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    451          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    452            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    453            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    454          ENDIF 
     479         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:) 
    455480 
    456481         ! 3. anti-diffusive flux : high order minus low order 
    … …  
    568593         END DO 
    569594 
    570          !                                 ! trend diagnostics (contribution of upstream fluxes) 
    571          IF( l_trd )  THEN  
     595        ! 
     596         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes) 
    572597            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed 
    573598            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
    574599            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed 
    575             ! 
    576             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )    
    577             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )   
    578             CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )  
    579             ! 
    580             CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz ) 
     600         ENDIF 
     601            ! 
     602         IF( l_trd ) THEN  
     603            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
     604            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
     605            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
     606            ! 
    581607         END IF 
    582          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    583          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    584            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 
    585            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 
     608         !                                             ! heat/salt transport 
     609         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 
     610 
     611         !                                            ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     612         IF( l_ptr ) THEN   
     613            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed 
     614            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 
    586615         ENDIF 
    587616         ! 
    588617      END DO 
    589618      ! 
    590       CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
    591       CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
    592       CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
     619                              CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav ) 
     620                              CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 
     621                              CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs ) 
     622      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 
     623      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry ) 
    593624      ! 
    594625      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts') 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

    r6140 r7403  
    2323   USE sbcrnf         ! river runoffs 
    2424   USE diaptr         ! poleward transport diagnostics 
     25   USE diaar5         ! AR5 diagnostics 
     26 
    2527   ! 
     28   USE iom 
    2629   USE wrk_nemo       ! Memory Allocation 
    2730   USE timing         ! Timing 
    … …  
    4043   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index 
    4144    
     45   LOGICAL  ::   l_trd   ! flag to compute trends 
     46   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
     47   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport 
     48 
    4249   !! * Substitutions 
    4350#  include "vectopt_loop_substitute.h90" 
    … …  
    116123      ENDIF  
    117124      !       
     125      l_trd = .FALSE. 
     126      l_hst = .FALSE. 
     127      l_ptr = .FALSE. 
     128      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
     129      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
     130      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     131         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
     132      ! 
    118133      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
    119134         ! 
    … …  
    192207         END DO         
    193208         !                                ! trend diagnostics 
    194          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   & 
    195             &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN 
     209         IF( l_trd )  THEN 
    196210            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 
    197211            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 
    198212         END IF 
    199          !                                ! "Poleward" heat and salt transports 
    200          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    201             IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    202             IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    203          ENDIF 
     213         !                                 ! "Poleward" heat and salt transports  
     214         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
     215         !                                 !  heat transport 
     216         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) ) 
    204217         ! 
    205218         !                          !* Vertical advective fluxes 
    … …  
    262275         END DO 
    263276         !                                ! send trends for diagnostic 
    264          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     & 
    265             &( cdtype == 'TRC' .AND. l_trdtrc )      )   & 
    266             CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
     277         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 
    267278         ! 
    268279      END DO                     ! end of tracer loop 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r6140 r7403  
    3434   PUBLIC   tra_adv_qck   ! routine called by step.F90 
    3535 
    36    LOGICAL  :: l_trd           ! flag to compute trends 
    3736   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio 
     37 
     38   LOGICAL  ::   l_trd   ! flag to compute trends 
     39   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
     40 
    3841 
    3942   !! * Substitutions 
    … …  
    103106      ! 
    104107      l_trd = .FALSE. 
    105       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE. 
     108      l_ptr = .FALSE. 
     109      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
     110      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
     111      ! 
    106112      ! 
    107113      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 
    … …  
    224230         END DO 
    225231         !                                 ! trend diagnostics 
    226          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
     232         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 
    227233         ! 
    228234      END DO 
    … …  
    347353         END DO 
    348354         !                                 ! trend diagnostics 
    349          IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
     355         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    350356         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    351          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    352            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
    353            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    354          ENDIF 
     357         IF( l_ptr )                     CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 
    355358         ! 
    356359      END DO 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r6140 r7403  
    1919   USE trdtra         ! trends manager: tracers  
    2020   USE diaptr         ! poleward transport diagnostics 
     21   USE diaar5         ! AR5 diagnostics 
     22 
    2123   ! 
     24   USE iom 
    2225   USE lib_mpp        ! I/O library 
    2326   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    … …  
    3235   PUBLIC   tra_adv_ubs   ! routine called by traadv module 
    3336 
    34    LOGICAL :: l_trd  ! flag to compute trends or not 
     37   LOGICAL :: l_trd   ! flag to compute trends 
     38   LOGICAL :: l_ptr   ! flag to compute poleward transport 
     39   LOGICAL :: l_hst   ! flag to compute heat transport 
     40 
    3541 
    3642   !! * Substitutions 
    … …  
    109115      ! 
    110116      l_trd = .FALSE. 
    111       IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
     117      l_hst = .FALSE. 
     118      l_ptr = .FALSE. 
     119      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
     120      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE.  
     121      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     122         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE. 
    112123      ! 
    113124      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers 
    … …  
    176187             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) ) 
    177188         END IF 
    178          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    179          IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
    180             IF( jn == jp_tem )  htr_adv(:) = ptr_sj( ztv(:,:,:) ) 
    181             IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) ) 
    182          ENDIF 
     189         !      
     190         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     191         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) ) 
     192         !                                !  heati/salt transport 
     193         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) ) 
     194         ! 
    183195         ! 
    184196         !                       !== vertical advective trend  ==! 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r6140 r7403  
    2424   USE ldfslp         ! iso-neutral slopes 
    2525   USE diaptr         ! poleward transport diagnostics 
     26   USE diaar5         ! AR5 diagnostics 
    2627   ! 
    2728   USE in_out_manager ! I/O manager 
    … …  
    3637 
    3738   PUBLIC   tra_ldf_iso   ! routine called by step.F90 
     39 
     40   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
     41   LOGICAL  ::   l_hst   ! flag to compute heat transport 
    3842 
    3943   !! * Substitutions 
    … …  
    107111      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      - 
    108112      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      - 
    109 #if defined key_diaar5 
    110       REAL(wp) ::   zztmp   ! local scalar 
    111 #endif 
    112113      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d 
    113114      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw  
    … …  
    127128         ah_wslp2(:,:,:) = 0._wp 
    128129      ENDIF 
    129       !                                               ! set time step size (Euler/Leapfrog) 
     130      !    
     131      l_hst = .FALSE. 
     132      l_ptr = .FALSE. 
     133      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE.  
     134      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     135         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
     136      ! 
     137      !                                            ! set time step size (Euler/Leapfrog) 
    130138      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    131139      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
    … …  
    369377            ! 
    370378            !                             ! "Poleward" diffusive heat or salt transports (T-S case only) 
    371             IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    372                ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    373                IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
    374                IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
    375             ENDIF 
    376             ! 
    377             IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
    378               ! 
    379               IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    380                   z2d(:,:) = zftu(ji,jj,1)  
    381                   DO jk = 2, jpkm1 
    382                      DO jj = 2, jpjm1 
    383                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    384                            z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
    385                         END DO 
    386                      END DO 
    387                   END DO 
    388 !!gm CAUTION I think there is an error of sign when using BLP operator.... 
    389 !!gm         a multiplication by zsign is required (to be checked twice !) 
    390                   z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    391                   CALL lbc_lnk( z2d, 'U', -1. ) 
    392                   CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    393                   ! 
    394                   z2d(:,:) = zftv(ji,jj,1)  
    395                   DO jk = 2, jpkm1 
    396                      DO jj = 2, jpjm1 
    397                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    398                            z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
    399                         END DO 
    400                      END DO 
    401                   END DO 
    402                   z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    403                   CALL lbc_lnk( z2d, 'V', -1. ) 
    404                   CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
    405                END IF 
    406                ! 
    407             ENDIF 
     379               ! note sign is reversed to give down-gradient diffusive transports ) 
     380            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:)  ) 
     381            !                          ! Diffusive heat transports 
     382            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) ) 
    408383            ! 
    409384         ENDIF                                                    !== end pass selection  ==! 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90

    r6140 r7403  
    1717   USE traldf_triad   ! iso-neutral lateral diffusion (triad    operator)     (tra_ldf_triad routine) 
    1818   USE diaptr         ! poleward transport diagnostics 
     19   USE diaar5         ! AR5 diagnostics 
    1920   USE trc_oce        ! share passive tracers/Ocean variables 
    2021   USE zpshde         ! partial step: hor. derivative     (zps_hde routine) 
    … …  
    2526   USE timing         ! Timing 
    2627   USE wrk_nemo       ! Memory allocation 
     28   USE iom 
    2729 
    2830   IMPLICIT NONE 
    … …  
    3941   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator 
    4042   INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator 
     43 
     44   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
     45   LOGICAL  ::   l_hst   ! flag to compute heat transport 
    4146 
    4247   !! * Substitutions 
    … …  
    95100      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev )  
    96101      ! 
     102      l_hst = .FALSE. 
     103      l_ptr = .FALSE. 
     104      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE.  
     105      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     106         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
     107      ! 
    97108      !                                !==  Initialization of metric arrays used for all tracers  ==! 
    98109      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    … …  
    150161         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==! 
    151162             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==! 
    152             IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    153                IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( -ztv(:,:,:) ) 
    154                IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( -ztv(:,:,:) ) 
    155             ENDIF 
     163 
     164            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  ) 
     165            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) ) 
    156166         ENDIF 
    157167         !                          ! ================== 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90

    r6140 r7403  
    2020   USE traldf_iso     ! lateral diffusion (Madec operator)         (tra_ldf_iso routine) 
    2121   USE diaptr         ! poleward transport diagnostics 
     22   USE diaar5         ! AR5 diagnostics 
    2223   USE zpshde         ! partial step: hor. derivative     (zps_hde routine) 
    2324   ! 
    … …  
    3536 
    3637   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels 
     38 
     39   LOGICAL  ::   l_ptr   ! flag to compute poleward transport 
     40   LOGICAL  ::   l_hst   ! flag to compute heat transport 
     41 
    3742 
    3843   !! * Substitutions 
    … …  
    8994      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 
    9095      REAL(wp) ::   zah, zah_slp, zaei_slp 
    91 #if defined key_diaar5 
    92       REAL(wp) ::   zztmp              ! local scalar 
    93 #endif 
    9496      REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d                                            ! 2D workspace 
    9597      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     - 
    … …  
    112114         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 
    113115      ENDIF 
    114       !                                               ! set time step size (Euler/Leapfrog) 
     116      !    
     117      l_hst = .FALSE. 
     118      l_ptr = .FALSE. 
     119      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE.  
     120      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     121         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
     122      ! 
     123      !                                                        ! set time step size (Euler/Leapfrog) 
    115124      IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    116125      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
    … …  
    416425            ! 
    417426            !                          ! "Poleward" diffusive heat or salt transports (T-S case only) 
    418             IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    419                IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names 
    420                IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) ) 
    421             ENDIF 
    422             ! 
    423             IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
    424               ! 
    425               IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    426                   z2d(:,:) = zftu(ji,jj,1)  
    427                   DO jk = 2, jpkm1 
    428                      DO jj = 2, jpjm1 
    429                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    430                            z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
    431                         END DO 
    432                      END DO 
    433                   END DO 
    434                   z2d(:,:) = rau0_rcp * z2d(:,:)  
    435                   CALL lbc_lnk( z2d, 'U', -1. ) 
    436                   CALL iom_put( "udiff_heattr", z2d )                  ! heat i-transport 
    437                   ! 
    438                   z2d(:,:) = zftv(ji,jj,1)  
    439                   DO jk = 2, jpkm1 
    440                      DO jj = 2, jpjm1 
    441                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    442                            z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
    443                         END DO 
    444                      END DO 
    445                   END DO 
    446                   z2d(:,:) = rau0_rcp * z2d(:,:)      
    447                   CALL lbc_lnk( z2d, 'V', -1. ) 
    448                   CALL iom_put( "vdiff_heattr", z2d )                  !  heat j-transport 
    449                ENDIF 
    450                ! 
    451             ENDIF 
     427            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  ) 
     428            !                          ! Diffusive heat transports 
     429            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) ) 
    452430            ! 
    453431         ENDIF                                                    !== end pass selection  ==! 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90

    r6140 r7403  
    3333# endif 
    3434   !                                                  !!!* Active tracers trends indexes 
    35    INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 14     !: Total trend nb: change it when adding/removing one indice below 
     35   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 20     !: Total trend nb: change it when adding/removing one indice below 
    3636   !                               ===============     !   
    3737   INTEGER, PUBLIC, PARAMETER ::   jptra_xad  =  1     !: x- horizontal advection 
    … …  
    3939   INTEGER, PUBLIC, PARAMETER ::   jptra_zad  =  3     !: z- vertical   advection 
    4040   INTEGER, PUBLIC, PARAMETER ::   jptra_sad  =  4     !: z- vertical   advection 
    41    INTEGER, PUBLIC, PARAMETER ::   jptra_ldf  =  5     !: lateral       diffusion 
    42    INTEGER, PUBLIC, PARAMETER ::   jptra_zdf  =  6     !: vertical      diffusion 
    43    INTEGER, PUBLIC, PARAMETER ::   jptra_zdfp =  7     !: "PURE" vert.  diffusion (ln_traldf_iso=T) 
    44    INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  =  8     !: Bottom Boundary Condition (geoth. heating)  
    45    INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  =  9     !: Bottom Boundary Layer (diffusive and/or advective) 
    46    INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 10     !: non-penetrative convection treatment 
    47    INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 11     !: internal restoring (damping) 
    48    INTEGER, PUBLIC, PARAMETER ::   jptra_qsr  = 12     !: penetrative solar radiation 
    49    INTEGER, PUBLIC, PARAMETER ::   jptra_nsr  = 13     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T) 
    50    INTEGER, PUBLIC, PARAMETER ::   jptra_atf  = 14     !: Asselin time filter 
     41   INTEGER, PUBLIC, PARAMETER ::   jptra_totad  =  5   !: total         advection 
     42   INTEGER, PUBLIC, PARAMETER ::   jptra_ldf  =  6     !: lateral       diffusion 
     43   INTEGER, PUBLIC, PARAMETER ::   jptra_zdf  =  7     !: vertical      diffusion 
     44   INTEGER, PUBLIC, PARAMETER ::   jptra_zdfp =  8     !: "PURE" vert.  diffusion (ln_traldf_iso=T) 
     45   INTEGER, PUBLIC, PARAMETER ::   jptra_evd  =  9     !: EVD term (convection) 
     46   INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  = 10     !: Bottom Boundary Condition (geoth. heating)  
     47   INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  = 11     !: Bottom Boundary Layer (diffusive and/or advective) 
     48   INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 12     !: non-penetrative convection treatment 
     49   INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 13     !: internal restoring (damping) 
     50   INTEGER, PUBLIC, PARAMETER ::   jptra_qsr  = 14     !: penetrative solar radiation 
     51   INTEGER, PUBLIC, PARAMETER ::   jptra_nsr  = 15     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T) 
     52   INTEGER, PUBLIC, PARAMETER ::   jptra_atf  = 16     !: Asselin time filter 
     53   INTEGER, PUBLIC, PARAMETER ::   jptra_tot  = 17     !: Model total trend 
    5154   ! 
    5255   !                                                  !!!* Passive tracers trends indices (use if "key_top" defined) 
    53    INTEGER, PUBLIC, PARAMETER ::   jptra_sms  = 15     !: sources m. sinks 
    54    INTEGER, PUBLIC, PARAMETER ::   jptra_radn = 16     !: corr. trn<0 in trcrad 
    55    INTEGER, PUBLIC, PARAMETER ::   jptra_radb = 17     !: corr. trb<0 in trcrad (like atf) 
     56   INTEGER, PUBLIC, PARAMETER ::   jptra_sms  = 18     !: sources m. sinks 
     57   INTEGER, PUBLIC, PARAMETER ::   jptra_radn = 19     !: corr. trn<0 in trcrad 
     58   INTEGER, PUBLIC, PARAMETER ::   jptra_radb = 20     !: corr. trb<0 in trcrad (like atf) 
    5659   ! 
    5760   !                                                  !!!* Momentum trends indices 
    58    INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 15     !: Total trend nb: change it when adding/removing one indice below 
     61   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 13     !: Total trend nb: change it when adding/removing one indice below 
    5962   !                               ===============     !   
    6063   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg  =  1     !: hydrostatic pressure gradient  

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90

    r6140 r7403  
    9090!!gm end 
    9191      ! 
    92       IF( .NOT.ln_linssh .AND. ( l_trdtra .OR. l_trddyn ) )  CALL ctl_stop( 'trend diagnostics with variable volume not validated' ) 
     92!      IF( .NOT.ln_linssh .AND. ( l_trdtra .OR. l_trddyn ) )  CALL ctl_stop( 'trend diagnostics with variable volume not validated' ) 
    9393       
    9494!!gm  : Potential BUG : 3D output only for vector invariant form!  add a ctl_stop or code the flux form case 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90

    r6140 r7403  
    2828   USE lib_mpp        ! MPP library 
    2929   USE wrk_nemo       ! Memory allocation 
     30   USE ldfslp         ! Isopycnal slopes 
    3031 
    3132   IMPLICIT NONE 
    … …  
    180181!                              CALL iom_put( "ketrd_bfri", zke2d ) 
    181182!         ENDIF 
    182          CASE( jpdyn_ken )   ;                                          ! kinetic energy 
    183                                  ! called in dynnxt.F90 before asselin time filter with putrd=ua and pvtrd=va 
    184                                  zke(:,:,:) = 0.5_wp * zke(:,:,:) 
    185                                  CALL iom_put( "KE", zke ) 
    186                                  ! 
    187                                  CALL ken_p2k( kt , zke ) 
    188                                  CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w 
     183        CASE( jpdyn_ken )   ;   ! kinetic energy 
     184                    ! called in dynnxt.F90 before asselin time filter 
     185                    ! with putrd=ua and pvtrd=va 
     186                    zke(:,:,:) = 0.5_wp * zke(:,:,:) 
     187                    CALL iom_put( "KE", zke ) 
     188                    ! 
     189                    CALL ken_p2k( kt , zke ) 
     190                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w 
    189191         ! 
    190192      END SELECT 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r6140 r7403  
    3939 
    4040   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt   ! use to store the temperature trends 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_evd  ! store avt_evd to calculate EVD trend 
    4142 
    4243   !! * Substitutions 
    … …  
    5455      !!                  ***  FUNCTION trd_tra_alloc  *** 
    5556      !!--------------------------------------------------------------------- 
    56       ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc ) 
     57      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc ) 
    5758      ! 
    5859      IF( lk_mpp             )   CALL mpp_sum ( trd_tra_alloc ) 
    … …  
    127128            zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp 
    128129            DO jk = 2, jpk 
    129                zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk) 
     130               zwt(:,:,jk) = avt_k(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk) 
    130131               zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk) 
    131132            END DO 
    … …  
    137138            END DO 
    138139            CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )   
     140            ! 
     141            !                         ! Also calculate EVD trend at this point.  
     142            zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes 
     143            DO jk = 2, jpk 
     144               zwt(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / e3w_n(:,:,jk) * tmask(:,:,jk) 
     145               zws(:,:,jk) = avt_evd(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / e3w_n(:,:,jk) * tmask(:,:,jk) 
     146            END DO 
     147            ! 
     148            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp 
     149            DO jk = 1, jpkm1 
     150               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t_n(:,:,jk) 
     151               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t_n(:,:,jk)  
     152            END DO 
     153            CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt )   
    139154            ! 
    140155            CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt ) 
    … …  
    311326                                  CALL wrk_dealloc( jpi, jpj, z2dx, z2dy ) 
    312327                               ENDIF 
     328      CASE( jptra_totad  ) ;   CALL iom_put( "ttrd_totad" , ptrdx )      ! total   advection 
     329                               CALL iom_put( "strd_totad" , ptrdy ) 
    313330      CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion 
    314331                               CALL iom_put( "strd_ldf" , ptrdy ) 
    … …  
    317334      CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution) 
    318335                               CALL iom_put( "strd_zdfp", ptrdy ) 
     336      CASE( jptra_evd )    ;   CALL iom_put( "ttrd_evd", ptrdx )         ! EVD trend (convection) 
     337                               CALL iom_put( "strd_evd", ptrdy ) 
    319338      CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping) 
    320339                               CALL iom_put( "strd_dmp" , ptrdy ) 
    … …  
    323342      CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing 
    324343                               CALL iom_put( "strd_npc" , ptrdy ) 
    325       CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx )        ! surface forcing + runoff (ln_rnf=T) 
    326                                CALL iom_put( "strd_cdt" , ptrdy ) 
     344      CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx(:,:,1) )        ! surface forcing + runoff (ln_rnf=T) 
     345                               CALL iom_put( "strd_cdt" , ptrdy(:,:,1) )        ! output as 2D surface fields 
    327346      CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature) 
    328347      CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature) 
    329348      CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter 
    330349                               CALL iom_put( "strd_atf" , ptrdy ) 
     350      CASE( jptra_tot  )   ;   CALL iom_put( "ttrd_tot" , ptrdx )        ! model total trend 
     351                               CALL iom_put( "strd_tot" , ptrdy ) 
    331352      END SELECT 
    332353      ! 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r5836 r7403  
    3535   INTEGER , PUBLIC ::   nn_npc      !: non penetrative convective scheme call  frequency 
    3636   INTEGER , PUBLIC ::   nn_npcp     !: non penetrative convective scheme print frequency 
     37   LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag 
    3738 
    3839 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r6140 r7403  
    1717   USE dom_oce         ! ocean space and time domain variables 
    1818   USE zdf_oce         ! ocean vertical physics variables 
     19   USE trd_oce         ! trends: ocean variables 
     20   USE trdtra          ! trends manager: tracers  
    1921   ! 
    2022   USE in_out_manager  ! I/O manager 
    … …  
    111113      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
    112114      CALL iom_put( "avt_evd", zavt_evd )              ! output this change 
     115      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 
    113116      ! 
    114117      CALL wrk_dealloc( jpi,jpj,jpk,   zavt_evd, zavm_evd )  

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r5836 r7403  
    5151      INTEGER ::   ioptio, ios       ! local integers 
    5252      !! 
    53       NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,   & 
    54          &              ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp 
     53      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  & 
     54         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp,       & 
     55         &        ln_zdfqiao 
    5556      !!---------------------------------------------------------------------- 
    5657 
    … …  
    8182         WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc 
    8283         WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp 
     84         WRITE(numout,*) '      Qiao formulation flag               ln_zdfqiao=', ln_zdfqiao 
    8385      ENDIF 
    8486 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r7048 r7403  
    205205         DO jj = 2, jpjm1 
    206206            DO ji = fs_2, fs_jpim1 
    207                IF( fsdept(ji,jj,jk) < ekm_dep(ji,jj) ) THEN 
     207               IF( gdepw_n(ji,jj,jk) < ekm_dep(ji,jj) ) THEN 
    208208                  avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix ) 
    209209                  avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix ) 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r6152 r7403  
    490490      IF( lk_floats     )   CALL     flo_init   ! drifting Floats 
    491491                            CALL dia_cfl_init   ! Initialise CFL diagnostics 
    492       IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag 
    493492                            CALL dia_ptr_init   ! Poleward TRansports initialization 
    494493      IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6464 r7403  
    2626   !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs 
    2727   !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state 
     28   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 
    2829   !!            3.7  !  2014-10  (G. Madec)  LDF simplication  
    2930   !!             -   !  2014-12  (G. Madec) remove KPP scheme 
    … …  
    7374      !!              -8- Outputs and diagnostics 
    7475      !!---------------------------------------------------------------------- 
    75       INTEGER ::   jk      ! dummy loop indice 
     76      INTEGER ::   ji,jj,jk ! dummy loop indice 
    7677      INTEGER ::   indic    ! error indicator if < 0 
    7778      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt) 
    … …  
    128129                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic) 
    129130      !                                               ! Vertical eddy viscosity and diffusivity coefficients 
    130       IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz 
    131       IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz 
    132       IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz 
    133       IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value) 
     131      IF( lk_zdfric  )   CALL zdf_ric ( kstp )             ! Richardson number dependent Kz 
     132      IF( lk_zdftke  )   CALL zdf_tke ( kstp )             ! TKE closure scheme for Kz 
     133      IF( lk_zdfgls  )   CALL zdf_gls ( kstp )             ! GLS closure scheme for Kz 
     134      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )             ! Qiao vertical mixing  
     135      ! 
     136      IF( lk_zdfcst  ) THEN                                ! Constant Kz (reset avt, avm[uv] to the background value) 
    134137         avt (:,:,:) = rn_avt0 * wmask (:,:,:) 
    135138         avmu(:,:,:) = rn_avm0 * wumask(:,:,:) 
    … …  
    207210                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
    208211                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
     212      IF( ln_wave .AND. ln_sdw .AND. ln_stcor)           & 
     213               &         CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    209214                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
    210215                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
    … …  
    234239      IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics 
    235240      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports 
    236       IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag 
     241                         CALL dia_ar5( kstp )         ! ar5 diag 
    237242      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    238243                         CALL dia_wri( kstp )         ! ocean model: outputs 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r6140 r7403  
    1919   USE sbcapr           ! surface boundary condition: atmospheric pressure 
    2020   USE sbctide          ! Tide initialisation 
     21   USE sbcwave          ! Wave intialisation 
    2122 
    2223   USE traqsr           ! solar radiation penetration      (tra_qsr routine) 
    … …  
    4142   USE dynzdf           ! vertical diffusion               (dyn_zdf routine) 
    4243   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
     44   USE dynstcor         ! simp. form of Stokes-Coriolis 
    4345 
    4446   USE dynnxt           ! time-stepping                    (dyn_nxt routine) 
    … …  
    7173   USE zdfric           ! Richardson vertical mixing       (zdf_ric routine) 
    7274   USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine) 
     75   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine) 
    7376 
    7477   USE step_diu        ! Time stepping for diurnal sst 

  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90

    r6140 r7403  
    2424   PUBLIC   trc_oce_alloc      ! function called by nemogcm.F90 
    2525 
     26   LOGICAL , PUBLIC                                      ::   l_co2cpl = .false.   !: atmospheric pco2 recieved from oasis 
    2627   INTEGER , PUBLIC                                      ::   nn_dttrc      !: frequency of step on passive tracers 
    2728   REAL(wp), PUBLIC                                      ::   r_si2         !: largest depth of extinction (blue & 0.01 mg.m-3)  (RGB) 
    2829   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   etot3         !: light absortion coefficient 
    29    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   facvol        !: volume for degraded regions 
     30   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   oce_co2   !: ocean carbon flux 
    3031 
    3132#if defined key_top  
    … …  
    7576      !!                  ***  trc_oce_alloc  *** 
    7677      !!---------------------------------------------------------------------- 
    77       INTEGER ::   ierr(2)        ! Local variables 
    78       !!---------------------------------------------------------------------- 
    79       ierr(:) = 0 
    80                      ALLOCATE( etot3 (jpi,jpj,jpk), STAT=ierr(1) ) 
    81       IF( lk_degrad) ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr(2) ) 
    82       trc_oce_alloc  = MAXVAL( ierr ) 
    83       ! 
     78      ALLOCATE( etot3(jpi,jpj,jpk), oce_co2(jpi,jpj), STAT=trc_oce_alloc ) 
    8479      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3 array') 
     80      ! 
    8581   END FUNCTION trc_oce_alloc 
    8682 
Note: See TracChangeset for help on using the changeset viewer.