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 5901 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 – NEMO

Ignore:
Timestamp:
2015-11-20T09:39:06+01:00 (8 years ago)
Author:
jamesharle
Message:

merging branch with head of the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r5038 r5901  
    2626   USE trc_oce         ! share ocean/biogeo variables 
    2727   USE phycst          ! physical constants 
     28   USE ldftra          ! lateral diffusivity coefficients 
    2829   USE trabbl          ! active tracer: bottom boundary layer 
    2930   USE ldfslp          ! lateral diffusion: iso-neutral slopes 
    30    USE ldfeiv          ! eddy induced velocity coef.  
    31    USE ldftra_oce      ! ocean tracer   lateral physics 
    3231   USE zdfmxl          ! vertical physics: mixed layer depth 
    3332   USE eosbn2          ! equation of state - Brunt Vaisala frequency 
     
    4039   USE fldread         ! read input fields  
    4140   USE timing          ! Timing 
     41   USE wrk_nemo 
    4242 
    4343   IMPLICIT NONE 
     
    5050   LOGICAL            ::   ln_dynwzv    !: vertical velocity read in a file (T) or computed from u/v (F) 
    5151   LOGICAL            ::   ln_dynbbl    !: bbl coef read in a file (T) or computed (F) 
    52    LOGICAL            ::   ln_degrad    !: degradation option enabled or not 
    5352   LOGICAL            ::   ln_dynrnf    !: read runoff data in file (T) or set to zero (F) 
    5453 
    55    INTEGER  , PARAMETER ::   jpfld = 21     ! maximum number of fields to read 
     54   INTEGER  , PARAMETER ::   jpfld = 15     ! maximum number of fields to read 
    5655   INTEGER  , SAVE      ::   jf_tem         ! index of temperature 
    5756   INTEGER  , SAVE      ::   jf_sal         ! index of salinity 
     
    6867   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef 
    6968   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef 
    70    INTEGER  , SAVE      ::   jf_ahu         ! index of u-diffusivity coef 
    71    INTEGER  , SAVE      ::   jf_ahv         ! index of v-diffusivity coef  
    72    INTEGER  , SAVE      ::   jf_ahw         ! index of w-diffusivity coef 
    73    INTEGER  , SAVE      ::   jf_eiu         ! index of u-eiv 
    74    INTEGER  , SAVE      ::   jf_eiv         ! index of v-eiv 
    75    INTEGER  , SAVE      ::   jf_eiw         ! index of w-eiv 
    7669   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux 
    7770 
     
    112105      !!             - interpolates data if needed 
    113106      !!---------------------------------------------------------------------- 
    114       ! 
    115       USE oce, ONLY:  zts    => tsa  
     107      USE oce, ONLY:  zts    => tsa 
    116108      USE oce, ONLY:  zuslp  => ua   , zvslp  => va 
    117       USE oce, ONLY:  zwslpi => rotb , zwslpj => rotn 
    118       USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => hdivb 
     109      USE oce, ONLY:  zwslpi => ua_sv , zwslpj => va_sv 
     110      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => rke 
    119111      ! 
    120112      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     113      ! 
     114!      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)  :: zts 
     115!      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zuslp, zvslp, zwslpi, zwslpj 
     116!      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zu, zv, zw 
     117      ! 
    121118      ! 
    122119      INTEGER  ::   ji, jj     ! dummy loop indices 
     
    138135         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==! 
    139136         ! 
    140          IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
     137         IF( l_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
    141138            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
    142139            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
     
    162159      ENDIF 
    163160      !  
    164       IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                        
     161      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                        
    165162         iswap_tem = 0 
    166163         IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1 
     
    245242      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    246243      ! 
    247       CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) )                                       ! In any case, we need rhop 
     244      ! 
     245      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     246      CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points 
     247      CALL bn2    ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl 
     248 
     249      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    248250      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl  
    249251      ! 
     
    259261      fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction 
    260262      qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation 
    261       IF ( ln_dynrnf ) & 
     263      IF( ln_dynrnf ) & 
    262264      rnf (:,:)        = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! river runoffs  
    263265 
     266      !                                               ! update eddy diffusivity coeff. and/or eiv coeff. at kt 
     267      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kt )  
    264268      !                                                      ! bbl diffusive coef 
    265269#if defined key_trabbl && ! defined key_c1d 
     
    271275         CALL bbl( kt, nit000, 'TRC') 
    272276      END IF 
    273 #endif 
    274 #if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d  
    275       aeiw(:,:)        = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1)    ! w-eiv 
    276       !                                                           ! Computes the horizontal values from the vertical value 
    277       DO jj = 2, jpjm1 
    278          DO ji = fs_2, fs_jpim1   ! vector opt. 
    279             aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )  ! Average the diffusive coefficient at u- v- points 
    280             aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )  ! at u- v- points 
    281          END DO 
    282       END DO 
    283       CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition 
    284 #endif 
    285        
    286 #if defined key_degrad && ! defined key_c1d  
    287       !                                          ! degrad option : diffusive and eiv coef are 3D 
    288       ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:) 
    289       ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:) 
    290       ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:) 
    291 #  if defined key_traldf_eiv  
    292       aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:) 
    293       aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:) 
    294       aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:) 
    295 #  endif 
    296277#endif 
    297278      ! 
     
    334315      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read 
    335316      TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf  ! informations about the fields to be read 
    336       TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl          !   "                                 " 
    337       TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf  !   "                                 " 
    338       !!---------------------------------------------------------------------- 
    339       ! 
    340       NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad, ln_dynrnf,    & 
     317      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf          !   "                                 " 
     318      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_dynrnf,    & 
    341319         &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf,  & 
    342          &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl,          & 
    343          &                sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf 
     320         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf   
     321      !!---------------------------------------------------------------------- 
    344322      ! 
    345323      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
     
    360338         WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv 
    361339         WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl 
    362          WRITE(numout,*) '      degradation option enabled (T) or not (F)            ln_degrad  = ', ln_degrad 
    363340         WRITE(numout,*) '      river runoff option enabled (T) or not (F)           ln_dynrnf  = ', ln_dynrnf 
    364341         WRITE(numout,*) 
    365342      ENDIF 
    366343      !  
    367       IF( ln_degrad .AND. .NOT.lk_degrad ) THEN 
    368          CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' ) 
    369          ln_degrad = .FALSE. 
    370       ENDIF 
    371344      IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN 
    372345         CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' ) 
     
    383356 
    384357      ! 
    385       IF ( ln_dynrnf ) THEN 
     358      IF( ln_dynrnf ) THEN 
    386359                jf_rnf = jfld + 1  ;  jfld  = jf_rnf 
    387360         slf_d(jf_rnf) = sn_rnf 
     
    390363      ENDIF 
    391364 
    392       ! 
    393       IF( .NOT.ln_degrad ) THEN     ! no degrad option 
    394          IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN        ! eiv & bbl 
    395                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;        jf_eiw  = jfld + 3   ;   jfld = jf_eiw 
    396            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl  ;   slf_d(jf_eiw) = sn_eiw 
    397          ENDIF 
    398          IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN   ! no eiv & bbl 
     365      IF( ln_dynbbl ) THEN         ! eiv & bbl 
    399366                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl 
    400367           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    401          ENDIF 
    402          IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN   ! eiv & no bbl 
    403            jf_eiw = jfld + 1 ; jfld = jf_eiw ; slf_d(jf_eiw) = sn_eiw 
    404          ENDIF 
    405       ELSE 
    406               jf_ahu  = jfld + 1 ;        jf_ahv  = jfld + 2 ;        jf_ahw  = jfld + 3  ;  jfld = jf_ahw 
    407         slf_d(jf_ahu) = sn_ahu  ;   slf_d(jf_ahv) = sn_ahv  ;   slf_d(jf_ahw) = sn_ahw 
    408         IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN         ! eiv & bbl 
    409                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ; 
    410            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    411                  jf_eiu  = jfld + 3 ;        jf_eiv  = jfld + 4 ;    jf_eiw  = jfld + 5   ;  jfld = jf_eiw  
    412            slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;    slf_d(jf_eiw) = sn_eiw 
    413         ENDIF 
    414         IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN    ! no eiv & bbl 
    415                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl 
    416            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    417         ENDIF 
    418         IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN    ! eiv & no bbl 
    419                  jf_eiu  = jfld + 1 ;         jf_eiv  = jfld + 2 ;    jf_eiw  = jfld + 3   ; jfld = jf_eiw  
    420            slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;   slf_d(jf_eiw) = sn_eiw 
    421         ENDIF 
    422       ENDIF 
    423    
     368      ENDIF 
     369 
     370 
    424371      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure 
    425372      IF( ierr > 0 ) THEN 
    426373         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN 
    427374      ENDIF 
     375      !                                         ! fill sf with slf_i and control print 
     376      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 
    428377      ! Open file for each variable to get his number of dimension 
    429378      DO ifpr = 1, jfld 
    430          CALL iom_open( TRIM( cn_dir )//TRIM( slf_d(ifpr)%clname ), inum ) 
    431          idv   = iom_varid( inum , slf_d(ifpr)%clvar )  ! id of the variable sdjf%clvar 
    432          idimv = iom_file ( inum )%ndims(idv)             ! number of dimension for variable sdjf%clvar 
    433          IF( inum /= 0 )   CALL iom_close( inum )       ! close file if already open 
     379         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 
     380         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar 
     381         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar 
     382         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 
     383         ierr1=0 
    434384         IF( idimv == 3 ) THEN    ! 2D variable 
    435385                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 ) 
     
    443393         ENDIF 
    444394      END DO 
    445       !                                         ! fill sf with slf_i and control print 
    446       CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 
    447       ! 
    448       IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes  
     395      ! 
     396      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes  
    449397         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation 
    450398            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    & 
     
    505453               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) 
    506454               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) 
    507                zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     455               zet = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    508456               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet  
    509457            END DO 
    510458         END DO 
    511459      END DO 
     460      !                              !  update the horizontal divergence with the runoff inflow 
     461      IF( ln_dynrnf )  zhdiv(:,:,1) = zhdiv(:,:,1) - rnf(:,:) * r1_rau0 / fse3t(:,:,1) 
     462      ! 
    512463      CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv 
    513       ! 
    514464      ! computation of vertical velocity from the bottom 
    515465      pw(:,:,jpk) = 0._wp 
     
    534484      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes 
    535485      !!--------------------------------------------------------------------- 
    536 #if defined key_ldfslp && ! defined key_c1d 
    537       CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points 
    538       CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala 
    539       IF( ln_zps )   & ! Partial steps: before Horizontal DErivative 
    540         &    CALL zps_hde( kt, jpts, pts, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    541         &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    542         &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )  ! of t, s, rd at the last ocean level 
    543         ! only gtsu, gtsv, rhd, gru , grv are used  
    544  
    545  
    546          !                                                            ! of t, s, rd at the bottom ocean level 
    547       CALL zdf_mxl( kt )            ! mixed layer depth 
    548       CALL ldf_slp( kt, rhd, rn2 )  ! slopes 
    549       puslp (:,:,:) = uslp (:,:,:)  
    550       pvslp (:,:,:) = vslp (:,:,:)  
    551       pwslpi(:,:,:) = wslpi(:,:,:)  
    552       pwslpj(:,:,:) = wslpj(:,:,:)  
    553 #else 
    554       puslp (:,:,:) = 0.            ! to avoid warning when compiling 
    555       pvslp (:,:,:) = 0. 
    556       pwslpi(:,:,:) = 0. 
    557       pwslpj(:,:,:) = 0. 
    558 #endif 
     486      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                        
     487         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) ) 
     488         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points 
     489         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala 
     490 
     491         ! Partial steps: before Horizontal DErivative 
     492         IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
     493            &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     494            &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     495         IF( ln_zps .AND.        ln_isfcav)                            & 
     496            &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     497            &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     498            &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
     499 
     500         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
     501         CALL zdf_mxl( kt )            ! mixed layer depth 
     502         CALL ldf_slp( kt, rhd, rn2 )  ! slopes 
     503         puslp (:,:,:) = uslp (:,:,:)  
     504         pvslp (:,:,:) = vslp (:,:,:)  
     505         pwslpi(:,:,:) = wslpi(:,:,:)  
     506         pwslpj(:,:,:) = wslpj(:,:,:)  
     507     ELSE 
     508         puslp (:,:,:) = 0.            ! to avoid warning when compiling 
     509         pvslp (:,:,:) = 0. 
     510         pwslpi(:,:,:) = 0. 
     511         pwslpj(:,:,:) = 0. 
     512     ENDIF 
    559513      ! 
    560514   END SUBROUTINE dta_dyn_slp 
Note: See TracChangeset for help on using the changeset viewer.