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 14200 for NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DOM/domvvl.F90 – NEMO

Ignore:
Timestamp:
2020-12-17T15:36:44+01:00 (4 years ago)
Author:
mcastril
Message:

Merging r14117 through r14199 into dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DOM/domvvl.F90

    r14086 r14200  
    22   !!====================================================================== 
    33   !!                       ***  MODULE domvvl   *** 
    4    !! Ocean :  
     4   !! Ocean : 
    55   !!====================================================================== 
    66   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
     
    2626   USE timing          ! Timing 
    2727 
    28 #if defined key_agrif 
    29    USE agrif_oce       ! initial state interpolation 
    30    USE agrif_oce_interp  
    31 #endif 
    32  
    3328   IMPLICIT NONE 
    3429   PRIVATE 
     
    5550   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    5651 
    57 #if defined key_qco 
     52#if defined key_qco   ||   defined key_linssh 
    5853   !!---------------------------------------------------------------------- 
    59    !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate 
     54   !!   'key_qco'                        Quasi-Eulerian vertical coordinate 
     55   !!       OR         EMPTY MODULE 
     56   !!   'key_linssh'                        Fix in time vertical coordinate 
    6057   !!---------------------------------------------------------------------- 
    6158#else 
     
    6360   !!   Default key      Old management of time varying vertical coordinate 
    6461   !!---------------------------------------------------------------------- 
    65     
     62 
    6663   !!---------------------------------------------------------------------- 
    6764   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
     
    7875   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
    7976   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    80     
     77 
    8178   !! * Substitutions 
    8279#  include "do_loop_substitute.h90" 
     
    114111      !!---------------------------------------------------------------------- 
    115112      !!                ***  ROUTINE dom_vvl_init  *** 
    116       !!                    
     113      !! 
    117114      !! ** Purpose :  Initialization of all scale factors, depths 
    118115      !!               and water column heights 
     
    123120      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
    124121      !!              - Regrid: e3[u/v](:,:,:,Kmm) 
    125       !!                        e3[u/v](:,:,:,Kmm)        
    126       !!                        e3w(:,:,:,Kmm)            
     122      !!                        e3[u/v](:,:,:,Kmm) 
     123      !!                        e3w(:,:,:,Kmm) 
    127124      !!                        e3[u/v]w_b 
    128       !!                        e3[u/v]w_n       
     125      !!                        e3[u/v]w_n 
    129126      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 
    130127      !!              - h(t/u/v)_0 
     
    156153      !!---------------------------------------------------------------------- 
    157154      !!                ***  ROUTINE dom_vvl_init  *** 
    158       !!                    
    159       !! ** Purpose :  Interpolation of all scale factors,  
     155      !! 
     156      !! ** Purpose :  Interpolation of all scale factors, 
    160157      !!               depths and water column heights 
    161158      !! 
     
    164161      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
    165162      !!              - Regrid: e3(u/v)_n 
    166       !!                        e3(u/v)_b        
    167       !!                        e3w_n            
    168       !!                        e3(u/v)w_b       
    169       !!                        e3(u/v)w_n       
     163      !!                        e3(u/v)_b 
     164      !!                        e3w_n 
     165      !!                        e3(u/v)w_b 
     166      !!                        e3(u/v)w_n 
    170167      !!                        gdept_n, gdepw_n and gde3w_n 
    171168      !!              - h(t/u/v)_0 
     
    185182      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
    186183      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
    187       CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
     184      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V 
    188185      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    189186      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
    190       !                                ! Vertical interpolation of e3t,u,v  
     187      !                                ! Vertical interpolation of e3t,u,v 
    191188      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
    192189      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
     
    210207         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    211208         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    212          !                             ! 0.5 where jk = mikt      
     209         !                             ! 0.5 where jk = mikt 
    213210!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    214211         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    215212         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    216213         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    217             &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     214            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
    218215         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    219216         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    220217         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    221             &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     218            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
    222219      END_3D 
    223220      ! 
     
    278275            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    279276               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
    280                   ii0 = 103 + nn_hls - 1   ;   ii1 = 111 + nn_hls - 1       
     277                  ii0 = 103 + nn_hls - 1   ;   ii1 = 111 + nn_hls - 1 
    281278                  ij0 = 128 + nn_hls       ;   ij1 = 135 + nn_hls 
    282279                  frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp 
     
    290287 
    291288 
    292    SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )  
     289   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
    293290      !!---------------------------------------------------------------------- 
    294291      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
    295       !!                    
     292      !! 
    296293      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt, 
    297294      !!                 tranxt and dynspg routines 
    298295      !! 
    299296      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness. 
    300       !!               - z_tilde_case: after scale factor increment =  
     297      !!               - z_tilde_case: after scale factor increment = 
    301298      !!                                    high frequency part of horizontal divergence 
    302299      !!                                  + retsoring towards the background grid 
     
    306303      !! 
    307304      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case 
    308       !!               - tilde_e3t_a: after increment of vertical scale factor  
     305      !!               - tilde_e3t_a: after increment of vertical scale factor 
    309306      !!                              in z_tilde case 
    310307      !!               - e3(t/u/v)_a 
     
    410407            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    411408               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    412             vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
     409            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
    413410               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    414411            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     
    455452               WRITE(numout, *) 'at i, j, k=', ijk_max 
    456453               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
    457                WRITE(numout, *) 'at i, j, k=', ijk_min             
     454               WRITE(numout, *) 'at i, j, k=', ijk_min 
    458455               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high') 
    459456            ENDIF 
     
    571568      !!---------------------------------------------------------------------- 
    572569      !!                ***  ROUTINE dom_vvl_sf_update  *** 
    573       !!                    
    574       !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors  
     570      !! 
     571      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors 
    575572      !!               compute all depths and related variables for next time step 
    576573      !!               write outputs and restart file 
     
    582579      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step 
    583580      !!               - Recompute: 
    584       !!                    e3(u/v)_b        
    585       !!                    e3w(:,:,:,Kmm)            
    586       !!                    e3(u/v)w_b       
    587       !!                    e3(u/v)w_n       
     581      !!                    e3(u/v)_b 
     582      !!                    e3w(:,:,:,Kmm) 
     583      !!                    e3(u/v)w_b 
     584      !!                    e3(u/v)w_n 
    588585      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w 
    589586      !!                    h(u/v) and h(u/v)r 
     
    616613            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 
    617614         ELSE 
    618             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &  
     615            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
    619616            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 
    620617         ENDIF 
     
    628625      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 
    629626      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    630        
     627 
    631628      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  ) 
    632        
     629 
    633630      ! Vertical scale factor interpolations 
    634631      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
     
    649646         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    650647         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    651              &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     648             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
    652649         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    653650      END_3D 
     
    768765      !!--------------------------------------------------------------------- 
    769766      !!                   ***  ROUTINE dom_vvl_rst  *** 
    770       !!                      
     767      !! 
    771768      !! ** Purpose :   Read or write VVL file in restart file 
    772769      !! 
     
    800797         IF( ln_rstart ) THEN                   !==  Read the restart file  ==! 
    801798            ! 
    802 #if defined key_agrif 
    803             IF ( (.NOT.Agrif_root()).AND.(ln_init_chfrpar) ) THEN 
    804                ! skip reading restart if initialized from parent: 
    805                id3 = -1 ; id4 = -1 ; id5 = -1 
    806             ELSE 
    807 #endif 
    808799            CALL rst_read_open                                          !*  open the restart file if necessary 
    809800            !                                         ! --------- ! 
     
    815806            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. ) 
    816807            ! 
    817 #if defined key_agrif 
    818             ENDIF 
    819 #endif 
    820808            !                                                           !*  scale factors 
    821809            IF(lwp) WRITE(numout,*)    '          Kmm scale factor read in the restart file' 
    822810            CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) ) 
    823             WHERE ( tmask(:,:,:) == 0.0_wp )  
     811            WHERE ( tmask(:,:,:) == 0.0_wp ) 
    824812               e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
    825813            END WHERE 
    826814            IF( l_1st_euler ) THEN                       ! euler 
    827815               IF(lwp) WRITE(numout,*) '          Euler first time step : e3t(Kbb) = e3t(Kmm)' 
    828                e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     816               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    829817            ELSE                                         ! leap frog 
    830818               IF(lwp) WRITE(numout,*) '          Kbb scale factor read in the restart file' 
    831819               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) ) 
    832                WHERE ( tmask(:,:,:) == 0.0_wp )  
     820               WHERE ( tmask(:,:,:) == 0.0_wp ) 
    833821                  e3t(:,:,:,Kbb) = e3t_0(:,:,:) 
    834822               END WHERE 
     
    854842                     CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 
    855843                  ENDIF 
    856                ELSE  
    857                  tilde_e3t_b(:,:,:) = 0.0_wp 
    858                  tilde_e3t_n(:,:,:) = 0.0_wp 
     844               ELSE 
     845                  tilde_e3t_b(:,:,:) = 0.0_wp 
     846                  tilde_e3t_n(:,:,:) = 0.0_wp 
    859847               ENDIF 
    860848               !                                      ! ------------ ! 
     
    864852                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:) ) 
    865853                  ELSE                ! array is missing 
    866                      hdiv_lf(:,:,:) = 0.0_wp  
     854                     hdiv_lf(:,:,:) = 0.0_wp 
    867855                  ENDIF 
    868856               ENDIF 
     
    885873      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN       !  Create restart file  ! 
    886874         !                                       !=======================! 
    887 #if defined key_agrif 
    888          IF ( .NOT.Agrif_root().AND.(ln_init_chfrpar) ) THEN 
    889             ! Interpolate initial ssh from parent: 
    890             CALL Agrif_istate_ssh( Kbb, Kmm ) 
    891             ! 
    892             DO jk = 1, jpk 
    893                e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) & 
    894                  &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    895                  &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    896             END DO 
    897             e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    898          ENDIF 
    899 #endif 
    900875         ! 
    901876         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 
     
    911886            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:)) 
    912887         END IF 
    913          !                                           ! -------------!     
     888         !                                           ! -------------! 
    914889         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    915890            !                                        ! ------------ ! 
     
    925900      !!--------------------------------------------------------------------- 
    926901      !!                  ***  ROUTINE dom_vvl_ctl  *** 
    927       !!                 
     902      !! 
    928903      !! ** Purpose :   Control the consistency between namelist options 
    929904      !!                for vertical coordinate 
     
    934909         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
    935910         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
    936       !!----------------------------------------------------------------------  
     911      !!---------------------------------------------------------------------- 
    937912      ! 
    938913      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 
Note: See TracChangeset for help on using the changeset viewer.