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 2257 for branches/nemo_v3_3_beta/NEMOGCM/NEMO – NEMO

Ignore:
Timestamp:
2010-10-13T17:58:28+02:00 (14 years ago)
Author:
cetlod
Message:

Apply the modified leap-frog scheme on runoff & correct a bug on time stepping for hpg implicit case

Location:
branches/nemo_v3_3_beta/NEMOGCM/NEMO
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r2236 r2257  
    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 over ocean and ice                 [Kg/m2/s] 
    54    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rnf               !: river runoff   [Kg/m2/s]   
     54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rnf    , rnf_b    !: river runoff   [Kg/m2/s]   
    5555   ! - ML - begin 
    56    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_n          !: sbc heat content trend now                   [K.m/s] 
    57    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_hc_b          !:  "   "      "      "   before                   " 
    58    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_n          !: sbc salt content trend now                   [psu.m/s] 
    59    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sbc_sc_b          !:  "   "      "      "   before                   " 
    60    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_n      !: heat content trend due to qsr flux now       [K.m/s] 
    61    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc_b      !:  "      "      "    "  "   "   "   before       " 
     56   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] 
     57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc , qsr_hc_b   !: heat content trend due to qsr flux     [K.m/s] 
    6258   ! - ML - end 
    6359   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s] 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r2239 r2257  
    3131   CHARACTER(len=100), PUBLIC ::   cn_dir       = './'    !: Root directory for location of ssr files 
    3232   LOGICAL           , PUBLIC ::   ln_rnf_depth = .false. !: depth       river runoffs attribute specified in a file 
    33    LOGICAL           , PUBLIC ::   ln_rnf_temp  = .false. !: temperature river runoffs attribute specified in a file  
     33   LOGICAL           , PUBLIC ::   ln_rnf_tem   = .false. !: temperature river runoffs attribute specified in a file  
    3434   LOGICAL           , PUBLIC ::   ln_rnf_sal   = .false. !: salinity    river runoffs attribute specified in a file  
    3535   LOGICAL           , PUBLIC ::   ln_rnf_emp   = .false. !: runoffs into a file to be read or already into precipitation 
     
    5656   INTEGER,  PUBLIC, DIMENSION(jpi,jpj) ::   nk_rnf       !: depth of runoff in model levels 
    5757 
    58    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2) ::  tsc_rnf  !: temperature & salinity content of river runoffs   [K.m/s & PSU.m/s] 
     58   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc  !: before and now  
     59   !                                                                 !: temp. & sal. content of river runoffs   [K.m/s & PSU.m/s] 
    5960 
    6061   !! * Substitutions   
     
    8889      IF( kt == nit000 )   CALL sbc_rnf_init                           ! Read namelist and allocate structures 
    8990 
     91      !                                            ! ---------------------------------------- ! 
     92      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          ! 
     93         !                                         ! ---------------------------------------- ! 
     94         rnf_b    (:,:  ) = rnf    (:,:  )               ! Swap the ocean forcing fields except at nit000 
     95         rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)               ! where before fields are set at the end of the routine 
     96         ! 
     97      ENDIF 
     98 
    9099      !                                                   !-------------------! 
    91100      IF( .NOT. ln_rnf_emp ) THEN                         !   Update runoff   ! 
    92101         !                                                !-------------------! 
    93102         ! 
    94                              CALL fld_read ( kt, nn_fsbc, sf_rnf     ! Read Runoffs data and provide it at kt  
    95          IF( ln_rnf_temp )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required 
     103                             CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt  
     104         IF( ln_rnf_tem  )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required 
    96105         IF( ln_rnf_sal  )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required 
    97106 
     
    107116 
    108117         ! C a u t i o n : runoff is negative and in kg/m2/s  
    109  
    110          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
     118         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    111119            rnf(:,:)  = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )   
    112120            ! 
    113121            z1_rau0 = 1.e0 / rau0 
    114122            !                                                              ! set temperature & salinity content of runoffs 
    115             IF( ln_rnf_temp )   THEN                                       ! use runoffs temperature data 
    116                tsc_rnf(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0 
     123            IF( ln_rnf_tem )   THEN                                        ! use runoffs temperature data 
     124               rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0 
    117125               WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 )                      ! if missing data value use SST as runoffs temperature   
    118                    tsc_rnf(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
     126                   rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
    119127               ENDWHERE 
    120128            ELSE                                                           ! use SST as runoffs temperature 
    121                tsc_rnf(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
     129               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
    122130            ENDIF   
    123131            !                                                              ! use runoffs salinity data  
    124             IF( ln_rnf_sal ) tsc_rnf(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0 
     132            IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0 
    125133            !                                                              ! else use S=0 for runoffs (done one for all in the init) 
    126134            ! 
    127             IF( ln_rnf_temp .OR. ln_rnf_sal ) THEN                         ! runoffs as outflow: use ocean SST and SSS 
     135            IF( ln_rnf_tem .OR. ln_rnf_sal ) THEN                         ! runoffs as outflow: use ocean SST and SSS 
    128136               WHERE( rnf(:,:) < 0.e0 )                                    ! example baltic model when flow is out of domain  
    129                   tsc_rnf(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
    130                   tsc_rnf(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * z1_rau0 
     137                  rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
     138                  rnf_tsc(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * z1_rau0 
    131139               ENDWHERE 
    132140            ENDIF 
    133  
     141            ! 
    134142            CALL iom_put( "runoffs", rnf )         ! output runoffs arrays 
    135143         ENDIF 
    136144         ! 
    137145      ENDIF 
     146      ! 
     147      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
     148         !                                             ! ---------------------------------------- ! 
     149         IF( ln_rstart .AND.    &                               !* Restart: read in restart file 
     150            & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN  
     151            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields red in the restart file' 
     152            CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b )     ! before runoff 
     153            CALL iom_get( numror, jpdom_autoglo, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) )   ! before heat content of runoff 
     154            CALL iom_get( numror, jpdom_autoglo, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) )   ! before salinity content of runoff 
     155         ELSE                                                   !* no restart: set from nit000 values 
     156            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000' 
     157             rnf_b    (:,:  ) = rnf    (:,:  )   
     158             rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)    
     159         ENDIF 
     160      ENDIF 
     161      !                                                ! ---------------------------------------- ! 
     162      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     ! 
     163         !                                             ! ---------------------------------------- ! 
     164         IF(lwp) WRITE(numout,*) 
     165         IF(lwp) WRITE(numout,*) 'sbcrnf : runoff forcing fields written in ocean restart file ',   & 
     166            &                    'at it= ', kt,' date= ', ndastp 
     167         IF(lwp) WRITE(numout,*) '~~~~' 
     168         CALL iom_rstput( kt, nitrst, numrow, 'rnf_b' , rnf ) 
     169         CALL iom_rstput( kt, nitrst, numrow, 'rnf_hc_b', rnf_tsc(:,:,jp_tem) ) 
     170         CALL iom_rstput( kt, nitrst, numrow, 'rnf_sc_b', rnf_tsc(:,:,jp_sal) ) 
     171      ENDIF 
     172 
    138173      ! 
    139174   END SUBROUTINE sbc_rnf 
     
    202237      INTEGER           ::   ierror, inum  ! temporary integer 
    203238      !!  
    204       NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, ln_rnf_depth, ln_rnf_temp, ln_rnf_sal,   & 
    205          &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf   , sn_dep_rnf,   &   
    206          &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf , rn_rfact   
     239      NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   & 
     240         &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &   
     241         &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact   
    207242      !!---------------------------------------------------------------------- 
    208243 
     
    243278         IF(lwp) WRITE(numout,*) 
    244279         IF(lwp) WRITE(numout,*) '          runoffs directly provided in the precipitations' 
    245          IF( ln_rnf_depth .OR. ln_rnf_temp .OR. ln_rnf_sal ) THEN 
     280         IF( ln_rnf_depth .OR. ln_rnf_tem .OR. ln_rnf_sal ) THEN 
    246281           CALL ctl_warn( 'runoffs already included in precipitations, so runoff (T,S, depth) attributes will not be used' )  
    247            ln_rnf_depth = .FALSE.   ;   ln_rnf_temp = .FALSE.   ;   ln_rnf_sal = .FALSE. 
     282           ln_rnf_depth = .FALSE.   ;   ln_rnf_tem = .FALSE.   ;   ln_rnf_sal = .FALSE. 
    248283         ENDIF 
    249284         ! 
     
    261296         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' ) 
    262297         ! 
    263          IF( ln_rnf_temp ) THEN                     ! Create (if required) sf_t_rnf structure 
     298         IF( ln_rnf_tem ) THEN                     ! Create (if required) sf_t_rnf structure 
    264299            IF(lwp) WRITE(numout,*) 
    265300            IF(lwp) WRITE(numout,*) '          runoffs temperatures read in a file' 
     
    294329            CALL iom_close( inum )                                      ! close file   
    295330   
    296             nk_rnf(:,:)=0                              ! set the number of level over which river runoffs are applied 
    297             DO jj=1,jpj   
    298               DO ji=1,jpi   
     331            nk_rnf(:,:) = 0                              ! set the number of level over which river runoffs are applied 
     332            DO jj = 1, jpj   
     333              DO ji = 1, jpi   
    299334                IF ( h_rnf(ji,jj) > 0.e0 ) THEN   
    300                   jk= 
    301                   DO WHILE ( jk/=(mbathy(ji,jj)-1) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) );  jk=jk+1;   ENDDO   
    302                   nk_rnf(ji,jj)=jk   
    303                 ELSE IF ( h_rnf(ji,jj) == -1   ) THEN   ;  nk_rnf(ji,jj)= 
    304                 ELSE IF ( h_rnf(ji,jj) == -999 ) THEN   ;  nk_rnf(ji,jj)=mbathy(ji,jj)-1 
     335                  jk =  
     336                  DO WHILE ( jk /= ( mbathy(ji,jj) - 1 ) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 ;  ENDDO   
     337                  nk_rnf(ji,jj) = jk   
     338                ELSE IF ( h_rnf(ji,jj) == -1   ) THEN   ;  nk_rnf(ji,jj) =  
     339                ELSE IF ( h_rnf(ji,jj) == -999 ) THEN   ;  nk_rnf(ji,jj) = mbathy(ji,jj) - 1 
    305340                ELSE IF ( h_rnf(ji,jj) /= 0 ) THEN   
    306341                  CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )   
     
    309344              ENDDO   
    310345            ENDDO   
    311             DO jj=1,jpj                               ! set the associated depth  
    312               DO ji=1,jpi  
    313                 h_rnf(ji,jj)=0.e0 
    314                 DO jk=1,nk_rnf(ji,jj)                         
    315                    h_rnf(ji,jj)=h_rnf(ji,jj)+fse3t(ji,jj,jk)   
     346            DO jj = 1, jpj                               ! set the associated depth  
     347              DO ji = 1, jpi  
     348                h_rnf(ji,jj) = 0.e0 
     349                DO jk = 1, nk_rnf(ji,jj)                         
     350                   h_rnf(ji,jj) = h_rnf(ji,jj)+fse3t(ji,jj,jk)   
    316351                ENDDO 
    317352              ENDDO 
    318353            ENDDO 
    319354         ELSE                                       ! runoffs applied at the surface  
    320             nk_rnf(:,:)= 
    321             h_rnf(:,:)=fse3t(:,:,1) 
     355            nk_rnf(:,:) =  
     356            h_rnf(:,:)  = fse3t(:,:,1) 
    322357         ENDIF   
    323358      !  
    324359      ENDIF 
    325360 
    326       tsc_rnf(:,:,:) = 0.e0                 ! runoffs temperature & salinty contents initilisation 
     361      rnf_tsc(:,:,:) = 0.e0                 ! runoffs temperature & salinty contents initilisation 
    327362      !                                   ! ======================== 
    328363      !                                   !   River mouth vicinity 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r2240 r2257  
    5757   REAL(wp)                 ::   rbcp            ! Brown & Campana parameters for semi-implicit hpg 
    5858   REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 
     59   LOGICAL  :: l_tra           ! active tracers or passive tracers 
    5960 
    6061   !! * Substitutions 
     
    9495      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    9596      !! 
    96       INTEGER  ::   jk    ! dummy loop indices 
    97       REAL(wp) ::   zfact ! local scalars 
     97      INTEGER  ::   jk, jn    ! dummy loop indices 
     98      REAL(wp) ::   zfact     ! local scalars 
    9899      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds 
    99100      !!---------------------------------------------------------------------- 
     
    142143 
    143144      ! Leap-Frog + Asselin filter time stepping 
    144       IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)      
    145       ELSE                  ;   CALL tra_nxt_fix( kt, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level  
    146       ENDIF 
     145      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     146         !                                             ! (only swap) 
     147         DO jn = 1, jpts 
     148            DO jk = 1, jpkm1 
     149               tsn(:,:,jk,jn) = tsa(:,:,jk,jn)     
     150            END DO 
     151         END DO 
     152         !                                               
     153      ELSE 
     154         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)      
     155         ELSE                 ;   CALL tra_nxt_fix( kt, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level  
     156         ENDIF 
     157      ENDIF  
    147158 
    148159#if defined key_agrif 
     
    202213      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend 
    203214      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices 
    204       REAL(wp) :: ztd, ztm         ! temporary scalars 
     215      REAL(wp) :: ztn, ztd, ztm         ! temporary scalars 
    205216      !!---------------------------------------------------------------------- 
    206217 
     
    211222      ENDIF 
    212223      ! 
    213       ! 
    214       IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
    215          !                                             ! (only swap) 
    216          DO jn = 1, kjpt 
    217             DO jk = 1, jpkm1 
    218                ptn(:,:,jk,jn) = pta(:,:,jk,jn)     ! ptb <-- ptn 
    219             END DO 
     224      IF( cdtype == 'TRA' )  THEN   ;   l_tra     = .TRUE.    ! active tracers case      
     225      ELSE                          ;   l_tra     = .FALSE.   ! passive tracers case 
     226      ENDIF 
     227      ! 
     228      DO jn = 1, kjpt 
     229         ! 
     230         DO jk = 1, jpkm1 
     231            DO jj = 1, jpj 
     232               DO ji = 1, jpi 
     233                  IF( l_tra .AND. ln_dynhpg_imp )  ztn = ptn(ji,jj,jk,jn)               ! implicit hpg: keep tn, sn in memory 
     234                  !  
     235                  ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)      !  time laplacian on tracers 
     236                  ! 
     237                  ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                      ! ptb <-- filtered ptn  
     238                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                   ! ptn <-- pta 
     239                  ! 
     240                  IF( l_tra .AND. ln_dynhpg_imp )  pta(ji,jj,jk,jn) = ztn + rbcp * ztd  ! pta <-- Brown & Campana average 
     241               END DO 
     242           END DO 
    220243         END DO 
    221          !                                               
    222       ELSE                                           ! general case (Leapfrog + Asselin filter 
    223244         ! 
    224          !                                                     ! ----------------------- ! 
    225          IF( ln_dynhpg_imp .AND. cdtype == 'TRA' ) THEN        ! semi-implicite hpg case ! 
    226             !                                                  ! ----------------------- ! 
    227             DO jn = 1, kjpt 
    228                DO jk = 1, jpkm1 
    229                   DO jj = 1, jpj 
    230                      DO ji = 1, jpi 
    231                         ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    !  time laplacian on tracers 
    232                         ztm = ptn(ji,jj,jk,jn) + rbcp * ztd                                 !  used for both Asselin and Brown & Campana filters 
    233                         ! 
    234                         ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn  
    235                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta 
    236                         pta(ji,jj,jk,jn) = ztm                                              ! pta <-- Brown & Campana average 
    237                      END DO 
    238                   END DO 
    239                END DO 
    240             END DO 
    241             !                                           ! ----------------------- ! 
    242          ELSE                                           !    explicit hpg case    ! 
    243             !                                           ! ----------------------- ! 
    244             DO jn = 1, kjpt 
    245                DO jk = 1, jpkm1 
    246                   DO jj = 1, jpj 
    247                      DO ji = 1, jpi 
    248                         ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn)    ! time laplacian on tracers 
    249                         !  
    250                         ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd                    ! ptb <-- filtered ptn  
    251                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                                 ! ptn <-- pta 
    252                      END DO 
    253                   END DO 
    254                END DO 
    255             END DO 
    256          ENDIF 
    257          ! 
    258       ENDIF 
     245      END DO 
    259246      ! 
    260247   END SUBROUTINE tra_nxt_fix 
     
    306293      ENDIF 
    307294      ! 
    308       ! 
    309       IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step 
    310          !                                           ! (only swap) 
    311          DO jn = 1, kjpt 
    312             DO jk = 1, jpkm1 
    313                ptn(:,:,jk,jn) = pta(:,:,jk,jn)       ! ptb <-- ptn 
     295      IF( cdtype == 'TRA' )  THEN   ;   l_tra     = .TRUE.    ! active tracers case      
     296      ELSE                          ;   l_tra     = .FALSE.   ! passive tracers case 
     297      ENDIF 
     298      ! 
     299      DO jn = 1, kjpt       
     300         DO jk = 1, jpkm1 
     301            zfact1 = atfp * rdttra(jk) 
     302            zfact2 = zfact1 / rau0 
     303            DO jj = 1, jpj 
     304               DO ji = 1, jpi 
     305                  ze3t_b = fse3t_b(ji,jj,jk) 
     306                  ze3t_n = fse3t_n(ji,jj,jk) 
     307                  ze3t_a = fse3t_a(ji,jj,jk) 
     308                  !                                         ! tracer content at Before, now and after 
     309                  ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
     310                  ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
     311                  ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
     312                  ! 
     313                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 
     314                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b 
     315                  ! 
     316                  ze3t_f = ze3t_n + atfp * ze3t_d 
     317                  ztc_f  = ztc_n  + atfp * ztc_d 
     318 
     319                  IF( l_tra .AND. jk == 1 ) THEN 
     320                      ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 
     321                      ztc_f  = ztc_f  - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) ) 
     322                  ENDIF 
     323                  IF( l_tra .AND. jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr )  & 
     324                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
     325 
     326                   ze3t_f = 1.e0 / ze3t_f 
     327                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
     328                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
     329                   ! 
     330                   IF( l_tra .AND. ln_dynhpg_imp ) THEN 
     331                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
     332                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
     333                   ENDIF 
     334               END DO 
    314335            END DO 
    315336         END DO 
    316          !                                               
    317       ELSE                                           ! general case (Leapfrog + Asselin filter) 
    318          ! 
    319          !                                                     ! ----------------------- ! 
    320          IF( ln_dynhpg_imp .AND. cdtype == 'TRA' ) THEN        ! semi-implicite hpg case ! 
    321             !                                                  ! ----------------------- ! 
    322             DO jn = 1, kjpt                           
    323                DO jk = 1, jpkm1 
    324                   DO jj = 1, jpj 
    325                      DO ji = 1, jpi 
    326                         ze3t_b = fse3t_b(ji,jj,jk) 
    327                         ze3t_n = fse3t_n(ji,jj,jk) 
    328                         ze3t_a = fse3t_a(ji,jj,jk) 
    329                         !                                         ! tracer content at Before, now and after 
    330                         ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b   
    331                         ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n   
    332                         ztc_a  = pta(ji,jj,jk,jn) * ze3t_a   
    333                         !                                         ! Time laplacian on tracer contents 
    334                         !                                         ! used for both Asselin and Brown & Campana filters 
    335                         ze3t_d =  ze3t_b - 2.* ze3t_n + ze3t_a  
    336                         ztc_d  =  ztc_b  - 2.* ztc_n  + ztc_a   
    337                         !                                         ! Asselin Filter on thicknesses and tracer contents 
    338                         ztc_f  = ztc_n  + atfp * ztc_d 
    339                         ztc_m  = ztc_n  + rbcp * ztc_d 
    340                         !                                          
    341                         ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d ) 
    342                         ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d ) 
    343                         !                                         ! swap of arrays 
    344                         ptb(ji,jj,jk,jn) = ztc_f * ze3t_f              ! ptb <-- ptn filtered 
    345                         pta(ji,jj,jk,jn) = ztc_m * ze3t_m              ! pta <-- Brown & Campana average 
    346                         ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)            ! ptn <-- pta 
    347                      END DO 
    348                   END DO 
    349                END DO 
    350             END DO 
    351             !                                           ! ----------------------- ! 
    352          ELSE                                           !    explicit hpg case    ! 
    353             !                                           ! ----------------------- ! 
    354             IF( cdtype == 'TRA' ) THEN 
    355                !               
    356                DO jn = 1, kjpt       
    357                   DO jk = 1, jpkm1 
    358                      zfact1 = atfp * rdttra(jk) 
    359                      zfact2 = zfact1 / rau0 
    360                      DO jj = 1, jpj 
    361                         DO ji = 1, jpi 
    362                            ze3t_b = fse3t_b(ji,jj,jk) 
    363                            ze3t_n = fse3t_n(ji,jj,jk) 
    364                            ze3t_a = fse3t_a(ji,jj,jk) 
    365                            !                                         ! tracer content at Before, now and after 
    366                            ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
    367                            ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
    368                            ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
    369                            ! 
    370                            ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 
    371                            ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b ) 
    372  
    373                            IF( jk == 1 ) THEN 
    374                                ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 
    375                                IF( jn == jp_tem )   ztc_f  = ztc_f  - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) ) 
    376                                IF( jn == jp_sal )   ztc_f  = ztc_f  - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) ) 
    377                            ENDIF 
    378                            IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr )  & 
    379                            &                        ztc_f  = ztc_f  - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
    380  
    381                            ze3t_f = 1.e0 / ze3t_f 
    382                            ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered 
    383                            ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta 
    384                         END DO 
    385                      END DO 
    386                   END DO 
    387                END DO 
    388                ! 
    389             ELSE IF( cdtype == 'TRC' ) THEN 
    390                !               
    391                DO jn = 1, kjpt       
    392                   DO jk = 1, jpkm1 
    393                      DO jj = 1, jpj 
    394                         DO ji = 1, jpi 
    395                            ze3t_b = fse3t_b(ji,jj,jk) 
    396                            ze3t_n = fse3t_n(ji,jj,jk) 
    397                            ze3t_a = fse3t_a(ji,jj,jk) 
    398                            !                                         ! tracer content at Before, now and after 
    399                            ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b 
    400                            ztc_n  = ptn(ji,jj,jk,jn) * ze3t_n 
    401                            ztc_a  = pta(ji,jj,jk,jn) * ze3t_a 
    402                            ! 
    403                            ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 
    404                            ztc_f  = ztc_n  + atfp * ( ztc_a  - 2. * ztc_n  + ztc_b  ) 
    405  
    406                            ze3t_f = 1.e0 / ze3t_f 
    407                            ptb(ji,jj,jk,jn) = ztc_f * ze3t_f                           ! tb <-- tn filtered 
    408                            ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                             ! tn <-- ta 
    409                         END DO 
    410                      END DO 
    411                   END DO 
    412                END DO 
    413                ! 
    414             END IF 
    415             ! 
    416          ENDIF 
    417          ! 
    418       ENDIF 
     337         !  
     338      END DO 
    419339      ! 
    420340   END SUBROUTINE tra_nxt_vvl 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r2236 r2257  
    132132         !                                        ! --------------------- 
    133133         zfact = 0.5e0 
    134          qsr_hc_b(:,:,:) = qsr_hc_n(:,:,:) 
     134         qsr_hc_b(:,:,:) = qsr_hc(:,:,:) 
    135135      ENDIF 
    136136      !                                        Compute now qsr tracer content field 
     
    143143            DO jj = 2, jpjm1 
    144144               DO ji = fs_2, fs_jpim1   ! vector opt. 
    145                   qsr_hc_n(ji,jj,jk) =  ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
     145                  qsr_hc(ji,jj,jk) =  ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
    146146               END DO 
    147147            END DO 
     
    200200               ! 
    201201               DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
    202                   qsr_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
     202                  qsr_hc(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
    203203               END DO 
    204204               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero 
     
    207207            ELSE                                                 !*  Constant Chlorophyll 
    208208               DO jk = 1, nksr 
    209                   qsr_hc_n(:,:,jk) =  etot3(:,:,jk) * qsr(:,:) 
     209                  qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:) 
    210210               END DO 
    211211            ENDIF 
     
    219219               DO jj = 2, jpjm1 
    220220                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    221                      qsr_hc_n(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) 
     221                     qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) 
    222222                  END DO 
    223223               END DO 
     
    233233            DO ji = fs_2, fs_jpim1   ! vector opt. 
    234234               z1_e3t = zfact / fse3t(ji,jj,jk) 
    235                tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc_n(ji,jj,jk) ) * z1_e3t 
     235               tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 
    236236            END DO 
    237237         END DO 
     
    244244            &                    'at it= ', kt,' date= ', ndastp 
    245245         IF(lwp) WRITE(numout,*) '~~~~' 
    246          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc_n ) 
     246         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc ) 
    247247         ! 
    248248      ENDIF 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r2239 r2257  
    139139            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file' 
    140140            zfact = 0.5e0 
    141             CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_hc_b )   ! before heat content sbc trend 
    142             CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_sc_b )   ! before salt content sbc trend 
     141            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend 
     142            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend 
    143143         ELSE                                         ! No restart or restart not found: Euler forward time stepping 
    144144            zfact = 1.e0 
    145             sbc_hc_b(:,:) = 0.e0 
    146             sbc_sc_b(:,:) = 0.e0 
     145            sbc_tsc_b(:,:,:) = 0.e0 
    147146         ENDIF 
    148147      ELSE                                         ! Swap of forcing fields 
    149148         !                                         ! ---------------------- 
    150149         zfact = 0.5e0 
    151          sbc_hc_b(:,:) = sbc_hc_n(:,:) 
    152          sbc_sc_b(:,:) = sbc_sc_n(:,:) 
     150         sbc_tsc_b(:,:,:) = sbc_tsc(:,:) 
    153151      ENDIF 
    154152      !                                          Compute now sbc tracer content fields 
     
    162160            DO ji = fs_2, fs_jpim1   ! vector opt. 
    163161               ! temperature : heat flux + cooling/heating effet of EMP flux 
    164                sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 
     162               sbc_tsc(ji,jj,jp_tem) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 
    165163               ! concent./dilut. effect due to sea-ice melt/formation and (possibly) SSS restoration 
    166                sbc_sc_n(ji,jj) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) 
     164               sbc_tsc(ji,jj,jp_sal) = ( emps(ji,jj) - emp(ji,jj) ) * zsrau * tsn(ji,jj,1,jp_sal) 
    167165            END DO 
    168166         END DO 
     
    171169            DO ji = fs_2, fs_jpim1   ! vector opt. 
    172170               ! temperature : heat flux 
    173                sbc_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 
     171               sbc_tsc(ji,jj,jp_tem) = ro0cpr * qns(ji,jj) 
    174172               ! salinity    : salt flux + concent./dilut. effect (both in emps) 
    175                sbc_sc_n(ji,jj) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal) 
     173               sbc_tsc(ji,jj,jp_sal) = zsrau * emps(ji,jj) * tsn(ji,jj,1,jp_sal) 
    176174            END DO 
    177175         END DO 
    178176      ENDIF 
    179177      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff   
    180       DO jj = 2, jpj 
    181          DO ji = fs_2, fs_jpim1   ! vector opt. 
    182             z1_e3t = zfact / fse3t(ji,jj,1) 
    183             tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + ( sbc_hc_b(ji,jj) + sbc_hc_n(ji,jj) ) * z1_e3t 
    184             tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + ( sbc_sc_b(ji,jj) + sbc_sc_n(ji,jj) ) * z1_e3t 
     178      DO jn = 1, jpts 
     179         DO jj = 2, jpj 
     180            DO ji = fs_2, fs_jpim1   ! vector opt. 
     181               z1_e3t = zfact / fse3t(ji,jj,1) 
     182               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t 
     183            END DO 
    185184         END DO 
    186185      END DO 
     
    192191            &                    'at it= ', kt,' date= ', ndastp 
    193192         IF(lwp) WRITE(numout,*) '~~~~' 
    194          CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_hc_n ) 
    195          CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_sc_n ) 
     193         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) 
     194         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) ) 
    196195      ENDIF 
    197196      !                             !==  Runoffs  ==! 
     
    200199         DO jj = 2, jpj  
    201200            DO ji = fs_2, fs_jpim1 
    202                zdep = 1. / h_rnf(ji,jj)   
     201               zdep = 1. / h_rnf(ji,jj) 
     202               zdep = zfact * zdep   
    203203               IF ( rnf(ji,jj) .ne. 0.0 ) THEN 
    204204                  DO jk = 1, nk_rnf(ji,jj) 
    205                                         tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + tsc_rnf(ji,jj,jp_tem) * zdep 
    206                      IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) + tsc_rnf(ji,jj,jp_sal) * zdep 
     205                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     206                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 
     207                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
     208                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep  
    207209                  ENDDO 
    208210               ENDIF 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90

    r2240 r2257  
    120120         ztrdt(:,:,:,:)  = trn(:,:,:,:) 
    121121      ENDIF 
    122  
    123122      ! Leap-Frog + Asselin filter time stepping 
    124       IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, 'TRC', trb, trn, tra, jptra )      ! variable volume level (vvl)  
    125       ELSE                ;   CALL tra_nxt_fix( kt, 'TRC', trb, trn, tra, jptra )      ! fixed    volume level  
     123      IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step 
     124         !                                             ! (only swap) 
     125         DO jn = 1, jptra 
     126            DO jk = 1, jpkm1 
     127               trn(:,:,jk,jn) = tra(:,:,jk,jn) 
     128            END DO 
     129         END DO 
     130         !                                               
     131      ELSE 
     132         ! Leap-Frog + Asselin filter time stepping 
     133         IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, 'TRC', trb, trn, tra, jptra )      ! variable volume level (vvl)  
     134         ELSE                ;   CALL tra_nxt_fix( kt, 'TRC', trb, trn, tra, jptra )      ! fixed    volume level  
     135         ENDIF 
    126136      ENDIF 
    127137 
Note: See TracChangeset for help on using the changeset viewer.