Changeset 1975


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

ticket: #663 MLF: first part

Location:
branches/DEV_r1837_MLF/NEMO/OPA_SRC
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90

    r1740 r1975  
    244244                     zve3b = vb(ji,jj,jk) * ze3v_b 
    245245                     ! 
    246                      zuf  = ( atfp  * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   & 
    247                         & / ( atfp  * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 
    248                      zvf  = ( atfp  * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   & 
    249                         & / ( atfp  * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 
     246                     zuf  = ( atfp * ( zue3b  + zue3a  ) + atfp1 * zue3n  )   & 
     247                        & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 
     248                     zvf  = ( atfp * ( zve3b  + zve3a  ) + atfp1 * zve3n  )   & 
     249                        & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 
    250250                     ! 
    251251                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90

    r1870 r1975  
    7171      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
    7272      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars 
    73       REAL(wp) ::   z2dt, zraur     ! temporary scalars 
     73      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars 
    7474      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace 
    7575      REAL(wp), DIMENSION(jpi,jpj) ::   z2d         ! 2D workspace 
     
    9494                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
    9595                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
    96                   sshf_b(ji,jj) = zcoeff * ( sshb(ji  ,jj) + sshb(ji  ,jj+1)                 & 
    97                      &                     + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 
    9896                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
    9997                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
    10098                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
    10199                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
    102                   sshf_n(ji,jj) = zcoeff * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)                 & 
    103                      &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 
    104100               END DO 
    105101            END DO 
    106102            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
    107103            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     104            DO jj = 1, jpjm1 
     105               DO ji = 1, jpim1      ! NO Vector Opt. 
     106                  sshf_b(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
     107                     &                 / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     108                     &                 * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_b(ji,jj  )     & 
     109                     &                   + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) ) 
     110               END DO 
     111            END DO 
     112            DO jj = 1, jpjm1 
     113               DO ji = 1, jpim1      ! NO Vector Opt. 
     114                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
     115                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     116                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
     117                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
     118               END DO 
     119            END DO 
    108120            CALL lbc_lnk( sshf_b, 'F', 1. )   ;   CALL lbc_lnk( sshf_n, 'F', 1. ) 
    109121         ENDIF 
     
    145157      !                                           !   After Sea Surface Height   ! 
    146158      !                                           !------------------------------! 
    147       zraur = 0.5 / rau0 
     159      z1_rau0 = 0.5 / rau0 
    148160      ! 
    149161      zhdiv(:,:) = 0.e0 
     
    152164      END DO 
    153165      !                                                ! Sea surface elevation time stepping 
    154       ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     166      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 
     167      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 
     168      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    155169 
    156170#if defined key_obc 
     
    170184                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     & 
    171185                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 
    172                sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1)                                 &  
    173                   &                                  * ( ssha(ji  ,jj) + ssha(ji  ,jj+1)                 & 
    174                   &                                    + ssha(ji+1,jj) + ssha(ji+1,jj+1) ) 
    175186            END DO 
    176187         END DO 
    177          CALL lbc_lnk( sshu_a, 'U', 1. )               ! Boundaries conditions 
     188         ! Boundaries conditions 
     189         CALL lbc_lnk( sshu_a, 'U', 1. ) 
    178190         CALL lbc_lnk( sshv_a, 'V', 1. ) 
     191         DO jj = 1, jpjm1 
     192            DO ji = 1, jpim1      ! NO Vector Opt. 
     193               sshf_a(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
     194                    &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     195                    &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_a(ji,jj  )     & 
     196                    &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_a(ji,jj+1) ) 
     197            END DO 
     198         END DO 
     199         ! Boundaries conditions 
    179200         CALL lbc_lnk( sshf_a, 'F', 1. ) 
    180201      ENDIF 
     
    183204      !                                           !     Now Vertical Velocity    ! 
    184205      !                                           !------------------------------! 
     206      z1_2dt = 1.e0 / z2dt 
    185207      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    186          wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)     & 
    187             &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
     208         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
     209         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
     210            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
     211            &                         * tmask(:,:,jk) * z1_2dt 
    188212      END DO 
    189213 
     
    248272            sshf_n(:,:) = sshf_a(:,:) 
    249273         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap 
     274            zec = atfp * rdt / rau0 
    250275            DO jj = 1, jpj 
    251276               DO ji = 1, jpi                                ! before <-- now filtered 
    252                   sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) ) 
    253                   sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 
    254                   sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 
    255                   sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 
     277                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )   & 
     278                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 
    256279                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after 
    257280                  sshu_n(ji,jj) = sshu_a(ji,jj) 
     
    260283               END DO 
    261284            END DO 
     285            DO jj = 1, jpjm1 
     286               DO ji = 1, jpim1      ! NO Vector Opt. 
     287                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   & 
     288                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     289                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
     290                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   & 
     291                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     292                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
     293               END DO 
     294            END DO 
     295            ! Boundaries conditions 
     296            CALL lbc_lnk( sshu_b, 'U', 1. ) 
     297            CALL lbc_lnk( sshv_b, 'V', 1. ) 
     298            DO jj = 1, jpjm1 
     299               DO ji = 1, jpim1      ! NO Vector Opt. 
     300                  sshf_b(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
     301                     &                 / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     302                     &                 * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_b(ji,jj  )     & 
     303                     &                   + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) ) 
     304               END DO 
     305            END DO 
     306            ! Boundaries conditions 
     307            CALL lbc_lnk( sshf_b, 'F', 1. ) 
    262308         ENDIF 
    263309         !                    !--------------------------! 
     
    268314            ! 
    269315         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap 
    270             zec = atfp * rdt / rau0 
    271316            DO jj = 1, jpj 
    272317               DO ji = 1, jpi                                ! before <-- now filtered 
    273                   sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   & 
    274                      &                      - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 
     318                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 
    275319                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after 
    276320               END DO 
     
    279323         ! 
    280324      ENDIF 
    281  
    282       ! time filter with global conservation correction and array swap 
    283       sshbb(:,:) = sshb(:,:) 
    284       sshb (:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + zssha(:,:) )    & 
    285            &       - zfact  *  
    286       sshn (:,:) = zssha(:,:) 
    287       empb (:,:) = emp(:,:) 
    288  
    289  
    290325      ! 
    291326      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90

    r1613 r1975  
    217217         hdivb(:,:,:) = hdivn(:,:,:) 
    218218         sshb (:,:)   = sshn (:,:) 
     219         ! - ML - sshbnc  
    219220      ENDIF 
    220221      ! 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r1870 r1975  
    5252   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps   , emps_b   !: freshwater budget: concentration/dillution   [Kg/m2/s] 
    5353   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P-R over ocean and ice               [Kg/m2/s] 
     54   ! - ML - beginning 
     55   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_hc_n      !: sbc heat content trend now                   [K/m/s] 
     56   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_hc_b      !:  "   "      "      "   before                   " 
     57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_sc_n      !: sbc salt content trend now                   [psu/m/s] 
     58   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_trd_sc_b      !:  "   "      "      "   before                   " 
     59   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_trd_hc_n  !: heat content trend due to qsr flux now       [K/m/s] 
     60   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_trd_hc_b  !:  "      "      "    "  "   "   "   before       " 
     61   ! - ML - end 
    5462   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s] 
    5563   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip           !: solid precipitation                          [Kg/m2/s] 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90

    r1870 r1975  
    193193         !                                         ! ---------------------------------------- ! 
    194194         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields 
    195          utau_b(:,:) = utau(:,:)                         ! (except at nitOOO where before fields 
     195         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields 
    196196         qns_b (:,:) = qns (:,:)                         !  are set the end of the routine) 
    197197         qsr_b (:,:) = qsr (:,:) 
    198198         emp_b (:,:) = emp (:,:) 
    199199         emps_b(:,:) = emps(:,:) 
     200         ! - ML - 
     201         sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:) 
     202         qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:) 
     203         IF ( .NOT. lk_vvl )  sbc_trd_sc_b(:,:) = sbc_trd_sc_n(:,:) 
     204          
    200205      ENDIF 
    201206 
     
    256261            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file' 
    257262            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point) 
    258             CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before j-stress  (V-point) 
     263            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point) 
    259264            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point) 
    260265            CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point) 
    261266            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point) 
    262267            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emp_b  )   ! before C/D freshwater flux (T-point) 
     268            ! - ML - 
     269            CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b )   ! before heat content sbc trend 
     270            CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b )   ! before heat content trend due to Qsr flux 
     271            IF ( .NOT. lk_vvl ) THEN 
     272               CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b )   ! before salt content sbc trend 
     273            ENDIF 
    263274            ! 
    264275         ELSE                                                   !* no restart: set from nit000 values 
    265276            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000' 
    266277            utau_b(:,:) = utau(:,:)  
    267             utau_b(:,:) = utau(:,:) 
     278            vtau_b(:,:) = vtau(:,:) 
    268279            qns_b (:,:) = qns (:,:) 
    269280            qsr_b (:,:) = qsr (:,:) 
     
    280291            &                    'at it= ', kt,' date= ', ndastp 
    281292         IF(lwp) WRITE(numout,*) '~~~~' 
    282          CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )    !  
    283          CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , vtau ) 
    284          CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  ) 
    285          CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  ) 
    286          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  ) 
    287          CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emp  ) 
     293         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau_b )    !  
     294         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau_b ) 
     295         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns_b  ) 
     296         CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr_b  ) 
     297         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp_b  ) 
     298         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps_b ) 
     299         ! - ML - 
     300         CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_b ) 
     301         CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_b ) 
     302         IF ( .NOT. lk_vvl ) THEN 
     303            CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_b ) 
     304         ENDIF 
    288305         ! 
    289306      ENDIF 
  • 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' ) 
  • branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90

    r1793 r1975  
    116116   USE restart         ! ocean restart                    (rst_wri routine) 
    117117   USE prtctl          ! Print control                    (prt_ctl routine) 
     118   ! - ML -  
     119   USE diatrb          ! global tracer conservation       (dia_trb  routine) 
    118120 
    119121#if defined key_agrif 
     
    317319 
    318320                               CALL ssh_nxt( kstp )         ! sea surface height at next time step 
     321                               CALL dia_trb( kstp )         ! - ML - global conservation diagnostics 
    319322 
    320323      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.