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 1975 for branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2010-06-28T19:22:14+02:00 (14 years ago)
Author:
mlelod
Message:

ticket: #663 MLF: first part

Location:
branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90

    r1870 r1975  
    2525   USE oce             ! ocean dynamics and tracers variables 
    2626   USE dom_oce         ! ocean space and time domain variables  
     27   USE sbc_oce         ! surface boundary condition: ocean 
    2728   USE zdf_oce         ! ??? 
    2829   USE domvvl          ! variable volume 
     
    3738   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3839   USE prtctl          ! Print control 
     40   USE traqsr          ! penetrative solar radiation (needed for nksr) 
    3941   USE agrif_opa_update 
    4042   USE agrif_opa_interp 
     
    4547   PUBLIC   tra_nxt    ! routine called by step.F90 
    4648 
    47    REAL(wp)                 ::   rbcp, r1_2bcp   ! Brown & Campana parameters for semi-implicit hpg 
     49   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg 
    4850   REAL(wp), DIMENSION(jpk) ::   r2dt_t          ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
    4951 
     
    98100         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    99101         ! 
    100          rbcp    = 0.25 * (1 + atfp) * (1 + atfp * atfp)       ! Brown & Campana parameter for semi-implicit hpg 
    101          r1_2bcp = 1.e0 - 2.e0 * rbcp 
     102         rbcp    = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)       ! Brown & Campana parameter for semi-implicit hpg 
    102103      ENDIF 
    103104      !                                      ! set time step size (Euler/Leapfrog) 
     
    164165      !!              - swap tracer fields to prepare the next time_step. 
    165166      !!                This can be summurized for temperature as: 
    166       !!             ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)       ln_dynhpg_imp = T 
     167      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T 
    167168      !!             ztm = 0                                   otherwise 
    168       !!                   with rbcp=(1+atfp)*(1+atfp*atfp)/4 
     169      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp) 
    169170      !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
    170171      !!             tn  = ta  
     
    177178      !! 
    178179      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
    179       REAL(wp) ::   ztm, ztf, zac   ! temporary scalars 
    180       REAL(wp) ::   zsm, zsf        !    -         - 
     180      REAL(wp) ::   zt_m, zs_m      ! temporary scalars 
     181      REAL(wp) ::   ztn, zsn        !    -         - 
    181182      !!---------------------------------------------------------------------- 
    182183 
     
    187188      ENDIF 
    188189      ! 
    189       !                                              ! ----------------------- ! 
    190       IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case ! 
    191          !                                           ! ----------------------- ! 
    192          ! 
    193          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    194             DO jk = 1, jpkm1                              ! (only swap) 
    195                tn(:,:,jk) = ta(:,:,jk)                                                          ! tb <-- tn 
    196                sn(:,:,jk) = sa(:,:,jk) 
    197             END DO 
    198          ELSE                                             ! general case (Leapfrog + Asselin filter) 
    199             DO jk = 1, jpkm1 
    200                DO jj = 1, jpj 
    201                   DO ji = 1, jpi 
    202                      ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk)      ! mean t 
    203                      zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk)  
    204                      ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk)   - 2.      * tn(ji,jj,jk) )    ! Asselin filter on t  
    205                      zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk)   - 2.      * sn(ji,jj,jk) ) 
    206                      tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                          ! tb <-- filtered tn  
    207                      sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
    208                      tn(ji,jj,jk) = ta(ji,jj,jk)                                                ! tn <-- ta 
    209                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    210                      ta(ji,jj,jk) = ztm                                                         ! ta <-- mean t 
    211                      sa(ji,jj,jk) = zsm 
    212                   END DO 
     190      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     191         DO jk = 1, jpkm1                                                 ! (only swap) 
     192            tn(:,:,jk) = ta(:,:,jk)                                       ! tn <-- ta 
     193            sn(:,:,jk) = sa(:,:,jk)                                       ! sn <-- sa 
     194         END DO 
     195      ELSE                                             ! General case (Leapfrog + Asselin filter) 
     196         DO jk = 1, jpkm1 
     197            DO jj = 1, jpj 
     198               DO ji = 1, jpi 
     199                  IF( ln_dynhpg_imp ) THEN                  ! implicit hpg: keep tn, sn in memory 
     200                     ztn = tn(ji,jj,jk) 
     201                     zsn = sn(ji,jj,jk) 
     202                  ENDIF 
     203                  !                                         ! time laplacian on tracers 
     204                  !                                         ! used for both Asselin and Brown & Campana filters 
     205                  zt_m = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk) 
     206                  zs_m = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk) 
     207                  ! 
     208                  !                                         ! swap of arrays 
     209                  tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_m               ! tb <-- tn filtered 
     210                  sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_m               ! sb <-- sn filtered 
     211                  tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta 
     212                  sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa 
     213                  !                                         ! semi imlicit hpg computation (Brown & Campana) 
     214                  IF( ln_dynhpg_imp ) THEN 
     215                     ta(ji,jj,jk) = ztn + rbcp * zt_m                     ! ta <-- Brown & Campana average 
     216                     sa(ji,jj,jk) = zsn + rbcp * zs_m                     ! sa <-- Brown & Campana average 
     217                  ENDIF 
    213218               END DO 
    214219            END DO 
    215          ENDIF 
    216          !                                           ! ----------------------- ! 
    217       ELSE                                           !    explicit hpg case    ! 
    218          !                                           ! ----------------------- ! 
    219          ! 
    220          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    221             DO jk = 1, jpkm1 
    222                tn(:,:,jk) = ta(:,:,jk)                                                    ! tn <-- ta 
    223                sn(:,:,jk) = sa(:,:,jk) 
    224             END DO 
    225          ELSE                                             ! general case (Leapfrog + Asselin filter) 
    226             DO jk = 1, jpkm1 
    227                DO jj = 1, jpj 
    228                   DO ji = 1, jpi 
    229                      ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )      ! Asselin filter on t  
    230                      zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 
    231                      tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                    ! tb <-- filtered tn  
    232                      sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 
    233                      tn(ji,jj,jk) = ta(ji,jj,jk)                                          ! tn <-- ta 
    234                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    235                   END DO 
    236                END DO 
    237             END DO 
    238          ENDIF 
    239       ENDIF 
    240       ! 
    241 !!gm from Matthieu : unchecked 
    242       IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter 
    243          zac = atfp * rdttra(1) 
    244          tb(:,:,1) = tb(:,:,1) - zac * ( qns (:,:) - qns_b (:,:) )         ! non solar surface flux 
    245          sb(:,:,1) = sn(:,:,1) - zac * ( emps(:,:) - emps_b(:,:) )         ! surface salt flux 
    246          ! 
    247          IF( ln_traqsr )                                                   ! solar penetrating flux 
    248             DO jk = 1, nksr 
    249                   DO jj = 1, jpj 
    250                      DO ji = 1, jpi 
    251                         IF( ln_sco ) THEN 
    252                            z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 
    253                         ELSEIF( ln_zps .OR. ln_zco ) THEN 
    254                            z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   & 
    255                                 &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 
    256                         ENDIF 
    257                         zt = zt - zfact1 * z_cor_qsr 
    258                      END DO 
    259                   END DO 
     220         END DO 
    260221      ENDIF 
    261222      ! 
     
    276237      !!              - swap tracer fields to prepare the next time_step. 
    277238      !!                This can be summurized for temperature as: 
    278       !!             ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb )   ln_dynhpg_imp = T 
    279       !!                  /( rbcp*e3t_a    + (2-rbcp)*e3t_n    + rbcp*e3t_b    )    
     239      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
     240      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )    
    280241      !!             ztm = 0                                                       otherwise 
    281242      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
     
    287248      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    288249      !!---------------------------------------------------------------------- 
    289       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     250      INTEGER, INTENT(in) ::   kt                  ! ocean time-step index 
    290251      !!      
    291       INTEGER  ::   ji, jj, jk             ! dummy loop indices 
    292       REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar 
    293       REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         - 
    294       REAL(wp) ::   ze3mr, ze3fr                           !    -         - 
    295       REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         - 
    296       !!---------------------------------------------------------------------- 
    297  
     252      INTEGER  ::   ji, jj, jk                     ! dummy loop indices 
     253      REAL     ::   ze3t_a, ze3t_n, ze3t_b         ! temporary scalars 
     254      REAL     ::   ztc_a, ztc_n, ztc_b            !    -         - 
     255      REAL     ::   zsc_a, zsc_n, zsc_b            !    -         - 
     256      REAL     ::   ztc_f, zsc_f, ztc_m, zsc_m     !    -         - 
     257      REAL     ::   ze3t_f, ze3t_m                 !    -         - 
     258      REAL     ::   zfact1, zfact2                 !    -         - 
     259      !!---------------------------------------------------------------------- 
     260       
    298261      IF( kt == nit000 ) THEN 
    299262         IF(lwp) WRITE(numout,*) 
     
    302265      ENDIF 
    303266 
    304       !                                              ! ----------------------- ! 
    305       IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !  (mean tracer computation) 
    306          !                                           ! ----------------------- ! 
    307          ! 
    308          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    309             DO jk = 1, jpkm1                                   ! (only swap) 
    310                tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta 
    311                sn(:,:,jk) = sa(:,:,jk) 
    312             END DO 
    313             !                                             ! general case (Leapfrog + Asselin filter) 
    314          ELSE                                             ! apply filter on thickness weighted tracer and swap 
    315             DO jk = 1, jpkm1 
    316                DO jj = 1, jpj 
    317                   DO ji = 1, jpi 
    318                      ze3t_b = fse3t_b(ji,jj,jk) 
    319                      ze3t_n = fse3t_n(ji,jj,jk) 
    320                      ze3t_a = fse3t_a(ji,jj,jk) 
    321                      !                                         ! tracer content at Before, now and after 
    322                      ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b 
    323                      ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n 
    324                      ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a 
    325                      ! 
    326                      !                                         ! Asselin filter on thickness and tracer content 
    327                      ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n ) 
    328                      ztc_f  = atfp * ( ztca   + ztcb   - 2.* ztcn   )  
    329                      zsc_f  = atfp * ( zsca   + zscb   - 2.* zscn   )  
    330                      ! 
    331                      !                                         ! filtered tracer including the correction  
    332                      ze3fr = 1.e0  / ( ze3t_n + ze3t_f ) 
    333                      ztf   = ze3fr * ( ztcn   + ztc_f  ) 
    334                      zsf   = ze3fr * ( zscn   + zsc_f  ) 
    335                      !                                         ! mean thickness and tracer 
    336                      ze3mr = 1.e0  / (  rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n ) 
    337                      ztm   = ze3mr * (  rbcp * ( ztca   + ztcb   ) + r1_2bcp * ztcn   ) 
    338                      zsm   = ze3mr * (  rbcp * ( zsca   + zscb   ) + r1_2bcp * zscn   ) 
    339 !!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !! 
    340 !!gm                 e3t_m(ji,jj,jk) = 1.e0 / ze3mr 
    341                      !                                         ! swap of arrays 
    342                      tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter 
    343                      sb(ji,jj,jk) = zsf 
    344                      tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta 
    345                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    346                      ta(ji,jj,jk) = ztm                             ! ta <-- mean t 
    347                      sa(ji,jj,jk) = zsm 
    348                   END DO 
     267      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     268         !                                                  ! (only swap) 
     269         DO jk = 1, jpkm1                                   ! tn <-- ta 
     270            tn(:,:,jk) = ta(:,:,jk)                         ! sn <-- sa 
     271            sn(:,:,jk) = sa(:,:,jk) 
     272         END DO 
     273         !                                             ! General case (Leapfrog + Asselin filter) 
     274      ELSE                                             ! apply filter on thickness weighted tracer and swap 
     275         DO jk = 1, jpkm1 
     276            zfact1 = atfp * r2dt_t(jk) 
     277            zfact2 = zfact1 / rau0 
     278            DO jj = 1, jpj 
     279               DO ji = 1, jpi 
     280                  !                                         ! scale factors at Before, now and after 
     281                  ze3t_b = fse3t_b(ji,jj,jk) 
     282                  ze3t_n = fse3t_n(ji,jj,jk) 
     283                  ze3t_a = fse3t_a(ji,jj,jk) 
     284                  !                                         ! tracer content at Before, now and after 
     285                  ztc_b  = tb(ji,jj,jk) * ze3t_b   ;   zsc_b = sb(ji,jj,jk) * ze3t_b 
     286                  ztc_n  = tn(ji,jj,jk) * ze3t_n   ;   zsc_n = sn(ji,jj,jk) * ze3t_n 
     287                  ztc_a  = ta(ji,jj,jk) * ze3t_a   ;   zsc_a = sa(ji,jj,jk) * ze3t_a 
     288                  ! 
     289                  !                                         ! Time laplacian on thickness and tracer content 
     290                  !                                         ! used for both Asselin and Brown & Campana filters 
     291                  ze3t_m = ze3t_a - 2. * ze3t_n + ze3t_b 
     292                  ztc_m  = ztc_a  - 2. * ztc_n  + ztc_b 
     293                  zsc_m  = zsc_a  - 2. * zsc_n  + zsc_b 
     294                  !                                         ! Asselin Filter + correction 
     295                  ze3t_f = ze3t_n + atfp * ze3t_m 
     296                  ztc_f  = ztc_n  + atfp * ztc_m 
     297                  zsc_f  = zsc_n  + atfp * zsc_m 
     298 
     299                  IF( jk == 1 ) THEN 
     300                     ze3t_f = ze3t_f - zfact2 * ( emp_b       (ji,jj)    - emp         (ji,jj)    ) 
     301                     ztc_f  = ztc_f  - zfact1 * ( sbc_trd_hc_n(ji,jj)    - sbc_trd_hc_b(ji,jj)    ) 
     302                  ENDIF 
     303                  IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN 
     304                     ztc_f  = ztc_f  - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) ) 
     305                  ENDIF 
     306                  !                                         ! swap of arrays 
     307                  ze3t_f = 1.e0 / ze3t_f 
     308                  tb(ji,jj,jk) = ztc_f * ze3t_f                           ! tb <-- tn filtered 
     309                  sb(ji,jj,jk) = zsc_f * ze3t_f                           ! sb <-- sn filtered 
     310                  tn(ji,jj,jk) = ta(ji,jj,jk)                             ! tn <-- ta 
     311                  sn(ji,jj,jk) = sa(ji,jj,jk)                             ! sn <-- sa 
     312                  !                                         ! semi imlicit hpg computation (Brown & Campana) 
     313                  IF( ln_dynhpg_imp ) THEN 
     314                     ze3t_m = 1.e0 / ( ze3t_n + rbcp * ze3t_m ) 
     315                     ta(ji,jj,jk) =  ( ztc_n  + rbcp * ztc_m  ) * ze3t_m  ! ta <-- Brown & Campana average 
     316                     sa(ji,jj,jk) =  ( zsc_n  + rbcp * zsc_m  ) * ze3t_m  ! sa <-- Brown & Campana average 
     317                  ENDIF 
    349318               END DO 
    350319            END DO 
    351          ENDIF 
    352          !                                           ! ----------------------- ! 
    353       ELSE                                           !    explicit hpg case    !  (NO mean tracer computation) 
    354          !                                           ! ----------------------- ! 
    355          ! 
    356          IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step 
    357             DO jk = 1, jpkm1                                   ! (only swap) 
    358                tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta 
    359                sn(:,:,jk) = sa(:,:,jk) 
    360             END DO 
    361             !                                             ! general case (Leapfrog + Asselin filter) 
    362          ELSE                                             ! apply filter on thickness weighted tracer and swap 
    363             DO jk = 1, jpkm1 
    364                DO jj = 1, jpj 
    365                   DO ji = 1, jpi 
    366                      ze3t_b = fse3t_b(ji,jj,jk) 
    367                      ze3t_n = fse3t_n(ji,jj,jk) 
    368                      ze3t_a = fse3t_a(ji,jj,jk) 
    369                      !                                         ! tracer content at Before, now and after 
    370                      ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b 
    371                      ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n 
    372                      ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a 
    373                      ! 
    374                      !                                         ! Asselin filter on thickness and tracer content 
    375                      ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 
    376                      ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   )  
    377                      zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   )  
    378                      ! 
    379                      !                                         ! filtered tracer including the correction  
    380                      ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 
    381                      ztf   =  ( ztcn  + ztc_f ) * ze3fr 
    382                      zsf   =  ( zscn  + zsc_f ) * ze3fr 
    383                      !                                         ! swap of arrays 
    384                      tb(ji,jj,jk) = ztf                             ! tb <-- tn filtered 
    385                      sb(ji,jj,jk) = zsf 
    386                      tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta 
    387                      sn(ji,jj,jk) = sa(ji,jj,jk) 
    388                   END DO 
    389                END DO 
    390             END DO 
    391          ENDIF 
    392       ENDIF 
    393 !!gm from Matthieu : unchecked 
    394       ! 
    395       IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter 
    396                IF( lk_vvl ) THEN                                        ! Varying levels 
    397                   DO jj = 1, jpj 
    398                      DO ji = 1, jpi 
    399                         ! - ML - modification for varaiance diagnostics 
    400                         zssh1 = tmask(ji,jj,jk) / ( fse3t(ji,jj,jk) + atfp * ( zfse3ta(ji,jj,jk) - 2*fse3t(ji,jj,jk) & 
    401                            &                                                 + fse3tb(ji,jj,jk) ) ) 
    402                         zt = zssh1 * ( tn(ji,jj,jk) +  atfp * ( ta(ji,jj,jk) - 2 * tn(ji,jj,jk) + tb(ji,jj,jk) ) ) 
    403                         zs = zssh1 * ( sn(ji,jj,jk) +  atfp * ( sa(ji,jj,jk) - 2 * sn(ji,jj,jk) + sb(ji,jj,jk) ) ) 
    404                         ! filter correction for global conservation 
    405                         !------------------------------------------ 
    406                         zfact1 = atfp * rdttra(1) * zssh1 
    407                         IF (jk == 1) THEN                        ! remove sbc trend from time filter 
    408                            zt = zt - zfact1 * ( qn(ji,jj) - qb(ji,jj) ) 
    409 !!???                           zs = zs - zfact1 * ( - emp( ji,jj) + empb(ji,jj) ) * zsrau * 35.5e0 
    410                         ENDIF 
    411                         ! remove qsr trend from time filter 
    412                         IF (jk <= nksr) THEN 
    413                            IF( ln_sco ) THEN 
    414                               z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 
    415                            ELSEIF( ln_zps .OR. ln_zco ) THEN 
    416                               z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   & 
    417                                    &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 
    418                            ENDIF 
    419                            zt = zt - zfact1 * z_cor_qsr 
    420                            IF (jk == nksr) qsrb(ji,jj) = qsr(ji,jj) 
    421                         ENDIF 
    422                         ! - ML - end of modification 
    423                         zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 
    424                         tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 
    425                         sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 
    426                      END DO 
    427                   END DO 
    428       ENDIF 
    429 !!gm end 
     320         END DO 
     321      ENDIF 
    430322      ! 
    431323   END SUBROUTINE tra_nxt_vvl 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90

    r1870 r1975  
    4747   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    4848    
    49    INTEGER                   ::   nksr    ! levels below which the light cannot penetrate (depth larger than 391 m) 
     49   INTEGER , PUBLIC          ::   nksr    ! levels below which the light cannot penetrate (depth larger than 391 m) 
    5050   REAL(wp), DIMENSION(3,61) ::   rkrgb   ! tabulated attenuation coefficients for RGB absorption 
    5151 
     
    9797      REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars 
    9898      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         - 
    99       REAL(wp) ::   zqsr                 !    -         - 
    10099      REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace 
    101100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace 
     
    123122               DO ji = fs_2, fs_jpim1   ! vector opt. 
    124123!!gm  how to stecify the mean of time step here : TOP versus OPA time stepping strategy not obvious 
    125                   ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
     124                  qsr_trd_hc_n(ji,jj,jk) = ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 
    126125               END DO 
    127126            END DO 
     
    155154               zsi0r = 1.e0 / rn_si0 
    156155               zcoef = ( 1. - rn_abs ) / 3.e0                           ! equi-partition in R-G-B 
    157 !!gm bug !!! zqsr only a constant not an array 
    158                zqsr  = 0.5 * ( qsr_b(:,:) + qsr(:,:) )                  ! mean over 2 time steps 
    159                ze0(:,:,1) = rn_abs  * zqsr 
    160                ze1(:,:,1) = zcoef   * zqsr 
    161                ze2(:,:,1) = zcoef   * zqsr 
    162                ze3(:,:,1) = zcoef   * zqsr 
    163                zea(:,:,1) =           zqsr 
     156 
     157               ze0(:,:,1) = rn_abs  * qsr(:,:) 
     158               ze1(:,:,1) = zcoef   * qsr(:,:) 
     159               ze2(:,:,1) = zcoef   * qsr(:,:) 
     160               ze3(:,:,1) = zcoef   * qsr(:,:) 
     161               zea(:,:,1) =           qsr(:,:) 
    164162               ! 
    165163               DO jk = 2, nksr+1                                     ! deeper values 
     
    183181               ! 
    184182               DO jk = 1, nksr                                       ! compute and add qsr trend to ta 
    185                   ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
     183                  qsr_trd_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
    186184               END DO 
    187185               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero 
     
    190188            ELSE                                                 !*  Constant Chlorophyll 
    191189               DO jk = 1, nksr 
    192                   ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) ) 
     190                  qsr_trd_hc_n(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 
    193191               END DO 
    194192            ENDIF 
     
    201199               DO jj = 2, jpjm1 
    202200                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    203                      ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) ) 
     201                     qsr_trd_hc_n(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) 
    204202                  END DO 
    205203               END DO 
     
    208206         ENDIF 
    209207         ! 
     208      ENDIF 
     209 
     210      ! Add qsr trend to ta in all cases 
     211      IF( neuler == 0 .AND. kt == nit000 ) THEN 
     212         DO jk = 1, nksr 
     213            DO jj = 2, jpjm1 
     214               DO ji = fs_2, fs_jpim1   ! vector opt. 
     215                  ta(ji,jj,jk) = ta(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) / fse3t(ji,jj,jk) 
     216               END DO 
     217            END DO 
     218         END DO 
     219      ELSE 
     220         DO jk = 1, nksr 
     221            DO jj = 2, jpjm1 
     222               DO ji = fs_2, fs_jpim1   ! vector opt. 
     223                  ta(ji,jj,jk) = ta(ji,jj,jk) + 0.5 * ( qsr_trd_hc_b(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) ) / fse3t(ji,jj,jk) 
     224               END DO 
     225            END DO 
     226         END DO 
    210227      ENDIF 
    211228 
     
    382399               ! 
    383400               DO jk = 1, nksr 
    384                   etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
     401                  etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
    385402               END DO 
    386403               etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero 
     
    404421                     zc0 = rn_abs * EXP( -fsdepw(ji,jj,jk  )*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk  )*zsi1r ) 
    405422                     zc1 = rn_abs * EXP( -fsdepw(ji,jj,jk+1)*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk+1)*zsi1r ) 
    406                      etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) / fse3t(ji,jj,jk) 
     423                     etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) 
    407424                  END DO 
    408425               END DO 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90

    r1870 r1975  
    104104      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
    105105      !! 
    106       INTEGER  ::   ji, jj                   ! dummy loop indices 
    107       REAL(wp) ::   zta, zsa, zsrau, zse3t   ! temporary scalars 
     106      INTEGER  ::   ji, jj            ! dummy loop indices 
     107      REAL(wp) ::   zsrau, zse3t      ! temporary scalars 
    108108      !!---------------------------------------------------------------------- 
    109109 
     
    129129         qsr(:,:) = 0.e0                     ! qsr set to zero 
    130130         IF( kt == nit000 ) THEN             ! idem on before field at nit000 
    131             qns_b(:,:) = qns_b(:,:) + qsr_b(:,:)   
    132131            qsr_b(:,:) = 0.e0                      
    133132         ENDIF 
    134133      ENDIF 
    135  
    136       ! Concentration dillution effect on (t,s) 
    137 !!      DO jj = 2, jpj 
    138 !!         DO ji = fs_2, fs_jpim1   ! vector opt. 
    139 !!#if ! defined key_zco 
    140 !!            zse3t = 1. / fse3t(ji,jj,1) 
    141 !!#endif 
    142 !!            IF( lk_vvl) THEN 
    143 !!               zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux 
    144 !!                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux 
    145 !!               zsa = 0.e0                                            ! No salinity concent./dilut. effect 
    146 !!            ELSE 
    147 !!               zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux 
    148 !!               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect 
    149 !!            ENDIF 
    150 !!            ta(ji,jj,1) = ta(ji,jj,1) + zta                          ! add the trend to the general tracer trend 
    151 !!            sa(ji,jj,1) = sa(ji,jj,1) + zsa 
    152 !!         END DO 
    153 !!      END DO 
    154134 
    155135      !                             ! ---------------------- ! 
    156136      IF( lk_vvl ) THEN             !  Variable Volume case  ! 
    157137         !                          ! ---------------------- ! 
    158          DO jj = 2, jpj 
    159             DO ji = fs_2, fs_jpim1   ! vector opt. 
    160                zse3t = 0.5 / fse3t(ji,jj,1) 
    161                !                           ! temperature: heat flux + cooling/heating effet of EMP flux 
    162                ta(ji,jj,1) = ta(ji,jj,1) + (  ro0cpr * ( qns(ji,jj) + qns_b(ji,jj) )                &  
    163                   &                         - zsrau  * ( emp(ji,jj) + emp_b(ji,jj) ) * tn(ji,jj,1)  ) * zse3t 
    164                !                           ! salinity: salt flux 
    165                sa(ji,jj,1) = sa(ji,jj,1) + ( emps(ji,jj) + emps_b(ji,jj) ) * zse3t 
    166  
    167 !!gm BUG : in key_vvl emps must be modified to only include the salt flux due to with sea-ice freezing/melting 
    168 !!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface. 
    169  
    170             END DO 
    171          END DO 
     138!!gm BUG : in key_vvl emps must be modified to only include the salt flux due to sea-ice freezing/melting 
     139!!gm       otherwise this flux will be missing  ==> modification required in limsbc,  limsbc_2 and CICE interface.s 
     140         IF ( neuler == 0 .AND. kt == nit000 ) THEN 
     141            DO jj = 2, jpj 
     142               DO ji = fs_2, fs_jpim1   ! vector opt. 
     143                  ! temperature : heat flux + cooling/heating effet of EMP flux 
     144                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1) 
     145#if ! defined key_zco 
     146                  zse3t = 1. / fse3t(ji,jj,1) 
     147#endif 
     148                  ta(ji,jj,1) = ta(ji,jj,1) + zse3t * sbc_trd_hc_n(ji,jj) 
     149                END DO 
     150            END DO 
     151         ELSE 
     152            DO jj = 2, jpj 
     153               DO ji = fs_2, fs_jpim1   ! vector opt. 
     154                  ! temperature : heat flux + cooling/heating effet of EMP flux 
     155                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1) 
     156#if ! defined key_zco 
     157                  zse3t = 1. / fse3t(ji,jj,1) 
     158#endif 
     159                  ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t 
     160               END DO 
     161            END DO 
     162         ENDIF 
    172163         !                          ! ---------------------- ! 
    173164      ELSE                          !  Constant Volume case  ! 
    174165         !                          ! ---------------------- ! 
    175          DO jj = 2, jpj 
    176             DO ji = fs_2, fs_jpim1   ! vector opt. 
    177                zse3t = 0.5 / fse3t(ji,jj,1) 
    178                !                           ! temperature: heat flux  
    179                ta(ji,jj,1) = ta(ji,jj,1) + ro0cpr * ( qns (ji,jj) + qns_b (ji,jj) )               * zse3t  
    180                !                           ! salinity: salt flux + concent./dilut. effect (both in emps) 
    181                sa(ji,jj,1) = sa(ji,jj,1) + zsrau  * ( emps(ji,jj) + emps_b(ji,jj) ) * sn(ji,jj,1) * zse3t 
    182             END DO 
    183          END DO 
     166         IF ( neuler == 0 .AND. kt == nit000 ) THEN 
     167            DO jj = 2, jpj 
     168               DO ji = fs_2, fs_jpim1   ! vector opt. 
     169                  ! temperature : heat flux 
     170                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 
     171                  ! salinity    : salt flux + concent./dilut. effect (both in emps) 
     172                  sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1) 
     173#if ! defined key_zco 
     174                  zse3t = 1. / fse3t(ji,jj,1) 
     175#endif 
     176                  ta(ji,jj,1) = ta(ji,jj,1) + sbc_trd_hc_n(ji,jj) * zse3t 
     177                  sa(ji,jj,1) = sa(ji,jj,1) + sbc_trd_sc_n(ji,jj) * zse3t 
     178               END DO 
     179            END DO 
     180         ELSE 
     181            DO jj = 2, jpj 
     182               DO ji = fs_2, fs_jpim1   ! vector opt. 
     183                  ! temperature : heat flux 
     184                  sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 
     185                  ! salinity    : salt flux + concent./dilut. effect (both in emps) 
     186                  sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1) 
     187#if ! defined key_zco 
     188                  zse3t = 1. / fse3t(ji,jj,1) 
     189#endif 
     190                  ! temperature : heat flux 
     191                  ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t  
     192                  sa(ji,jj,1) = sa(ji,jj,1) + 0.5 * ( sbc_trd_sc_b(ji,jj) + sbc_trd_sc_n(ji,jj) ) * zse3t 
     193               END DO 
     194            END DO 
     195         ENDIF 
    184196         ! 
    185197      ENDIF 
    186198 
    187       IF( l_trdtra ) THEN        ! save the sbc trends for diagnostic 
     199      IF( l_trdtra ) THEN           ! save the sbc trends for diagnostic 
    188200         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 
    189201         ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 
    190202         CALL trd_mod(ztrdt, ztrds, jptra_trd_nsr, 'TRA', kt) 
    191203      ENDIF 
    192       !                          ! control print 
     204      !                             ! control print 
    193205      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask,   & 
    194206         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
Note: See TracChangeset for help on using the changeset viewer.