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 13151 – NEMO

Changeset 13151


Ignore:
Timestamp:
2020-06-24T14:38:26+02:00 (4 years ago)
Author:
gm
Message:

result from merge with qco r12983

Location:
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src
Files:
8 added
163 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/ABL/ablrst.F90

    r11945 r13151  
    7474            ENDIF 
    7575            ! 
    76             CALL iom_open( TRIM(clpath)//TRIM(clname), numraw, ldwrt = .TRUE., kdlev = jpka ) 
     76            CALL iom_open( TRIM(clpath)//TRIM(clname), numraw, ldwrt = .TRUE., kdlev = jpka, cdcomp = 'ABL' ) 
    7777            lrst_abl = .TRUE. 
    7878         ENDIF 
    … …  
    146146      ENDIF 
    147147 
    148       CALL iom_open ( TRIM(cn_ablrst_indir)//'/'//cn_ablrst_in, numrar, kdlev = jpka ) 
     148      CALL iom_open ( TRIM(cn_ablrst_indir)//'/'//cn_ablrst_in, numrar ) 
    149149 
    150150      ! Time info 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/ABL/sbcabl.F90

    r12489 r13151  
    7575      !!--------------------------------------------------------------------- 
    7676 
    77       REWIND( numnam_ref )              ! Namelist namsbc_abl in reference namelist : ABL parameters 
     77      ! Namelist namsbc_abl in reference namelist : ABL parameters 
    7878      READ  ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 ) 
    7979901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' ) 
    80       ! 
    81       REWIND( numnam_cfg )              ! Namelist namsbc_abl in configuration namelist : ABL parameters 
     80      ! Namelist namsbc_abl in configuration namelist : ABL parameters 
    8281      READ  ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 
    8382902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/ICE/icectl.F90

    r12489 r13151  
    331331      IF(lwp) WRITE(numout,*)                 
    332332 
    333       CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
     333      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 
    334334       
    335335      CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 )    ! ice mass spurious lost/gain 
    … …  
    725725       
    726726      CALL prt_ctl_info(' ') 
    727       CALL prt_ctl_info(' - Heat / FW fluxes : ') 
    728       CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ') 
    729       CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ') 
    730       CALL prt_ctl(tab2d_1=qsr    , clinfo1= ' qsr   : ', tab2d_2=qns       , clinfo2= ' qns       : ') 
    731       CALL prt_ctl(tab2d_1=emp    , clinfo1= ' emp   : ', tab2d_2=sfx       , clinfo2= ' sfx       : ') 
    732        
    733       CALL prt_ctl_info(' ') 
    734727      CALL prt_ctl_info(' - Stresses : ') 
    735728      CALL prt_ctl_info('   ~~~~~~~~~~ ') 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/ICE/icedyn_rhg_evp.F90

    r12489 r13151  
    4949   !! * Substitutions 
    5050#  include "do_loop_substitute.h90" 
     51#  include "domzgr_substitute.h90" 
    5152   !!---------------------------------------------------------------------- 
    5253   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/ICE/iceistate.F90

    r12489 r13151  
    1818   USE oce            ! dynamics and tracers variables 
    1919   USE dom_oce        ! ocean domain 
    20    USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd  
     20   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 
    2121   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b 
    2222   USE eosbn2         ! equation of state 
    … …  
    6060   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m) 
    6161   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read) 
    62    !    
     62 
    6363   !! * Substitutions 
    6464#  include "do_loop_substitute.h90" 
    … …  
    7777      !! 
    7878      !! ** Method  :   This routine will put some ice where ocean 
    79       !!                is at the freezing point, then fill in ice  
    80       !!                state variables using prescribed initial  
    81       !!                values in the namelist             
     79      !!                is at the freezing point, then fill in ice 
     80      !!                state variables using prescribed initial 
     81      !!                values in the namelist 
    8282      !! 
    8383      !! ** Steps   :   1) Set initial surface and basal temperatures 
    … …  
    9191      !!              where there is no ice 
    9292      !!-------------------------------------------------------------------- 
    93       INTEGER, INTENT(in) :: kt            ! time step  
     93      INTEGER, INTENT(in) :: kt            ! time step 
    9494      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 
    9595      ! 
    … …  
    102102      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    103103      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file 
    104       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays 
     104      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !locak arrays 
    105105      !! 
    106106      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 
    … …  
    117117      ! basal temperature (considered at freezing point)   [Kelvin] 
    118118      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
    119       t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
     119      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
    120120      ! 
    121121      ! surface temperature and conductivity 
    … …  
    142142      e_i (:,:,:,:) = 0._wp 
    143143      e_s (:,:,:,:) = 0._wp 
    144        
     144 
    145145      ! general fields 
    146146      a_i (:,:,:) = 0._wp 
    … …  
    213213            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 
    214214               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 
    215                &                              * si(jp_ati)%fnow(:,:,1)  
     215               &                              * si(jp_ati)%fnow(:,:,1) 
    216216            ! 
    217217            ! pond depth 
    … …  
    227227            ! 
    228228            ! change the switch for the following 
    229             WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1)  
     229            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
    230230            ELSEWHERE                         ;   zswitch(:,:) = 0._wp 
    231231            END WHERE 
    … …  
    234234            !                          !---------------! 
    235235            ! no ice if (sst - Tfreez) >= thresold 
    236             WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp  
     236            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
    237237            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1) 
    238238            END WHERE 
    … …  
    247247               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 
    248248               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
    249                zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
     249               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
    250250               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
    251251            ELSEWHERE 
    … …  
    268268            zhpnd_ini(:,:) = 0._wp 
    269269         ENDIF 
    270           
     270 
    271271         !-------------! 
    272272         ! fill fields ! 
    … …  
    295295         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 
    296296            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 
    297           
     297 
    298298         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
    299299         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   & 
    … …  
    341341         DO jl = 1, jpl 
    342342            DO_3D_11_11( 1, nlay_i ) 
    343                t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)  
     343               t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
    344344               ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 
    345345               e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 
    … …  
    357357         END WHERE 
    358358         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    359            
     359 
    360360         ! specific temperatures for coupled runs 
    361361         tn_ice(:,:,:) = t_su(:,:,:) 
    … …  
    377377         ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0 
    378378         ! 
    379          IF( .NOT.ln_linssh ) THEN 
    380             ! 
    381             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:) 
    382             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE 
    383             ! 
    384             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                 
    385                e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:) 
    386                e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
    387                e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm) 
    388             END DO 
    389             ! 
    390             ! Reconstruction of all vertical scale factors at now and before time-steps 
    391             ! ========================================================================= 
    392             ! Horizontal scale factor interpolations 
    393             ! -------------------------------------- 
    394             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 
    395             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 
    396             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
    397             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    398             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 
    399             ! Vertical scale factor interpolations 
    400             ! ------------------------------------ 
    401             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  ) 
    402             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
    403             CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
    404             CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
    405             CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    406             ! t- and w- points depth 
    407             ! ---------------------- 
    408             !!gm not sure of that.... 
    409             gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
    410             gdepw(:,:,1,Kmm) = 0.0_wp 
    411             gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    412             DO jk = 2, jpk 
    413                gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm) 
    414                gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 
    415                gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm) 
    416             END DO 
    417          ENDIF 
     379         IF( .NOT.ln_linssh )   CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column 
     380! !!st 
     381!          IF( .NOT.ln_linssh ) THEN 
     382!             ! 
     383!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:) 
     384!             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE 
     385!             ! 
     386!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
     387!                e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:) 
     388!                e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
     389!                e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm) 
     390!             END DO 
     391!             ! 
     392!             ! Reconstruction of all vertical scale factors at now and before time-steps 
     393!             ! ========================================================================= 
     394!             ! Horizontal scale factor interpolations 
     395!             ! -------------------------------------- 
     396!             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 
     397!             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 
     398!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     399!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     400!             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 
     401!             ! Vertical scale factor interpolations 
     402!             ! ------------------------------------ 
     403!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  ) 
     404!             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     405!             CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
     406!             CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     407!             CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
     408!             ! t- and w- points depth 
     409!             ! ---------------------- 
     410!             !!gm not sure of that.... 
     411!             gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     412!             gdepw(:,:,1,Kmm) = 0.0_wp 
     413!             gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
     414!             DO jk = 2, jpk 
     415!                gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm) 
     416!                gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 
     417!                gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm) 
     418!             END DO 
     419!          ENDIF 
    418420      ENDIF 
    419        
     421 
    420422      !------------------------------------ 
    421423      ! 4) store fields at before time-step 
    … …  
    432434      v_ice_b(:,:)     = v_ice(:,:) 
    433435      ! total concentration is needed for Lupkes parameterizations 
    434       at_i_b (:,:)     = at_i (:,:)  
     436      at_i_b (:,:)     = at_i (:,:) 
    435437 
    436438!!clem: output of initial state should be written here but it is impossible because 
    437439!!      the ocean and ice are in the same file 
    438 !!      CALL dia_wri_state( 'output.init' ) 
     440!!      CALL dia_wri_state( Kmm, 'output.init' ) 
    439441      ! 
    440442   END SUBROUTINE ice_istate 
    … …  
    444446      !!------------------------------------------------------------------- 
    445447      !!                   ***  ROUTINE ice_istate_init  *** 
    446       !!         
    447       !! ** Purpose :   Definition of initial state of the ice  
    448       !! 
    449       !! ** Method  :   Read the namini namelist and check the parameter  
     448      !! 
     449      !! ** Purpose :   Definition of initial state of the ice 
     450      !! 
     451      !! ** Method  :   Read the namini namelist and check the parameter 
    450452      !!              values called at the first timestep (nit000) 
    451453      !! 
    … …  
    453455      !! 
    454456      !!----------------------------------------------------------------------------- 
    455       INTEGER ::   ios   ! Local integer output status for namelist read 
    456       INTEGER ::   ifpr, ierror 
     457      INTEGER ::   ios, ifpr, ierror   ! Local integers 
     458 
    457459      ! 
    458460      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files 
    … …  
    488490         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst 
    489491         IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 
    490             WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s  
     492            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
    491493            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s 
    492494            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/ICE/icerst.F90

    r12377 r13151  
    8080            ENDIF 
    8181            ! 
    82             CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kdlev = jpl ) 
     82            CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 
    8383            lrst_ice = .TRUE. 
    8484         ENDIF 
    … …  
    185185      ENDIF 
    186186 
    187       CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kdlev = jpl ) 
     187      CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir ) 
    188188 
    189189      ! test if v_i exists  

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/ASM/asminc.F90

    r12489 r13151  
    99   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR 
    1010   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2 
    11    !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init  
     11   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init 
    1212   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization 
    1313   !!                 ! 2014-09  (D. Lea)  Local calc_date removed use routine from OBS 
    … …  
    3131   USE zpshde          ! Partial step : Horizontal Derivative 
    3232   USE asmpar          ! Parameters for the assmilation interface 
    33    USE asmbkg          !  
     33   USE asmbkg          ! 
    3434   USE c1d             ! 1D initialization 
    3535   USE sbc_oce         ! Surface boundary condition variables. 
    … …  
    4545   IMPLICIT NONE 
    4646   PRIVATE 
    47     
     47 
    4848   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights 
    4949   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments 
    … …  
    7272   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components 
    7373   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S 
    74    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components  
     74   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components 
    7575   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step 
    7676#if defined key_asminc 
    … …  
    8080   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term 
    8181   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization 
    82    INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval  
     82   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
    8383   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval 
    84    !  
     84   ! 
    8585   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting 
    86    !                                 !: = 1   Linear hat-like, centred in middle of IAU interval  
     86   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval 
    8787   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix) 
    8888 
    … …  
    9595   !! * Substitutions 
    9696#  include "do_loop_substitute.h90" 
     97#  include "domzgr_substitute.h90" 
    9798   !!---------------------------------------------------------------------- 
    9899   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    105106      !!---------------------------------------------------------------------- 
    106107      !!                    ***  ROUTINE asm_inc_init  *** 
    107       !!           
     108      !! 
    108109      !! ** Purpose : Initialize the assimilation increment and IAU weights. 
    109110      !! 
    110111      !! ** Method  : Initialize the assimilation increment and IAU weights. 
    111112      !! 
    112       !! ** Action  :  
     113      !! ** Action  : 
    113114      !!---------------------------------------------------------------------- 
    114115      INTEGER, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices 
    … …  
    262263         ! 
    263264         !                                !--------------------------------------------------------- 
    264          IF( niaufn == 0 ) THEN           ! Constant IAU forcing  
     265         IF( niaufn == 0 ) THEN           ! Constant IAU forcing 
    265266            !                             !--------------------------------------------------------- 
    266267            DO jt = 1, iiauper 
    … …  
    268269            END DO 
    269270            !                             !--------------------------------------------------------- 
    270          ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval  
     271         ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval 
    271272            !                             !--------------------------------------------------------- 
    272273            ! Compute the normalization factor 
    273274            znorm = 0._wp 
    274275            IF( MOD( iiauper, 2 ) == 0 ) THEN   ! Even number of time steps in IAU interval 
    275                imid = iiauper / 2  
     276               imid = iiauper / 2 
    276277               DO jt = 1, imid 
    277278                  znorm = znorm + REAL( jt ) 
    … …  
    279280               znorm = 2.0 * znorm 
    280281            ELSE                                ! Odd number of time steps in IAU interval 
    281                imid = ( iiauper + 1 ) / 2         
     282               imid = ( iiauper + 1 ) / 2 
    282283               DO jt = 1, imid - 1 
    283284                  znorm = znorm + REAL( jt ) 
    … …  
    306307             DO jt = 1, icycper 
    307308                ztotwgt = ztotwgt + wgtiau(jt) 
    308                 WRITE(numout,*) '         ', jt, '       ', wgtiau(jt)  
    309              END DO    
     309                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
     310             END DO 
    310311             WRITE(numout,*) '         ===================================' 
    311312             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt 
    312313             WRITE(numout,*) '         ===================================' 
    313314          ENDIF 
    314           
     315 
    315316      ENDIF 
    316317 
    … …  
    337338         CALL iom_open( c_asminc, inum ) 
    338339         ! 
    339          CALL iom_get( inum, 'time'       , zdate_inc   )  
     340         CALL iom_get( inum, 'time'       , zdate_inc   ) 
    340341         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb ) 
    341342         CALL iom_get( inum, 'z_inc_datef', z_inc_datef ) 
    … …  
    344345         ! 
    345346         IF(lwp) THEN 
    346             WRITE(numout,*)  
     347            WRITE(numout,*) 
    347348            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef 
    348349            WRITE(numout,*) '~~~~~~~~~~~~' 
    … …  
    358359            &                ' not agree with Direct Initialization time' ) 
    359360 
    360          IF ( ln_trainc ) THEN    
     361         IF ( ln_trainc ) THEN 
    361362            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
    362363            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
    … …  
    370371         ENDIF 
    371372 
    372          IF ( ln_dyninc ) THEN    
    373             CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )               
    374             CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )               
     373         IF ( ln_dyninc ) THEN 
     374            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 ) 
     375            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 ) 
    375376            ! Apply the masks 
    376377            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:) 
    … …  
    381382            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0 
    382383         ENDIF 
    383          
     384 
    384385         IF ( ln_sshinc ) THEN 
    385386            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 ) 
    … …  
    407408      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter 
    408409         !                                         !-------------------------------------- 
    409          ALLOCATE( zhdiv(jpi,jpj) )  
     410         ALLOCATE( zhdiv(jpi,jpj) ) 
    410411         ! 
    411412         DO jt = 1, nn_divdmp 
    … …  
    417418                     &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    & 
    418419                     &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    & 
    419                      &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) / e3t(ji,jj,jk,Kmm) 
     420                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) & 
     421                     &            / e3t(ji,jj,jk,Kmm) 
    420422               END_2D 
    421423               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change) 
    … …  
    425427                     &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
    426428                  v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         & 
    427                      &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)  
     429                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
    428430               END_2D 
    429431            END DO 
    … …  
    431433         END DO 
    432434         ! 
    433          DEALLOCATE( zhdiv )  
     435         DEALLOCATE( zhdiv ) 
    434436         ! 
    435437      ENDIF 
    … …  
    452454         CALL iom_open( c_asmdin, inum ) 
    453455         ! 
    454          CALL iom_get( inum, 'rdastp', zdate_bkg )  
     456         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
    455457         ! 
    456458         IF(lwp) THEN 
    457             WRITE(numout,*)  
     459            WRITE(numout,*) 
    458460            WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg 
    459461            WRITE(numout,*) 
    … …  
    464466            &                ' not agree with Direct Initialization time' ) 
    465467         ! 
    466          IF ( ln_trainc ) THEN    
     468         IF ( ln_trainc ) THEN 
    467469            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) 
    468470            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg ) 
    … …  
    471473         ENDIF 
    472474         ! 
    473          IF ( ln_dyninc ) THEN    
     475         IF ( ln_dyninc ) THEN 
    474476            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) 
    475477            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg ) 
    … …  
    499501      ! 
    500502   END SUBROUTINE asm_inc_init 
    501     
    502     
     503 
     504 
    503505   SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 
    504506      !!---------------------------------------------------------------------- 
    505507      !!                    ***  ROUTINE tra_asm_inc  *** 
    506       !!           
     508      !! 
    507509      !! ** Purpose : Apply the tracer (T and S) assimilation increments 
    508510      !! 
    509511      !! ** Method  : Direct initialization or Incremental Analysis Updating 
    510512      !! 
    511       !! ** Action  :  
     513      !! ** Action  : 
    512514      !!---------------------------------------------------------------------- 
    513515      INTEGER                                  , INTENT(in   ) :: kt             ! Current time step 
    … …  
    521523      !!---------------------------------------------------------------------- 
    522524      ! 
    523       ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)  
    524       ! used to prevent the applied increments taking the temperature below the local freezing point  
     525      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 
     526      ! used to prevent the applied increments taking the temperature below the local freezing point 
    525527      DO jk = 1, jpkm1 
    526528        CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 
    … …  
    537539            ! 
    538540            IF(lwp) THEN 
    539                WRITE(numout,*)  
     541               WRITE(numout,*) 
    540542               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 
    541543               WRITE(numout,*) '~~~~~~~~~~~~' 
    … …  
    547549                  ! Do not apply negative increments if the temperature will fall below freezing 
    548550                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 
    549                      &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )  
    550                      pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt   
     551                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
     552                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 
    551553                  END WHERE 
    552554               ELSE 
    553                   pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt   
     555                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt 
    554556               ENDIF 
    555557               IF (ln_salfix) THEN 
    … …  
    557559                  ! minimum value salfixmin 
    558560                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 
    559                      &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )  
     561                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
    560562                     pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 
    561563                  END WHERE 
    … …  
    574576      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization 
    575577         !                             !-------------------------------------- 
    576          !             
     578         ! 
    577579         IF ( kt == nitdin_r ) THEN 
    578580            ! 
    … …  
    582584            IF (ln_temnofreeze) THEN 
    583585               ! Do not apply negative increments if the temperature will fall below freezing 
    584                WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
    585                   pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     586               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
     587                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 
    586588               END WHERE 
    587589            ELSE 
    588                pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     590               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:) 
    589591            ENDIF 
    590592            IF (ln_salfix) THEN 
    591593               ! Do not apply negative increments if the salinity will fall below a specified 
    592594               ! minimum value salfixmin 
    593                WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin )  
    594                   pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     595               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin ) 
     596                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 
    595597               END WHERE 
    596598            ELSE 
    597                pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     599               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:) 
    598600            ENDIF 
    599601 
    … …  
    617619            DEALLOCATE( s_bkg    ) 
    618620         ENDIF 
    619          !   
     621         ! 
    620622      ENDIF 
    621623      ! Perhaps the following call should be in step 
    … …  
    628630      !!---------------------------------------------------------------------- 
    629631      !!                    ***  ROUTINE dyn_asm_inc  *** 
    630       !!           
     632      !! 
    631633      !! ** Purpose : Apply the dynamics (u and v) assimilation increments. 
    632634      !! 
    633635      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    634636      !! 
    635       !! ** Action  :  
     637      !! ** Action  : 
    636638      !!---------------------------------------------------------------------- 
    637639      INTEGER                             , INTENT( in )  ::  kt             ! ocean time-step index 
    … …  
    654656            ! 
    655657            IF(lwp) THEN 
    656                WRITE(numout,*)  
     658               WRITE(numout,*) 
    657659               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 
    658660               WRITE(numout,*) '~~~~~~~~~~~~' 
    … …  
    674676      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization 
    675677         !                          !----------------------------------------- 
    676          !          
     678         ! 
    677679         IF ( kt == nitdin_r ) THEN 
    678680            ! 
    … …  
    681683            ! Initialize the now fields with the background + increment 
    682684            puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 
    683             pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:)   
     685            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
    684686            ! 
    685687            puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields 
    … …  
    700702      !!---------------------------------------------------------------------- 
    701703      !!                    ***  ROUTINE ssh_asm_inc  *** 
    702       !!           
     704      !! 
    703705      !! ** Purpose : Apply the sea surface height assimilation increment. 
    704706      !! 
    705707      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    706708      !! 
    707       !! ** Action  :  
     709      !! ** Action  : 
    708710      !!---------------------------------------------------------------------- 
    709711      INTEGER, INTENT(IN) :: kt         ! Current time step 
    … …  
    725727            ! 
    726728            IF(lwp) THEN 
    727                WRITE(numout,*)  
     729               WRITE(numout,*) 
    728730               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', & 
    729731                  &  kt,' with IAU weight = ', wgtiau(it) 
    … …  
    758760            ! 
    759761            ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields 
     762#if ! defined key_qco 
    760763            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     764#endif 
    761765!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 
    762766            ! 
    … …  
    775779      !!                  ***  ROUTINE ssh_asm_div  *** 
    776780      !! 
    777       !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence           
     781      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence 
    778782      !!                across all the water column 
    779783      !! 
    … …  
    791795      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array 
    792796      !!---------------------------------------------------------------------- 
    793       !  
     797      ! 
    794798#if defined key_asminc 
    795799      CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments) 
    796800      ! 
    797       IF( ln_linssh ) THEN  
     801      IF( ln_linssh ) THEN 
    798802         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 
    799       ELSE  
     803      ELSE 
    800804         ALLOCATE( ztim(jpi,jpj) ) 
    801805         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 
    802          DO jk = 1, jpkm1                                  
    803             phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)  
     806         DO jk = 1, jpkm1 
     807            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
    804808         END DO 
    805809         ! 
    … …  
    814818      !!---------------------------------------------------------------------- 
    815819      !!                    ***  ROUTINE seaice_asm_inc  *** 
    816       !!           
     820      !! 
    817821      !! ** Purpose : Apply the sea ice assimilation increment. 
    818822      !! 
    819823      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    820824      !! 
    821       !! ** Action  :  
     825      !! ** Action  : 
    822826      !! 
    823827      !!---------------------------------------------------------------------- 
    … …  
    840844            ! 
    841845            it = kt - nit000 + 1 
    842             zincwgt = wgtiau(it)      ! IAU weight for the current time step  
     846            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
    843847            ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 
    844848            ! 
    845849            IF(lwp) THEN 
    846                WRITE(numout,*)  
     850               WRITE(numout,*) 
    847851               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 
    848852               WRITE(numout,*) '~~~~~~~~~~~~' 
    … …  
    862866            ! 
    863867            ! Nudge sea ice depth to bring it up to a required minimum depth 
    864             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin )  
    865                zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt     
     868            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
     869               zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt 
    866870            ELSEWHERE 
    867871               zhicifinc(:,:) = 0.0_wp 
    … …  
    896900         IF ( kt == nitdin_r ) THEN 
    897901            ! 
     902<<<<<<< .working 
    898903            l_1st_euler = 0              ! Force Euler forward step 
     904======= 
     905            l_1st_euler = .TRUE.              ! Force Euler forward step 
     906>>>>>>> .merge-right.r13092 
    899907            ! 
    900908            ! Sea-ice : SI3 case 
    … …  
    903911            zofrld (:,:) = 1._wp - at_i(:,:) 
    904912            zohicif(:,:) = hm_i(:,:) 
    905             !  
     913            ! 
    906914            ! Initialize the now fields the background + increment 
    907915            at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 
    908             at_i_b(:,:) = at_i(:,:)  
     916            at_i_b(:,:) = at_i(:,:) 
    909917            fr_i(:,:) = at_i(:,:)        ! adjust ice fraction 
    910918            ! 
    … …  
    912920            ! 
    913921            ! Nudge sea ice depth to bring it up to a required minimum depth 
    914             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin )  
     922            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin ) 
    915923               zhicifinc(:,:) = zhicifmin - hm_i(:,:) 
    916924            ELSEWHERE 
    … …  
    942950!#if defined defined key_si3 || defined key_cice 
    943951! 
    944 !            IF (ln_seaicebal ) THEN        
     952!            IF (ln_seaicebal ) THEN 
    945953!             !! balancing salinity increments 
    946954!             !! simple case from limflx.F90 (doesn't include a mass flux) 
    … …  
    954962! 
    955963!             DO jj = 1, jpj 
    956 !               DO ji = 1, jpi  
     964!               DO ji = 1, jpi 
    957965!           ! calculate change in ice and snow mass per unit area 
    958966!           ! positive values imply adding salt to the ocean (results from ice formation) 
    … …  
    965973! 
    966974!           ! prevent small mld 
    967 !           ! less than 10m can cause salinity instability  
     975!           ! less than 10m can cause salinity instability 
    968976!                 IF (mld < 10) mld=10 
    969977! 
    970 !           ! set to bottom of a level  
     978!           ! set to bottom of a level 
    971979!                 DO jk = jpk-1, 2, -1 
    972 !                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN  
    973 !                     mld=gdepw(ji,jj,jk+1) 
     980!                   IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN 
     981!                     mld=gdepw(ji,jj,jk+1,Kmm) 
    974982!                     jkmax=jk 
    975983!                   ENDIF 
    … …  
    977985! 
    978986!            ! avoid applying salinity balancing in shallow water or on land 
    979 !            !  
     987!            ! 
    980988! 
    981989!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 
    … …  
    988996! 
    989997!           ! put increments in for levels in the mixed layer 
    990 !           ! but prevent salinity below a threshold value  
    991 ! 
    992 !                   DO jk = 1, jkmax               
    993 ! 
    994 !                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN  
     998!           ! but prevent salinity below a threshold value 
     999! 
     1000!                   DO jk = 1, jkmax 
     1001! 
     1002!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 
    9951003!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 
    9961004!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 
    … …  
    10031011!      ! 
    10041012!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean 
    1005 !      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt   
    1006 !      !!                
    1007 !      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)  
     1013!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
     1014!      !! 
     1015!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d) 
    10081016!      !!                                                     ! E-P (kg m-2 s-2) 
    10091017!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2) 
    … …  
    10181026      ! 
    10191027   END SUBROUTINE seaice_asm_inc 
    1020     
     1028 
    10211029   !!====================================================================== 
    10221030END MODULE asminc 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/BDY/bdydta.F90

    r12396 r13151  
    7070   !! * Substitutions 
    7171#  include "do_loop_substitute.h90" 
     72#  include "domzgr_substitute.h90" 
    7273   !!---------------------------------------------------------------------- 
    7374   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    9293      INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers 
    9394      INTEGER,   DIMENSION(jpbgrd)     ::   ilen1  
    94       INTEGER,   DIMENSION(:), POINTER ::   nblen, nblenrim  ! short cuts 
    9595      TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut 
    9696      TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias 
    … …  
    108108         DO jbdy = 1, nb_bdy 
    109109            ! 
    110             nblen    => idx_bdy(jbdy)%nblen 
    111             nblenrim => idx_bdy(jbdy)%nblenrim 
    112             ! 
    113110            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN  
    114                ilen1(:) = nblen(:) 
    115111               IF( dta_bdy(jbdy)%lneed_ssh ) THEN  
    116112                  igrd = 1 
    117                   DO ib = 1, ilen1(igrd) 
     113                  DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)   ! ssh is allocated and used only on the rim 
    118114                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    119115                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    … …  
    121117                  END DO 
    122118               ENDIF 
    123                IF( dta_bdy(jbdy)%lneed_dyn2d) THEN  
     119               IF( dta_bdy(jbdy)%lneed_dyn2d .AND. ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer 
    124120                  igrd = 2 
    125                   DO ib = 1, ilen1(igrd) 
     121                  DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)   ! u2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init 
    126122                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    127123                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    … …  
    129125                  END DO 
    130126                  igrd = 3 
    131                   DO ib = 1, ilen1(igrd) 
     127                  DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)   ! v2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init 
    132128                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    133129                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    … …  
    138134            ! 
    139135            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN  
    140                ilen1(:) = nblen(:) 
    141136               IF( dta_bdy(jbdy)%lneed_dyn3d ) THEN  
    142137                  igrd = 2  
    143                   DO ib = 1, ilen1(igrd) 
     138                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    144139                     DO ik = 1, jpkm1 
    145140                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    … …  
    149144                  END DO 
    150145                  igrd = 3  
    151                   DO ib = 1, ilen1(igrd) 
     146                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    152147                     DO ik = 1, jpkm1 
    153148                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    … …  
    160155 
    161156            IF( nn_tra_dta(jbdy) == 0 ) THEN  
    162                ilen1(:) = nblen(:) 
    163157               IF( dta_bdy(jbdy)%lneed_tra ) THEN 
    164158                  igrd = 1  
    165                   DO ib = 1, ilen1(igrd) 
     159                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    166160                     DO ik = 1, jpkm1 
    167161                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    … …  
    176170#if defined key_si3 
    177171            IF( nn_ice_dta(jbdy) == 0 ) THEN    ! set ice to initial values 
    178                ilen1(:) = nblen(:) 
    179172               IF( dta_bdy(jbdy)%lneed_ice ) THEN 
    180173                  igrd = 1    
    181174                  DO jl = 1, jpl 
    182                      DO ib = 1, ilen1(igrd) 
     175                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    183176                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    184177                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    … …  
    236229         ! tidal harmonic forcing ONLY: initialise arrays 
    237230         IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d  
    238             IF( dta_alias%lneed_ssh   ) dta_alias%ssh(:) = 0._wp 
    239             IF( dta_alias%lneed_dyn2d ) dta_alias%u2d(:) = 0._wp 
    240             IF( dta_alias%lneed_dyn2d ) dta_alias%v2d(:) = 0._wp 
     231            IF( dta_alias%lneed_ssh   .AND. ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp 
     232            IF( dta_alias%lneed_dyn2d .AND. ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp 
     233            IF( dta_alias%lneed_dyn2d .AND. ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp 
    241234         ENDIF 
    242235 
    … …  
    245238            ! 
    246239            igrd = 2                       ! zonal velocity 
    247             dta_alias%u2d(:) = 0._wp       ! compute barotrope zonal velocity and put it in u2d 
    248240            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    249241               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    250242               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     243               dta_alias%u2d(ib) = 0._wp   ! compute barotrope zonal velocity and put it in u2d 
    251244               DO ik = 1, jpkm1 
    252                   dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) 
     245                  dta_alias%u2d(ib) = dta_alias%u2d(ib)   & 
     246                     &              + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) 
    253247               END DO 
    254248               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm) 
    … …  
    258252            END DO 
    259253            igrd = 3                       ! meridional velocity 
    260             dta_alias%v2d(:) = 0._wp       ! compute barotrope meridional velocity and put it in v2d 
    261254            DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    262255               ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    263256               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     257               dta_alias%v2d(ib) = 0._wp   ! compute barotrope meridional velocity and put it in v2d 
    264258               DO ik = 1, jpkm1 
    265                   dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) 
     259                  dta_alias%v2d(ib) = dta_alias%v2d(ib)   & 
     260                     &              + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) 
    266261               END DO 
    267262               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm) 
    … …  
    283278 
    284279#if defined key_si3 
    285          IF( dta_alias%lneed_ice ) THEN 
     280         IF( dta_alias%lneed_ice .AND. idx_bdy(jbdy)%nblen(1) > 0 ) THEN 
    286281            ! fill temperature and salinity arrays 
    287282            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem (jbdy) 
    … …  
    338333            DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop 
    339334               IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 
    340                   nblen => idx_bdy(jbdy)%nblen 
    341                   nblenrim => idx_bdy(jbdy)%nblenrim 
    342                   IF( cn_dyn2d(jbdy) == 'frs' ) THEN   ;   ilen1(:)=nblen(:) 
    343                   ELSE                                 ;   ilen1(:)=nblenrim(:) 
     335                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN   ;   ilen1(:)=idx_bdy(jbdy)%nblen(:) 
     336                  ELSE                                 ;   ilen1(:)=idx_bdy(jbdy)%nblenrim(:) 
    344337                  ENDIF 
    345338                  IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/BDY/bdydyn.F90

    r12377 r13151  
    2929 
    3030   PUBLIC   bdy_dyn    ! routine called in dyn_nxt 
    31  
     31    
     32   !! * Substitutions 
     33#  include "domzgr_substitute.h90" 
    3234   !!---------------------------------------------------------------------- 
    3335   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/C1D/step_c1d.F90

    r12377 r13151  
    8383      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )  ! after vertical scale factors  
    8484 
    85       IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity  
     85      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, Naa, ww )  ! now cross-level velocity  
    8686      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    8787      ! diagnostics and outputs        

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/CRS/crsfld.F90

    r12377 r13151  
    3333   !! * Substitutions 
    3434#  include "do_loop_substitute.h90" 
     35#  include "domzgr_substitute.h90" 
    3536   !!---------------------------------------------------------------------- 
    3637   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6869 
    6970      ! Depth work arrrays 
    70       ze3t(:,:,:) = e3t(:,:,:,Kmm) 
    71       ze3u(:,:,:) = e3u(:,:,:,Kmm) 
    72       ze3v(:,:,:) = e3v(:,:,:,Kmm) 
    73       ze3w(:,:,:) = e3w(:,:,:,Kmm) 
     71      DO jk = 1 , jpk  
     72         ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 
     73         ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 
     74         ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 
     75         ze3w(:,:,jk) = e3w(:,:,jk,Kmm) 
     76      END DO 
    7477 
    7578      IF( kt == nit000  ) THEN 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/CRS/crsini.F90

    r12377 r13151  
    2828   PUBLIC   crs_init   ! called by nemogcm.F90 module 
    2929 
     30   !! * Substitutions 
     31#  include "domzgr_substitute.h90" 
    3032   !!---------------------------------------------------------------------- 
    3133   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    174176      
    175177     ! 
    176      ze3t(:,:,:) = e3t(:,:,:,Kmm) 
    177      ze3u(:,:,:) = e3u(:,:,:,Kmm) 
    178      ze3v(:,:,:) = e3v(:,:,:,Kmm) 
    179      ze3w(:,:,:) = e3w(:,:,:,Kmm) 
     178     DO jk = 1, jpk 
     179        ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 
     180        ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 
     181        ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 
     182        ze3w(:,:,jk) = e3w(:,:,jk,Kmm) 
     183     END DO   
    180184 
    181185     !    3.d.2   Surfaces  

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diaar5.F90

    r12489 r13151  
    3232   REAL(wp)                         ::   vol0         ! ocean volume (interior domain) 
    3333   REAL(wp)                         ::   area_tot     ! total ocean surface (interior domain) 
    34    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   area         ! cell surface (interior domain) 
    3534   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain) 
    3635   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity 
    … …  
    4039   !! * Substitutions 
    4140#  include "do_loop_substitute.h90" 
     41#  include "domzgr_substitute.h90" 
    4242   !!---------------------------------------------------------------------- 
    4343   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5454      !!---------------------------------------------------------------------- 
    5555      ! 
    56       ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc ) 
     56      ALLOCATE( thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc ) 
    5757      ! 
    5858      CALL mpp_sum ( 'diaar5', dia_ar5_alloc ) 
    … …  
    7777      ! 
    7878      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace  
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zpe, z2d                   ! 2D workspace  
    80       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd , zrhop, ztpot   ! 3D workspace 
     79      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace  
     80      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd , zrhop, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute) 
    8181      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace 
    8282 
    … …  
    9090         ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) ) 
    9191         ALLOCATE( ztsn(jpi,jpj,jpk,jpts) ) 
    92          zarea_ssh(:,:) = area(:,:) * ssh(:,:,Kmm) 
    93       ENDIF 
    94       ! 
    95       CALL iom_put( 'e2u'      , e2u (:,:) ) 
    96       CALL iom_put( 'e1v'      , e1v (:,:) ) 
    97       CALL iom_put( 'areacello', area(:,:) ) 
     92         zarea_ssh(:,:) = e1e2t(:,:) * ssh(:,:,Kmm) 
     93      ENDIF 
     94      ! 
     95      CALL iom_put( 'e2u'      , e2u  (:,:) ) 
     96      CALL iom_put( 'e1v'      , e1v  (:,:) ) 
     97      CALL iom_put( 'areacello', e1e2t(:,:) ) 
    9898      ! 
    9999      IF( iom_use( 'volcello' ) .OR. iom_use( 'masscello' )  ) THEN   
    100100         zrhd(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace 
    101101         DO jk = 1, jpkm1 
    102             zrhd(:,:,jk) = area(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     102            zrhd(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    103103         END DO 
     104         DO jk = 1, jpk 
     105            z3d(:,:,jk) =  rho0 * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     106         END DO  
    104107         CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000 
    105          CALL iom_put( 'masscello' , rho0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass 
     108         CALL iom_put( 'masscello' , z3d (:,:,:) )   ! ocean mass 
    106109      ENDIF  
    107110      ! 
    … …  
    129132         ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh 
    130133         ztsn(:,:,:,jp_sal) = sn0(:,:,:) 
    131          CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity 
     134         DO jk = 1, jpk 
     135            zgdept(:,:,jk) = gdept(:,:,jk,Kmm) 
     136         END DO 
     137         CALL eos( ztsn, zrhd, zgdept)                       ! now in situ density using initial salinity 
    132138         ! 
    133139         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    … …  
    151157         END IF 
    152158         !                                          
    153          zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) )  
     159         zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) )  
    154160         zssh_steric = - zarho / area_tot 
    155161         CALL iom_put( 'sshthster', zssh_steric ) 
    156162       
    157163         !                                         ! steric sea surface height 
    158          CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density 
     164         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, zgdept )                 ! now in situ and potential density 
    159165         zrhop(:,:,jpk) = 0._wp 
    160166         CALL iom_put( 'rhop', zrhop ) 
    … …  
    177183         END IF 
    178184         !     
    179          zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) )  
     185         zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) )  
    180186         zssh_steric = - zarho / area_tot 
    181187         CALL iom_put( 'sshsteric', zssh_steric ) 
    … …  
    191197          ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity 
    192198          DO_3D_11_11( 1, jpkm1 ) 
    193              zztmp = area(ji,jj) * e3t(ji,jj,jk,Kmm) 
     199             zztmp = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) 
    194200             ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm) 
    195201             ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * ts(ji,jj,jk,jp_sal,Kmm) 
    … …  
    237243               z2d(:,:) = 0._wp 
    238244               DO jk = 1, jpkm1 
    239                  z2d(:,:) = z2d(:,:) + area(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk) 
     245                 z2d(:,:) = z2d(:,:) + e1e2t(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk) 
    240246               END DO 
    241247               ztemp = glob_sum( 'diaar5', z2d(:,:)  )  
    … …  
    244250             ! 
    245251             IF( iom_use( 'ssttot' ) ) THEN   ! Output potential temperature in case we use TEOS-10 
    246                zsst = glob_sum( 'diaar5',  area(:,:) * ztpot(:,:,1)  )  
     252               zsst = glob_sum( 'diaar5',  e1e2t(:,:) * ztpot(:,:,1)  )  
    247253               CALL iom_put( 'ssttot', zsst / area_tot ) 
    248254             ENDIF 
    … …  
    259265      ELSE        
    260266         IF( iom_use('ssttot') ) THEN   ! Output sst in case we use EOS-80 
    261             zsst  = glob_sum( 'diaar5', area(:,:) * ts(:,:,1,jp_tem,Kmm) ) 
     267            zsst  = glob_sum( 'diaar5', e1e2t(:,:) * ts(:,:,1,jp_tem,Kmm) ) 
    262268            CALL iom_put('ssttot', zsst / area_tot ) 
    263269         ENDIF 
    … …  
    375381         IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' ) 
    376382 
    377          area(:,:) = e1e2t(:,:) 
    378          area_tot  = glob_sum( 'diaar5', area(:,:) ) 
     383         area_tot  = glob_sum( 'diaar5', e1e2t(:,:) ) 
    379384 
    380385         ALLOCATE( zvol0(jpi,jpj) ) 
    … …  
    383388         DO_3D_11_11( 1, jpkm1 ) 
    384389            idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk) 
    385             zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj) 
     390            zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * e1e2t(ji,jj) 
    386391            thick0(ji,jj) = thick0(ji,jj) +  idep     
    387392         END_3D 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diacfl.F90

    r12489 r13151  
    3434   !! * Substitutions 
    3535#  include "do_loop_substitute.h90" 
     36#  include "domzgr_substitute.h90" 
    3637   !!---------------------------------------------------------------------- 
    3738   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diadct.F90

    r12489 r13151  
    1111   !!            3.4  ! 09/2011 (C Bricaud) 
    1212   !!---------------------------------------------------------------------- 
    13    !! does not work with agrif 
    14 #if ! defined key_agrif 
     13#if ! defined key_agrif        
     14   !!                        ==>>  CAUTION: does not work with agrif 
    1515   !!---------------------------------------------------------------------- 
    1616   !!   dia_dct      :  Compute the transport through a sec. 
    … …  
    6666   TYPE SECTION 
    6767      CHARACTER(len=60)                            :: name              ! name of the sec 
    68       LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and 
    69                                                                        ! heat transports 
     68      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and heat transports 
    7069      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation 
    7170      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line 
    … …  
    7473      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section 
    7574      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class 
    76       REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want) 
    77                                                       zsigp           ,&! potential density classes    (99 if you don't want) 
    78                                                       zsal            ,&! salinity classes   (99 if you don't want) 
    79                                                       ztem            ,&! temperature classes(99 if you don't want) 
    80                                                       zlay              ! level classes      (99 if you don't want) 
     75      REAL(wp), DIMENSION(nb_class_max)            :: zsigi             ! in-situ   density classes    (99 if you don't want) 
     76      REAL(wp), DIMENSION(nb_class_max)            :: zsigp             ! potential density classes    (99 if you don't want) 
     77      REAL(wp), DIMENSION(nb_class_max)            :: zsal              ! salinity classes   (99 if you don't want) 
     78      REAL(wp), DIMENSION(nb_class_max)            :: ztem              ! temperature classes(99 if you don't want) 
     79      REAL(wp), DIMENSION(nb_class_max)            :: zlay              ! level classes      (99 if you don't want) 
    8180      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
    8281      REAL(wp)                                         :: slopeSection  ! slope of the section 
    … …  
    9089   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
    9190 
     91 
     92   !! * Substitutions 
     93#  include "domzgr_substitute.h90" 
    9294   !!---------------------------------------------------------------------- 
    9395   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    9597   !! Software governed by the CeCILL license (see ./LICENSE) 
    9698   !!---------------------------------------------------------------------- 
     99 
    97100CONTAINS 
    98101  
    … …  
    11191122  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1) 
    11201123  !!    |               |                  |       zbis =  
    1121   !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ] 
    1122   !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]  
     1124  !!    |               |                  |      [ e3w_n(I+1,J,K,NOW)*ptab(I,J,K) + ( e3w_n(I,J,K,NOW) - e3w_n(I+1,J,K,NOW) ) * ptab(I,J,K-1) ] 
     1125  !!    |               |                  |     /[ e3w_n(I+1,J,K,NOW)             +   e3w_n(I,J,K,NOW) - e3w_n(I+1,J,K,NOW) ]  
    11231126  !!    |               |                  |  
    11241127  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point 
    … …  
    12121215  ELSE       ! full step or partial step case  
    12131216 
    1214      ze3t  = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm)  
    1215      zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm) 
    1216      zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm) 
     1217     ze3t  =   e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm)  
     1218     zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) )   & 
     1219        &    / e3w(ii2,ij2,kk,Kmm) 
     1220     zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) )   & 
     1221        &    / e3w(ii1,ij1,kk,Kmm) 
    12171222 
    12181223     IF(kk .NE. 1)THEN 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diahsb.F90

    r12489 r13151  
    5050   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini 
    5151 
     52   !! * Substitutions 
     53#  include "domzgr_substitute.h90" 
    5254   !!---------------------------------------------------------------------- 
    5355   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    156158      ! 
    157159      DO jk = 1, jpkm1           ! volume variation (calculated with scale factors) 
    158          zwrk(:,:,jk) = surf(:,:)*e3t(:,:,jk,Kmm)*tmask(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk)*tmask_ini(:,:,jk) 
     160         zwrk(:,:,jk) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   & 
     161            &           - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk) 
    159162      END DO 
    160163      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different 
    161164      DO jk = 1, jpkm1           ! heat content variation 
    162          zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) 
     165         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   & 
     166            &           - surf_ini(:,:) *         hc_loc_ini(:,:,jk) ) 
    163167      END DO 
    164168      zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
    165169      DO jk = 1, jpkm1           ! salt content variation 
    166          zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) 
     170         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   & 
     171            &           - surf_ini(:,:) *         sc_loc_ini(:,:,jk) ) 
    167172      END DO 
    168173      zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
    … …  
    287292            DO jk = 1, jpk 
    288293              ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance). 
    289                e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors 
     294               e3t_ini   (:,:,jk) =                       e3t(:,:,jk,Kmm)*tmask(:,:,jk)  ! initial vertical scale factors 
    290295               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask 
    291                hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content 
    292                sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content 
     296               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm)*e3t(:,:,jk,Kmm)*tmask(:,:,jk)  ! initial heat content 
     297               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm)*e3t(:,:,jk,Kmm)*tmask(:,:,jk)  ! initial salt content 
    293298            END DO 
    294299            frc_v = 0._wp                                           ! volume       trend due to forcing 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diahth.F90

    r12489 r13151  
    4242   !! * Substitutions 
    4343#  include "do_loop_substitute.h90" 
     44#  include "domzgr_substitute.h90" 
    4445   !!---------------------------------------------------------------------- 
    4546   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    361362         ik = ilevel(ji,jj) 
    362363         zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep 
    363          phtc(ji,jj)   = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 
     364         phtc(ji,jj)   = phtc(ji,jj)    & 
     365            &           + pt (ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 
    364366                                                       * tmask(ji,jj,ik+1) 
    365367      END_2D 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diamlr.F90

    r12377 r13151  
    44   !! Management of the IOM context for multiple-linear-regression analysis 
    55   !!====================================================================== 
    6    !! History :       !  2019  (S. Mueller) 
     6   !! History :  4.0  !  2019  (S. Mueller)   Original code 
    77   !!---------------------------------------------------------------------- 
    88 
    99   USE par_oce        , ONLY :   wp, jpi, jpj 
    1010   USE phycst         , ONLY :   rpi 
     11   USE dom_oce        , ONLY :   adatrj 
     12   USE tide_mod 
     13   ! 
    1114   USE in_out_manager , ONLY :   lwp, numout, ln_timing 
    1215   USE iom            , ONLY :   iom_put, iom_use, iom_update_file_name 
    13    USE dom_oce        , ONLY :   adatrj 
    1416   USE timing         , ONLY :   timing_start, timing_stop 
    1517#if defined key_iomput 
    1618   USE xios 
    1719#endif 
    18    USE tide_mod 
    1920 
    2021   IMPLICIT NONE 
    2122   PRIVATE 
    2223 
    23    LOGICAL, PUBLIC ::   lk_diamlr = .FALSE. 
     24   LOGICAL, PUBLIC ::   lk_diamlr = .FALSE.   !:         ===>>>   NOT a DOCTOR norm name :  use l_diamlr 
     25   !                                                              lk_  is used only for logical controlled by a CPP key 
    2426 
    2527   PUBLIC ::   dia_mlr_init, dia_mlr_iom_init, dia_mlr 
    … …  
    3335   !!---------------------------------------------------------------------- 
    3436CONTAINS 
    35     
     37 
    3638   SUBROUTINE dia_mlr_init 
    3739      !!---------------------------------------------------------------------- 
    3840      !!                 ***  ROUTINE dia_mlr_init  *** 
    3941      !! 
    40       !! ** Purpose : initialisation of IOM context management for  
     42      !! ** Purpose : initialisation of IOM context management for 
    4143      !!              multiple-linear-regression analysis 
    4244      !! 
    4345      !!---------------------------------------------------------------------- 
    44  
     46      ! 
    4547      lk_diamlr = .TRUE. 
    46  
     48      ! 
    4749      IF(lwp) THEN 
    4850         WRITE(numout, *) 
    … …  
    5052         WRITE(numout, *) '~~~~~~~~~~~~   multiple-linear-regression analysis' 
    5153      END IF 
    52  
     54      ! 
    5355   END SUBROUTINE dia_mlr_init 
     56 
    5457 
    5558   SUBROUTINE dia_mlr_iom_init 
    … …  
    8487      INTEGER                                     ::   itide                       ! Number of available tidal components 
    8588      REAL(wp)                                    ::   ztide_phase                 ! Tidal-constituent phase at adatrj=0 
    86       CHARACTER (LEN=4), DIMENSION(jpmax_harmo)   ::   ctide_selected = ' n/a ' 
     89      CHARACTER (LEN=4), DIMENSION(jpmax_harmo)   ::   ctide_selected = 'n/a ' 
    8790      TYPE(tide_harmonic), DIMENSION(:), POINTER  ::   stideconst 
    8891 
    … …  
    145148            ! Retrieve information (frequency, phase, nodal correction) about all 
    146149            ! available tidal constituents for placeholder substitution below 
    147             ctide_selected(1:34) = (/ 'Mf', 'Mm', 'Ssa', 'Mtm', 'Msf',    & 
    148                &                      'Msqm', 'Sa', 'K1', 'O1', 'P1',     & 
    149                &                      'Q1', 'J1', 'S1', 'M2', 'S2', 'N2', & 
    150                &                      'K2', 'nu2', 'mu2', '2N2', 'L2',    & 
    151                &                      'T2', 'eps2', 'lam2', 'R2', 'M3',   & 
    152                &                      'MKS2', 'MN4', 'MS4', 'M4', 'N4',   & 
    153                &                      'S4', 'M6', 'M8' /) 
     150            ! Warning: we must use the same character length in an array constructor (at least for gcc compiler) 
     151            ctide_selected(1:34) = (/ 'Mf  ', 'Mm  ', 'Ssa ', 'Mtm ', 'Msf ',         & 
     152               &                      'Msqm', 'Sa  ', 'K1  ', 'O1  ', 'P1  ',         & 
     153               &                      'Q1  ', 'J1  ', 'S1  ', 'M2  ', 'S2  ', 'N2  ', & 
     154               &                      'K2  ', 'nu2 ', 'mu2 ', '2N2 ', 'L2  ',         & 
     155               &                      'T2  ', 'eps2', 'lam2', 'R2  ', 'M3  ',         & 
     156               &                      'MKS2', 'MN4 ', 'MS4 ', 'M4  ', 'N4  ',         & 
     157               &                      'S4  ', 'M6  ', 'M8  ' /) 
    154158            CALL tide_init_harmonics(ctide_selected, stideconst) 
    155159            itide = size(stideconst) 
    … …  
    157161            itide = 0 
    158162         ENDIF 
    159           
     163 
    160164         DO jm = 1, jpscanmax 
    161165            WRITE (cl3i, '(i3.3)') jm 
    … …  
    236240               ! If enabled, keep handle in list of fields selected for analysis 
    237241               IF ( llxatt_enabled ) THEN 
    238                    
     242 
    239243                  ! Set name attribute (and overwrite possible pre-configured name) 
    240244                  ! with field id to enable id string retrieval from stored handle 
    … …  
    323327            CALL xios_set_attr  ( slxhdl_fld, standard_name=TRIM( clxatt_comment ), long_name=TRIM( clxatt_expr ),   & 
    324328               &                  operation="average" ) 
    325                 
     329 
    326330            ! iii) set up the output of scalar products with itself and with 
    327331            !      other active regressors 
    … …  
    396400   END SUBROUTINE dia_mlr_iom_init 
    397401 
     402 
    398403   SUBROUTINE dia_mlr 
    399404      !!---------------------------------------------------------------------- 
    … …  
    403408      !! 
    404409      !!---------------------------------------------------------------------- 
    405  
    406410      REAL(wp), DIMENSION(jpi,jpj) ::   zadatrj2d 
     411      !!---------------------------------------------------------------------- 
    407412 
    408413      IF( ln_timing )   CALL timing_start('dia_mlr') 
    … …  
    411416      ! (value of adatrj converted to time in units of seconds) 
    412417      ! 
    413       ! A 2-dimensional field of constant value is sent, and subsequently used 
    414       ! directly or transformed to a scalar or a constant 3-dimensional field as 
    415       ! required. 
     418      ! A 2-dimensional field of constant value is sent, and subsequently used directly  
     419      ! or transformed to a scalar or a constant 3-dimensional field as required. 
    416420      zadatrj2d(:,:) = adatrj*86400.0_wp 
    417421      IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d) 
    418        
     422      ! 
    419423      IF( ln_timing )   CALL timing_stop('dia_mlr') 
    420  
     424      ! 
    421425   END SUBROUTINE dia_mlr 
    422426 
     427   !!====================================================================== 
    423428END MODULE diamlr 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diaptr.F90

    r12489 r13151  
    6060 
    6161   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini) 
     62    
    6263   !! * Substitutions 
    6364#  include "do_loop_substitute.h90" 
     65#  include "domzgr_substitute.h90" 
    6466   !!---------------------------------------------------------------------- 
    6567   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diawri.F90

    r12493 r13151  
    8585   !! * Substitutions 
    8686#  include "do_loop_substitute.h90" 
     87#  include "domzgr_substitute.h90" 
    8788   !!---------------------------------------------------------------------- 
    8889   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    136137      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
    137138      ! 
    138       CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 
    139       CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 
    140       CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 
    141       CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 
    142       IF( iom_use("e3tdef") )   & 
    143          CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    144  
    145       IF( ll_wd ) THEN 
    146          CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
     139      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t 
     140         DO jk = 1, jpk 
     141            z3d(:,:,jk) =  e3t(:,:,jk,Kmm) 
     142         END DO 
     143         CALL iom_put( "e3t"     ,     z3d(:,:,:) ) 
     144         CALL iom_put( "e3tdef"  , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )  
     145      ENDIF  
     146      IF ( iom_use("e3u") ) THEN                         ! time-varying e3u 
     147         DO jk = 1, jpk 
     148            z3d(:,:,jk) =  e3u(:,:,jk,Kmm) 
     149         END DO  
     150         CALL iom_put( "e3u" , z3d(:,:,:) ) 
     151      ENDIF 
     152      IF ( iom_use("e3v") ) THEN                         ! time-varying e3v 
     153         DO jk = 1, jpk 
     154            z3d(:,:,jk) =  e3v(:,:,jk,Kmm) 
     155         END DO  
     156         CALL iom_put( "e3v" , z3d(:,:,:) ) 
     157      ENDIF 
     158      IF ( iom_use("e3w") ) THEN                         ! time-varying e3w 
     159         DO jk = 1, jpk 
     160            z3d(:,:,jk) =  e3w(:,:,jk,Kmm) 
     161         END DO  
     162         CALL iom_put( "e3w" , z3d(:,:,:) ) 
     163      ENDIF 
     164 
     165      IF( ll_wd ) THEN                                   ! sea surface height (brought back to the reference used for wetting and drying) 
     166         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) 
    147167      ELSE 
    148168         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    … …  
    208228 
    209229      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
    210       ! 
    211230      CALL iom_put( "woce", ww )                   ! vertical velocity 
     231 
    212232      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
    213233         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    … …  
    415435      ! 
    416436      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    417       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     437      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept       ! 3D workspace 
    418438      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace 
    419439      !!---------------------------------------------------------------------- 
    … …  
    455475      it = kt 
    456476      itmod = kt - nit000 + 1 
     477 
     478      ! store e3t for subsitute 
     479      DO jk = 1, jpk 
     480         ze3t  (:,:,jk) =  e3t  (:,:,jk,Kmm) 
     481         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm) 
     482      END DO 
    457483 
    458484 
    … …  
    569595         DEALLOCATE(zw3d_abl) 
    570596         ENDIF 
     597         ! 
    571598 
    572599         ! Declare all the output fields as NETCDF variables 
    … …  
    578605            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    579606         IF(  .NOT.ln_linssh  ) THEN 
    580             CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm) 
     607            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t n 
    581608            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    582             CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm) 
     609            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t n 
    583610            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    584             CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm) 
     611            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t n 
    585612            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    586613         ENDIF 
    … …  
    766793 
    767794      IF( .NOT.ln_linssh ) THEN 
    768          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content 
    769          CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content 
    770          CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content 
    771          CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content 
     795         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content 
     796         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content 
     797         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content 
     798         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content 
    772799      ELSE 
    773800         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature 
    … …  
    777804      ENDIF 
    778805      IF( .NOT.ln_linssh ) THEN 
    779          zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    780          CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness 
    781          CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth 
     806         zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
     807         CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)    , ndim_T , ndex_T  )   ! level thickness 
     808         CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T  )   ! t-point depth  
    782809         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation 
    783810      ENDIF 
    … …  
    918945      !! 
    919946      INTEGER :: inum, jk 
     947      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept      ! 3D workspace !!st patch to use substitution 
    920948      !!---------------------------------------------------------------------- 
    921949      !  
    922       IF(lwp) WRITE(numout,*) 
    923       IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state' 
    924       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
    925       IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
    926  
    927 #if defined key_si3 
    928      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
    929 #else 
    930      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
    931 #endif 
    932  
     950      IF(lwp) THEN 
     951         WRITE(numout,*) 
     952         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state' 
     953         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
     954         WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
     955      ENDIF  
     956      ! 
     957      DO jk = 1, jpk 
     958         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm) 
     959         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm) 
     960      END DO 
     961      ! 
     962      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
     963      ! 
    933964      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature 
    934965      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity 
    935       CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height 
    936       CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity 
    937       CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity 
     966      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:        ,Kmm) )    ! sea surface height 
     967      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:       ,Kmm) )    ! now i-velocity 
     968      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:       ,Kmm) )    ! now j-velocity 
    938969      IF( ln_zad_Aimp ) THEN 
    939970         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity 
    … …  
    942973      ENDIF 
    943974      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity 
    944       CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height 
    945  
     975      CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height 
     976      ! 
    946977      IF ( ln_isf ) THEN 
    947978         IF (ln_isfcav_mlt) THEN 
    … …  
    949980            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity 
    950981            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity 
    951             CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity 
    952             CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity 
    953             CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 ) 
     982            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity 
     983            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity 
     984            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 
    954985         END IF 
    955986         IF (ln_isfpar_mlt) THEN 
    956             CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,8) )    ! now k-velocity 
     987            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity 
    957988            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity 
    958989            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity 
    959990            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity 
    960             CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity 
    961             CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity 
    962             CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 ) 
     991            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity 
     992            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity 
     993            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 
    963994         END IF 
    964995      END IF 
    965  
     996      ! 
    966997      IF( ALLOCATED(ahtu) ) THEN 
    967998         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
    … …  
    9781009      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress 
    9791010      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress 
    980       IF(  .NOT.ln_linssh  ) THEN              
    981          CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth  
    982          CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness   
     1011      IF(  .NOT.ln_linssh  ) THEN 
     1012         CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept        )    !  T-cell depth  
     1013         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t          )    !  T-cell thickness   
    9831014      END IF 
    9841015      IF( ln_wave .AND. ln_sdw ) THEN 
    … …  
    9931024         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity 
    9941025      ENDIF 
    995   
     1026      ! 
     1027      CALL iom_close( inum ) 
     1028      !  
    9961029#if defined key_si3 
    9971030      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid 
     1031         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 
    9981032         CALL ice_wri_state( inum ) 
     1033         CALL iom_close( inum ) 
    9991034      ENDIF 
    10001035#endif 
    1001       ! 
    1002       CALL iom_close( inum ) 
    1003       !  
     1036 
    10041037   END SUBROUTINE dia_wri_state 
    10051038 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/dom_oce.F90

    r12489 r13151  
    22   !!====================================================================== 
    33   !!                       ***  MODULE dom_oce  *** 
    4    !!        
    54   !! ** Purpose :   Define in memory all the ocean space domain variables 
    65   !!====================================================================== 
    7    !! History :  1.0  ! 2005-10  (A. Beckmann, G. Madec)  reactivate s-coordinate  
     6   !! History :  1.0  ! 2005-10  (A. Beckmann, G. Madec)  reactivate s-coordinate 
    87   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
    98   !!            3.4  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
    … …  
    1312   !!             -   ! 2015-11  (G. Madec, A. Coward)  time varying zgr by default 
    1413   !!            4.1  ! 2019-08  (A. Coward, D. Storkey) rename prognostic variables in preparation for new time scheme. 
     14   !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
    1515   !!---------------------------------------------------------------------- 
    1616 
    … …  
    7171   !                                !  = 6 cyclic East-West AND North fold F-point pivot 
    7272   !                                !  = 7 bi-cyclic East-West AND North-South 
    73    LOGICAL, PUBLIC ::   l_Iperio, l_Jperio   !   should we explicitely take care I/J periodicity  
    74  
    75    !                                 ! domain MPP decomposition parameters 
     73   LOGICAL, PUBLIC ::   l_Iperio, l_Jperio   !   should we explicitely take care I/J periodicity 
     74 
     75   !                             !: domain MPP decomposition parameters 
    7676   INTEGER             , PUBLIC ::   nimpp, njmpp     !: i- & j-indexes for mpp-subdomain left bottom 
    7777   INTEGER             , PUBLIC ::   nreci, nrecj     !: overlap region in i and j 
    … …  
    8181   INTEGER, ALLOCATABLE, PUBLIC ::   nbondi_bdy(:)    !: mark i-direction local boundaries for BDY open boundaries 
    8282   INTEGER, ALLOCATABLE, PUBLIC ::   nbondj_bdy(:)    !: mark j-direction local boundaries for BDY open boundaries 
    83    INTEGER, ALLOCATABLE, PUBLIC ::   nbondi_bdy_b(:)  !: mark i-direction of neighbours local boundaries for BDY open boundaries   
    84    INTEGER, ALLOCATABLE, PUBLIC ::   nbondj_bdy_b(:)  !: mark j-direction of neighbours local boundaries for BDY open boundaries   
     83   INTEGER, ALLOCATABLE, PUBLIC ::   nbondi_bdy_b(:)  !: mark i-direction of neighbours local boundaries for BDY open boundaries 
     84   INTEGER, ALLOCATABLE, PUBLIC ::   nbondj_bdy_b(:)  !: mark j-direction of neighbours local boundaries for BDY open boundaries 
    8585 
    8686   INTEGER, PUBLIC ::   npolj             !: north fold mark (0, 3 or 4) 
    … …  
    126126   LOGICAL, PUBLIC ::   ln_zps       !: z-coordinate - partial step 
    127127   LOGICAL, PUBLIC ::   ln_sco       !: s-coordinate or hybrid z-s coordinate 
    128    LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF  
     128   LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF 
    129129   !                                                        !  reference scale factors 
    130130   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0   !: t- vert. scale factor [m] 
    … …  
    136136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0   !: vw-vert. scale factor [m] 
    137137   !                                                        !  time-dependent scale factors 
     138#if ! defined key_qco 
    138139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e3t, e3u, e3v, e3w, e3uw, e3vw  !: vert. scale factor [m] 
    139140   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   e3f                             !: F-point vert. scale factor [m] 
     141#endif 
     142   !                                                        !  time-dependent ratio ssh / h_0 
     143   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   r3t, r3u, r3v                   !: time-dependent    ratio at t-, u- and v-point [-] 
     144   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3f                             !: mid-time-level    ratio at f-point            [-] 
     145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   r3t_f, r3u_f, r3v_f             !: now time-filtered ratio at t-, u- and v-point [-] 
    140146 
    141147   !                                                        !  reference depths of cells 
    142    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0  !: t- depth              [m] 
    143    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0  !: w- depth              [m] 
    144    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0  !: w- depth (sum of e3w) [m] 
     148   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdept_0  !: t- depth              [m] 
     149   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gdepw_0  !: w- depth              [m] 
     150   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w_0  !: w- depth (sum of e3w) [m] 
    145151   !                                                        !  time-dependent depths of cells 
    146    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  gdept, gdepw   
    147    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  gde3w   
    148     
    149    !                                                      !  reference heights of water column 
    150    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m] 
    151    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m] 
    152    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m] 
    153                                                           ! time-dependent heights of water column 
    154    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht                     !: height of water column at T-points [m] 
    155    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hu, hv, r1_hu, r1_hv   !: height of water column [m] and reciprocal [1/m] 
     152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   gdept, gdepw 
     153   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   gde3w 
     154 
     155   !                                                        !  reference heights of ocean water column and its inverse 
     156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht_0, r1_ht_0   !: t-depth        [m] and [1/m] 
     157   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hu_0, r1_hu_0   !: u-depth        [m] and [1/m] 
     158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hv_0, r1_hv_0   !: v-depth        [m] and [1/m] 
     159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hf_0, r1_hf_0   !: f-depth        [m] and [1/m] 
     160   !                                                        ! time-dependent heights of ocean water column 
     161#if ! defined key_qco 
     162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ht          !: t-points           [m] 
     163#endif 
     164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hu, r1_hu   !: u-depth            [m] and [1/m] 
     165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hv, r1_hv   !: v-depth            [m] and [1/m] 
    156166 
    157167   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
    158    INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1)  
     168   INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1) 
    159169 
    160170   !! 1D reference  vertical coordinate 
    … …  
    169179   !! --------------------------------------------------------------------- 
    170180!!gm Proposition of new name for top/bottom vertical indices 
    171 !   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mtk_t, mtk_u, mtk_v   !: top first wet T-, U-, V-, F-level (ISF) 
    172 !   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbk_t, mbk_u, mbk_v   !: bottom last wet T-, U- and V-level 
     181!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mtk_t, mtk_u, mtk_v   !: top    first wet T-, U-, and V-level (ISF) 
     182!   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbk_t, mbk_u, mbk_v   !: bottom last  wet T-, U-, and V-level 
    173183!!gm 
    174184   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt, mbku, mbkv   !: bottom last wet T-, U- and V-level 
    … …  
    178188   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level           (ISF) 
    179189 
    180    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts 
    181    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    182    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
     190   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   ssmask, ssumask, ssvmask, ssfmask   !: surface mask at T-,U-, V- and F-pts 
     191   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask            !: land/ocean mask at T-, U-, V- and F-pts 
     192   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask                 !: land/ocean mask at WT-, WU- and WV-pts 
    183193 
    184194   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    … …  
    198208   INTEGER , PUBLIC ::   nsec_monday   !: seconds between 00h         of the last Monday   and half of the current time step 
    199209   INTEGER , PUBLIC ::   nsec_day      !: seconds between 00h         of the current   day and half of the current time step 
    200    REAL(wp), PUBLIC ::   fjulday       !: current julian day  
     210   REAL(wp), PUBLIC ::   fjulday       !: current julian day 
    201211   REAL(wp), PUBLIC ::   fjulstartyear !: first day of the current year in julian days 
    202212   REAL(wp), PUBLIC ::   adatrj        !: number of elapsed days since the begining of the whole simulation 
    203213   !                                   !: (cumulative duration of previous runs that may have used different time-step size) 
    204    INTEGER , PUBLIC, DIMENSION(  0: 2) ::   nyear_len     !: length in days of the previous/current/next year 
    205    INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_len    !: length in days of the months of the current year 
    206    INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_beg    !: second since Jan 1st 0h of the current year and the half of the months 
    207    INTEGER , PUBLIC                  ::   nsec1jan000     !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year 
    208    INTEGER , PUBLIC                  ::   nsec000_1jan000   !: second since Jan 1st 0h of nit000 year and nit000 
    209    INTEGER , PUBLIC                  ::   nsecend_1jan000   !: second since Jan 1st 0h of nit000 year and nitend 
     214   INTEGER , PUBLIC, DIMENSION(  0: 2) ::   nyear_len         !: length in days of the previous/current/next year 
     215   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_len        !: length in days of the months of the current year 
     216   INTEGER , PUBLIC, DIMENSION(-11:25) ::   nmonth_beg        !: second since Jan 1st 0h of the current year and the half of the months 
     217   INTEGER , PUBLIC                    ::   nsec1jan000       !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year 
     218   INTEGER , PUBLIC                    ::   nsec000_1jan000   !: second since Jan 1st 0h of nit000 year and nit000 
     219   INTEGER , PUBLIC                    ::   nsecend_1jan000   !: second since Jan 1st 0h of nit000 year and nitend 
    210220 
    211221   !!---------------------------------------------------------------------- 
    … …  
    220230   !!---------------------------------------------------------------------- 
    221231   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    222    !! $Id$  
     232   !! $Id$ 
    223233   !! Software governed by the CeCILL license (see ./LICENSE) 
    224234   !!---------------------------------------------------------------------- 
    … …  
    234244 
    235245   CHARACTER(len=3) FUNCTION Agrif_CFixed() 
    236       Agrif_CFixed = '0'  
     246      Agrif_CFixed = '0' 
    237247   END FUNCTION Agrif_CFixed 
    238248#endif 
    … …  
    240250   INTEGER FUNCTION dom_oce_alloc() 
    241251      !!---------------------------------------------------------------------- 
    242       INTEGER, DIMENSION(12) :: ierr 
     252      INTEGER                ::   ii 
     253      INTEGER, DIMENSION(30) :: ierr 
    243254      !!---------------------------------------------------------------------- 
    244       ierr(:) = 0 
     255      ii = 0   ;   ierr(:) = 0 
    245256      ! 
    246       ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(1) ) 
    247          ! 
    248       ALLOCATE( mi0(jpiglo)   , mi1 (jpiglo),  mj0(jpjglo)   , mj1 (jpjglo) ,     & 
    249          &      tpol(jpiglo)  , fpol(jpiglo)                                , STAT=ierr(2) ) 
    250          ! 
     257      ii = ii+1 
     258      ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(ii) ) 
     259         ! 
     260      ii = ii+1 
     261      ALLOCATE( mi0 (jpiglo) , mi1 (jpiglo),  mj0(jpjglo) , mj1 (jpjglo) ,     & 
     262         &      tpol(jpiglo) , fpol(jpiglo)                              , STAT=ierr(ii) ) 
     263         ! 
     264      ii = ii+1 
    251265      ALLOCATE( glamt(jpi,jpj) ,    glamu(jpi,jpj) ,  glamv(jpi,jpj) ,  glamf(jpi,jpj) ,     & 
    252266         &      gphit(jpi,jpj) ,    gphiu(jpi,jpj) ,  gphiv(jpi,jpj) ,  gphif(jpi,jpj) ,     & 
    … …  
    259273         &      e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj)                   ,     & 
    260274         &      e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj)                                     ,     & 
    261          &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) ) 
    262          ! 
     275         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(ii) ) 
     276         ! 
     277      ii = ii+1 
    263278      ALLOCATE( gdept_0(jpi,jpj,jpk)     , gdepw_0(jpi,jpj,jpk)     , gde3w_0(jpi,jpj,jpk) ,      & 
    264          &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(4) ) 
    265          ! 
    266       ALLOCATE( e3t_0(jpi,jpj,jpk)     , e3u_0(jpi,jpj,jpk)     , e3v_0(jpi,jpj,jpk)     , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk)     ,   & 
    267          &      e3t  (jpi,jpj,jpk,jpt) , e3u  (jpi,jpj,jpk,jpt) , e3v  (jpi,jpj,jpk,jpt) , e3f  (jpi,jpj,jpk) , e3w  (jpi,jpj,jpk,jpt) ,   &  
    268          &      e3uw_0(jpi,jpj,jpk)     , e3vw_0(jpi,jpj,jpk)     ,         & 
    269          &      e3uw  (jpi,jpj,jpk,jpt) , e3vw  (jpi,jpj,jpk,jpt) ,    STAT=ierr(5) )                        
    270          ! 
    271       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj)    , hv_0(jpi,jpj)     ,                                             & 
    272          &      ht  (jpi,jpj) , hu(  jpi,jpj,jpt), hv(  jpi,jpj,jpt) , r1_hu(jpi,jpj,jpt) , r1_hv(jpi,jpj,jpt) ,   & 
    273          &                      STAT=ierr(6)  ) 
    274          ! 
    275       ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  )  
    276          ! 
    277       ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) ) 
    278          ! 
    279       ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &  
    280          &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     & 
    281          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
    282          ! 
    283       ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) 
    284          ! 
    285       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &  
    286          &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) ) 
    287          ! 
    288       ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
     279         &      gdept  (jpi,jpj,jpk,jpt) , gdepw  (jpi,jpj,jpk,jpt) , gde3w  (jpi,jpj,jpk) , STAT=ierr(ii) ) 
     280         ! 
     281      ii = ii+1 
     282      ALLOCATE(  e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , e3v_0 (jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) ,      & 
     283         &       e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk)                      ,  STAT=ierr(ii) ) 
     284         ! 
     285#if ! defined key_qco 
     286      ii = ii+1 
     287      ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      & 
     288         &      e3w(jpi,jpj,jpk,jpt) , e3uw(jpi,jpj,jpk,jpt) , e3vw(jpi,jpj,jpk,jpt)                    ,  STAT=ierr(ii) ) 
     289#endif   
     290         ! 
     291      ii = ii+1 
     292      ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,  & 
     293         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )        
     294         ! 
     295      ii = ii+1 
     296      ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)    ,    hv_0(jpi,jpj)     , hf_0(jpi,jpj) ,       & 
     297         &   r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) ,    r1_hv_0(jpi,jpj),   r1_hf_0(jpi,jpj) ,   STAT=ierr(ii)  ) 
     298         ! 
     299#if ! defined key_qco 
     300      ii = ii+1 
     301      ALLOCATE( ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       & 
     302         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  ) 
     303#else 
     304      ii = ii+1 
     305      ALLOCATE(                    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       & 
     306         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  ) 
     307#endif 
     308         ! 
     309      ii = ii+1 
     310      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii)  )  
     311         ! 
     312      ii = ii+1 
     313      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(ii) ) 
     314         ! 
     315      ii = ii+1 
     316      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                           & 
     317         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) ,     & 
     318         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) ,                    STAT=ierr(ii) ) 
     319         ! 
     320      ii = ii+1 
     321      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(ii) ) 
     322         ! 
     323      ii = ii+1 
     324      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
     325         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(ii) ) 
     326         ! 
     327      ii = ii+1 
     328      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(ii) ) 
    289329      ! 
    290330      dom_oce_alloc = MAXVAL(ierr) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/domain.F90

    r12489 r13151  
    66   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code 
    77   !!                 !  1992-01  (M. Imbard) insert time step initialization 
    8    !!                 !  1996-06  (G. Madec) generalized vertical coordinate  
     8   !!                 !  1996-06  (G. Madec) generalized vertical coordinate 
    99   !!                 !  1997-02  (G. Madec) creation of domwri.F 
    1010   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea 
    … …  
    1515   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default 
    1616   !!            4.0  !  2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface 
     17   !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
    1718   !!---------------------------------------------------------------------- 
    18     
     19 
    1920   !!---------------------------------------------------------------------- 
    2021   !!   dom_init      : initialize the space and time domain 
    … …  
    3435   USE dommsk         ! domain: set the mask system 
    3536   USE domwri         ! domain: write the meshmask file 
     37#if ! defined key_qco 
    3638   USE domvvl         ! variable volume 
     39#else 
     40   USE domqco          ! variable volume 
     41#endif 
    3742   USE c1d            ! 1D configuration 
    3843   USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine) 
    … …  
    6166      !!---------------------------------------------------------------------- 
    6267      !!                  ***  ROUTINE dom_init  *** 
    63       !!                     
    64       !! ** Purpose :   Domain initialization. Call the routines that are  
    65       !!              required to create the arrays which define the space  
     68      !! 
     69      !! ** Purpose :   Domain initialization. Call the routines that are 
     70      !!              required to create the arrays which define the space 
    6671      !!              and time domain of the ocean model. 
    6772      !! 
    … …  
    7681      CHARACTER (len=*), INTENT(in) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables 
    7782      ! 
    78       INTEGER ::   ji, jj, jk, ik   ! dummy loop indices 
     83      INTEGER ::   ji, jj, jk, jt   ! dummy loop indices 
    7984      INTEGER ::   iconf = 0    ! local integers 
    80       CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"  
     85      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))" 
    8186      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
    8287      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0 
    … …  
    110115         CASE( 7 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north-south)' 
    111116         CASE DEFAULT 
    112             CALL ctl_stop( 'jperio is out of range' ) 
     117            CALL ctl_stop( 'dom_init:   jperio is out of range' ) 
    113118         END SELECT 
    114119         WRITE(numout,*)     '      Ocean model configuration used:' 
    … …  
    140145      IF( ln_closea ) CALL dom_clo      ! Read in masks to define closed seas and lakes 
    141146 
    142       CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry 
     147      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry (return top and bottom ocean t-level indices) 
    143148 
    144149      CALL dom_msk( ik_top, ik_bot )    ! Masks 
    … …  
    147152      hu_0(:,:) = 0._wp 
    148153      hv_0(:,:) = 0._wp 
     154      hf_0(:,:) = 0._wp 
    149155      DO jk = 1, jpk 
    150156         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
    151157         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 
    152158         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 
     159         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * fmask(:,:,jk) 
    153160      END DO 
    154161      ! 
     162      r1_ht_0(:,:) = ssmask (:,:) / ( ht_0(:,:) + 1._wp -  ssmask (:,:) ) 
     163      r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp -  ssumask(:,:) ) 
     164      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp -  ssvmask(:,:) ) 
     165      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp -  ssfmask(:,:) ) 
     166 
     167      ! 
     168#if defined key_qco 
     169      !           !==  initialisation of time varying coordinate  ==!   Quasi-Euerian coordinate case 
     170      ! 
     171      IF( .NOT.l_offline )   CALL dom_qco_init( Kbb, Kmm, Kaa ) 
     172      ! 
     173      IF( ln_linssh )        CALL ctl_stop('STOP','domain: key_qco and ln_linssh = T are incompatible') 
     174      ! 
     175#else 
    155176      !           !==  time varying part of coordinate system  ==! 
    156177      ! 
    157178      IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all 
    158       ! 
    159          !       before        !          now          !       after         ! 
    160             gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   ;   gdept(:,:,:,Kaa) = gdept_0   ! depth of grid-points 
    161             gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   ;   gdepw(:,:,:,Kaa) = gdepw_0   ! 
    162                                    gde3w = gde3w_0   !        ---          ! 
    163          !                                                                   
    164               e3t(:,:,:,Kbb) =   e3t_0  ;     e3t(:,:,:,Kmm) =   e3t_0   ;   e3t(:,:,:,Kaa) =  e3t_0    ! scale factors 
    165               e3u(:,:,:,Kbb) =   e3u_0  ;     e3u(:,:,:,Kmm) =   e3u_0   ;   e3u(:,:,:,Kaa) =  e3u_0    ! 
    166               e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    ! 
    167                                      e3f =   e3f_0   !        ---          ! 
    168               e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   ;    e3w(:,:,:,Kaa) =   e3w_0   !  
    169              e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   ;   e3uw(:,:,:,Kaa) =  e3uw_0   !   
    170              e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   ;   e3vw(:,:,:,Kaa) =  e3vw_0   ! 
    171          ! 
    172          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF 
    173          z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 
    174          ! 
    175          !        before       !          now          !       after         ! 
    176                                       ht =    ht_0   !                     ! water column thickness 
    177                hu(:,:,Kbb) =    hu_0  ;      hu(:,:,Kmm) =    hu_0   ;    hu(:,:,Kaa) =    hu_0   !  
    178                hv(:,:,Kbb) =    hv_0  ;      hv(:,:,Kmm) =    hv_0   ;    hv(:,:,Kaa) =    hv_0   ! 
    179             r1_hu(:,:,Kbb) = z1_hu_0  ;   r1_hu(:,:,Kmm) = z1_hu_0   ; r1_hu(:,:,Kaa) = z1_hu_0   ! inverse of water column thickness 
    180             r1_hv(:,:,Kbb) = z1_hv_0  ;   r1_hv(:,:,Kmm) = z1_hv_0   ; r1_hv(:,:,Kaa) = z1_hv_0   ! 
    181          ! 
     179         ! 
     180         DO jt = 1, jpt                         ! depth of t- and w-grid-points 
     181            gdept(:,:,:,jt) = gdept_0(:,:,:) 
     182            gdepw(:,:,:,jt) = gdepw_0(:,:,:) 
     183         END DO 
     184            gde3w(:,:,:)    = gde3w_0(:,:,:)    ! = gdept as the sum of e3t 
     185         ! 
     186         DO jt = 1, jpt                         ! vertical scale factors 
     187            e3t(:,:,:,jt) =  e3t_0(:,:,:) 
     188            e3u(:,:,:,jt) =  e3u_0(:,:,:) 
     189            e3v(:,:,:,jt) =  e3v_0(:,:,:) 
     190            e3w(:,:,:,jt) =  e3w_0(:,:,:) 
     191            e3uw(:,:,:,jt) = e3uw_0(:,:,:) 
     192            e3vw(:,:,:,jt) = e3vw_0(:,:,:) 
     193         END DO 
     194            e3f(:,:,:)    =  e3f_0(:,:,:) 
     195         ! 
     196         DO jt = 1, jpt                         ! water column thickness and its inverse 
     197            hu(:,:,jt)    =    hu_0(:,:) 
     198            hv(:,:,jt)    =    hv_0(:,:) 
     199            r1_hu(:,:,jt) = r1_hu_0(:,:) 
     200            r1_hv(:,:,jt) = r1_hv_0(:,:) 
     201         END DO 
     202            ht(:,:) =    ht_0(:,:) 
    182203         ! 
    183204      ELSE                       != time varying : initialize before/now/after variables 
    184205         ! 
    185          IF( .NOT.l_offline )  CALL dom_vvl_init( Kbb, Kmm, Kaa ) 
    186          ! 
    187       ENDIF 
     206         IF( .NOT.l_offline )   CALL dom_vvl_init( Kbb, Kmm, Kaa ) 
     207         ! 
     208      ENDIF 
     209#endif 
     210 
    188211      ! 
    189212      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point 
    … …  
    198221         WRITE(numout,*) 'dom_init :   ==>>>   END of domain initialization' 
    199222         WRITE(numout,*) '~~~~~~~~' 
    200          WRITE(numout,*)  
     223         WRITE(numout,*) 
    201224      ENDIF 
    202225      ! 
    … …  
    210233      !! ** Purpose :   initialization of global domain <--> local domain indices 
    211234      !! 
    212       !! ** Method  :    
     235      !! ** Method  : 
    213236      !! 
    214237      !! ** Action  : - mig , mjg : local  domain indices ==> global domain indices 
    … …  
    226249      END DO 
    227250      !                              ! global domain indices ==> local domain indices 
    228       !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the  
    229       !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.  
     251      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the 
     252      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain. 
    230253      DO ji = 1, jpiglo 
    231254        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) ) 
    … …  
    273296      !!---------------------------------------------------------------------- 
    274297      !!                     ***  ROUTINE dom_nam  *** 
    275       !!                     
     298      !! 
    276299      !! ** Purpose :   read domaine namelists and print the variables. 
    277300      !! 
    … …  
    355378      l_1st_euler = ln_1st_euler 
    356379      IF( .NOT. l_1st_euler .AND. .NOT. ln_rstart ) THEN 
    357          IF(lwp) WRITE(numout,*)   
     380         IF(lwp) WRITE(numout,*) 
    358381         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)' 
    359382         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : l_1st_euler is forced to .true. '    
    … …  
    383406      IF(lwp) WRITE(numout,*) 
    384407      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL 
    385       CASE (  1 )  
     408      CASE (  1 ) 
    386409         CALL ioconf_calendar('gregorian') 
    387410         IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "gregorian", i.e. leap year' 
    … …  
    419442      IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    420443         lrxios = ln_xios_read.AND.ln_rstart 
    421 !set output file type for XIOS based on NEMO namelist  
    422          IF (nn_wxios > 0) lwxios = .TRUE.  
     444!set output file type for XIOS based on NEMO namelist 
     445         IF (nn_wxios > 0) lwxios = .TRUE. 
    423446         nxioso = nn_wxios 
    424447      ENDIF 
    … …  
    463486      !!---------------------------------------------------------------------- 
    464487      INTEGER, DIMENSION(2) ::   imi1, imi2, ima1, ima2 
    465       INTEGER, DIMENSION(2) ::   iloc   !  
     488      INTEGER, DIMENSION(2) ::   iloc   ! 
    466489      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max 
    467490      !!---------------------------------------------------------------------- 
    … …  
    473496         CALL mpp_maxloc( 'domain', e2t(:,:), tmask_i(:,:), ze2max, ima2 ) 
    474497      ELSE 
    475          ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
    476          ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )     
    477          ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
    478          ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )     
     498         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
     499         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 
     500         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
     501         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    479502         ! 
    480503         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    … …  
    507530      !!---------------------------------------------------------------------- 
    508531      !!                     ***  ROUTINE dom_nam  *** 
    509       !!                     
     532      !! 
    510533      !! ** Purpose :   read the domain size in domain configuration file 
    511534      !! 
    … …  
    514537      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name 
    515538      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution 
    516       INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes  
    517       INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c.  
     539      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
     540      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
    518541      ! 
    519542      INTEGER ::   inum   ! local integer 
    … …  
    547570         cd_cfg = 'UNKNOWN' 
    548571         kk_cfg = -9999999 
    549                                           !- or they may be present as global attributes  
    550                                           !- (netcdf only)   
     572                                          !- or they may be present as global attributes 
     573                                          !- (netcdf only) 
    551574         CALL iom_getatt( inum, 'cn_cfg', cd_cfg )  ! returns   !  if not found 
    552575         CALL iom_getatt( inum, 'nn_cfg', kk_cfg )  ! returns -999 if not found 
    … …  
    570593         WRITE(numout,*) '      type of global domain lateral boundary   jperio = ', kperio 
    571594      ENDIF 
    572       !         
     595      ! 
    573596   END SUBROUTINE domain_cfg 
    574     
    575     
     597 
     598 
    576599   SUBROUTINE cfg_write 
    577600      !!---------------------------------------------------------------------- 
    578601      !!                  ***  ROUTINE cfg_write  *** 
    579       !!                    
    580       !! ** Purpose :   Create the "cn_domcfg_out" file, a NetCDF file which  
    581       !!              contains all the ocean domain informations required to  
     602      !! 
     603      !! ** Purpose :   Create the "cn_domcfg_out" file, a NetCDF file which 
     604      !!              contains all the ocean domain informations required to 
    582605      !!              define an ocean configuration. 
    583606      !! 
    … …  
    585608      !!              ocean configuration. 
    586609      !! 
    587       !! ** output file :   domcfg_out.nc : domain size, characteristics, horizontal  
     610      !! ** output file :   domcfg_out.nc : domain size, characteristics, horizontal 
    588611      !!                       mesh, Coriolis parameter, and vertical scale factors 
    589612      !!                    NB: also contain ORCA family information 
    … …  
    603626      !                       !  create 'domcfg_out.nc' file  ! 
    604627      !                       ! ============================= ! 
    605       !          
     628      ! 
    606629      clnam = cn_domcfg_out  ! filename (configuration information) 
    607630      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE. ) 
    608        
     631 
    609632      ! 
    610633      !                             !==  ORCA family specificities  ==! 
    611634      IF( cn_cfg == "ORCA" ) THEN 
    612635         CALL iom_rstput( 0, 0, inum, 'ORCA'      , 1._wp            , ktype = jp_i4 ) 
    613          CALL iom_rstput( 0, 0, inum, 'ORCA_index', REAL( nn_cfg, wp), ktype = jp_i4 )          
     636         CALL iom_rstput( 0, 0, inum, 'ORCA_index', REAL( nn_cfg, wp), ktype = jp_i4 ) 
    614637      ENDIF 
    615638      ! 
    … …  
    643666      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 ) 
    644667      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 ) 
    645       !                                 
     668      ! 
    646669      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )   ! longitude 
    647670      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 ) 
    648671      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 ) 
    649672      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 ) 
    650       !                                 
     673      ! 
    651674      CALL iom_rstput( 0, 0, inum, 'e1t'  , e1t  , ktype = jp_r8 )   ! i-scale factors (e1.) 
    652675      CALL iom_rstput( 0, 0, inum, 'e1u'  , e1u  , ktype = jp_r8 ) 
    … …  
    663686      ! 
    664687      !                             !==  vertical mesh  ==! 
    665       !                                                      
     688      ! 
    666689      CALL iom_rstput( 0, 0, inum, 'e3t_1d'  , e3t_1d , ktype = jp_r8 )   ! reference 1D-coordinate 
    667690      CALL iom_rstput( 0, 0, inum, 'e3w_1d'  , e3w_1d , ktype = jp_r8 ) 
    … …  
    674697      CALL iom_rstput( 0, 0, inum, 'e3uw_0'  , e3uw_0 , ktype = jp_r8 ) 
    675698      CALL iom_rstput( 0, 0, inum, 'e3vw_0'  , e3vw_0 , ktype = jp_r8 ) 
    676       !                                          
     699      ! 
    677700      !                             !==  wet top and bottom level  ==!   (caution: multiplied by ssmask) 
    678701      ! 
    … …  
    694717      ! 
    695718      !                                ! ============================ 
    696       !                                !        close the files  
     719      !                                !        close the files 
    697720      !                                ! ============================ 
    698721      CALL iom_close( inum ) 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/dommsk.F90

    r12377 r13151  
    22   !!====================================================================== 
    33   !!                       ***  MODULE dommsk   *** 
    4    !! Ocean initialization : domain land/sea mask  
     4   !! Ocean initialization : domain land/sea mask 
    55   !!====================================================================== 
    66   !! History :  OPA  ! 1987-07  (G. Madec)  Original code 
    … …  
    1818   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask 
    1919   !!            4.0  ! 2016-06  (G. Madec, S. Flavoni)  domain configuration / user defined interface 
     20   !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
    2021   !!---------------------------------------------------------------------- 
    2122 
    … …  
    4041   !                            !!* Namelist namlbc : lateral boundary condition * 
    4142   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity 
    42    LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition  
     43   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition 
    4344   !                                            with analytical eqs. 
    4445 
    … …  
    4748   !!---------------------------------------------------------------------- 
    4849   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    49    !! $Id$  
     50   !! $Id$ 
    5051   !! Software governed by the CeCILL license (see ./LICENSE) 
    5152   !!---------------------------------------------------------------------- 
    … …  
    5960      !!      zontal velocity points (u & v), vorticity points (f) points. 
    6061      !! 
    61       !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top  
    62       !!      and ko_bot, the indices of the fist and last ocean t-levels which  
     62      !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top 
     63      !!      and ko_bot, the indices of the fist and last ocean t-levels which 
    6364      !!      are either defined in usrdef_zgr or read in zgr_read. 
    64       !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)  
     65      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask) 
    6566      !!      are deduced from a product of the two neighboring tmask. 
    6667      !!                The vorticity mask (fmask) is deduced from tmask taking 
    … …  
    7778      !!                due to cyclic or North Fold boundaries as well as MPP halos. 
    7879      !! 
    79       !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask  
     80      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask 
    8081      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.) 
    81       !!               fmask   : land/ocean mask at f-point (=0., or =1., or  
     82      !!               fmask   : land/ocean mask at f-point (=0., or =1., or 
    8283      !!                         =rn_shlat along lateral boundaries) 
    83       !!               tmask_i : interior ocean mask  
     84      !!               tmask_i : interior ocean mask 
    8485      !!               tmask_h : halo mask 
    8586      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask 
    … …  
    108109902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlbc in configuration namelist' ) 
    109110      IF(lwm) WRITE ( numond, namlbc ) 
    110        
     111 
    111112      IF(lwp) THEN                  ! control print 
    112113         WRITE(numout,*) 
    … …  
    115116         WRITE(numout,*) '   Namelist namlbc' 
    116117         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat 
    117          WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat  
     118         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
    118119      ENDIF 
    119120      ! 
    … …  
    140141      ! 
    141142      ! the following call is mandatory 
    142       ! it masks boundaries (bathy=0) where needed depending on the configuration (closed, periodic...)   
     143      ! it masks boundaries (bathy=0) where needed depending on the configuration (closed, periodic...) 
    143144      CALL lbc_lnk( 'dommsk', tmask  , 'T', 1._wp )      ! Lateral boundary conditions 
    144  
     145       
    145146     ! Mask corrections for bdy (read in mppini2) 
    146147      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
    … …  
    157158         END_3D 
    158159      ENDIF 
    159           
     160 
    160161      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask) 
    161162      ! ---------------------------------------- 
    … …  
    174175      END DO 
    175176      CALL lbc_lnk_multi( 'dommsk', umask, 'U', 1., vmask, 'V', 1., fmask, 'F', 1. )      ! Lateral boundary conditions 
    176   
     177 
    177178      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask) 
    178179      !----------------------------------------- 
    … …  
    182183      DO jk = 2, jpk                   ! interior values 
    183184         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 
    184          wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)    
     185         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1) 
    185186         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 
    186187      END DO 
    … …  
    192193      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 ) 
    193194      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 
     195      ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 ) 
    194196 
    195197 
    … …  
    201203      ! 
    202204      !                          ! halo mask : 0 on the halo and 1 elsewhere 
    203       tmask_h(:,:) = 1._wp                   
     205      tmask_h(:,:) = 1._wp 
    204206      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns 
    205207      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
    … …  
    208210      ! 
    209211      !                          ! north fold mask 
    210       tpol(1:jpiglo) = 1._wp  
     212      tpol(1:jpiglo) = 1._wp 
    211213      fpol(1:jpiglo) = 1._wp 
    212214      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot 
    … …  
    225227      ENDIF 
    226228      ! 
    227       !                          ! interior mask : 2D ocean mask x halo mask  
     229      !                          ! interior mask : 2D ocean mask x halo mask 
    228230      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:) 
    229231 
    230232 
    231233      ! Lateral boundary conditions on velocity (modify fmask) 
    232       ! ---------------------------------------   
     234      ! --------------------------------------- 
    233235      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition 
    234236         ! 
    … …  
    236238         ! 
    237239         DO jk = 1, jpk 
    238             zwf(:,:) = fmask(:,:,jk)          
     240            zwf(:,:) = fmask(:,:,jk) 
    239241            DO_2D_00_00 
    240242               IF( fmask(ji,jj,jk) == 0._wp ) THEN 
    … …  
    250252                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
    251253               ENDIF 
    252             END DO          
     254            END DO 
    253255            DO ji = 2, jpim1 
    254256               IF( fmask(ji,1,jk) == 0._wp ) THEN 
    … …  
    259261               ENDIF 
    260262            END DO 
    261 #if defined key_agrif  
    262             IF( .NOT. AGRIF_Root() ) THEN  
    263                IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east  
    264                IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west  
    265                IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north  
    266                IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south  
    267             ENDIF  
    268 #endif  
     263#if defined key_agrif 
     264            IF( .NOT. AGRIF_Root() ) THEN 
     265               IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east 
     266               IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west 
     267               IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north 
     268               IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south 
     269            ENDIF 
     270#endif 
    269271         END DO 
    270272         ! 
    … …  
    276278         ! 
    277279      ENDIF 
    278        
     280 
    279281      ! User defined alteration of fmask (use to reduce ocean transport in specified straits) 
    280       ! --------------------------------  
     282      ! -------------------------------- 
    281283      ! 
    282284      CALL usr_def_fmask( cn_cfg, nn_cfg, fmask ) 
    283285      ! 
    284286   END SUBROUTINE dom_msk 
    285     
     287 
    286288   !!====================================================================== 
    287289END MODULE dommsk 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/domvvl.F90

    r12489 r13151  
    22   !!====================================================================== 
    33   !!                       ***  MODULE domvvl   *** 
    4    !! Ocean :  
     4   !! Ocean : 
    55   !!====================================================================== 
    66   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
    … …  
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    1010   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 
     11   !!            4.x  ! 2020-02  (G. Madec, S. Techene) introduce ssh to h0 ratio 
    1112   !!---------------------------------------------------------------------- 
    1213 
    13    !!---------------------------------------------------------------------- 
    14    !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
    15    !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
    16    !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid 
    17    !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
    18    !!   dom_vvl_rst      : read/write restart file 
    19    !!   dom_vvl_ctl      : Check the vvl options 
    20    !!---------------------------------------------------------------------- 
    2114   USE oce             ! ocean dynamics and tracers 
    2215   USE phycst          ! physical constant 
    … …  
    3528   IMPLICIT NONE 
    3629   PRIVATE 
    37  
    38    PUBLIC  dom_vvl_init       ! called by domain.F90 
    39    PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90 
    40    PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    41    PUBLIC  dom_vvl_sf_update  ! called by step.F90 
    42    PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    43  
     30    
    4431   !                                                      !!* Namelist nam_vvl 
    4532   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate 
    … …  
    6350   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6451 
     52#if defined key_qco 
     53   !!---------------------------------------------------------------------- 
     54   !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate 
     55   !!---------------------------------------------------------------------- 
     56#else 
     57   !!---------------------------------------------------------------------- 
     58   !!   Default key      Old management of time varying vertical coordinate 
     59   !!---------------------------------------------------------------------- 
     60    
     61   !!---------------------------------------------------------------------- 
     62   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
     63   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
     64   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid 
     65   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
     66   !!   dom_vvl_rst      : read/write restart file 
     67   !!   dom_vvl_ctl      : Check the vvl options 
     68   !!---------------------------------------------------------------------- 
     69 
     70   PUBLIC  dom_vvl_init       ! called by domain.F90 
     71   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90 
     72   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
     73   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
     74   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
     75    
    6576   !! * Substitutions 
    6677#  include "do_loop_substitute.h90" 
    … …  
    98109      !!---------------------------------------------------------------------- 
    99110      !!                ***  ROUTINE dom_vvl_init  *** 
    100       !!                    
     111      !! 
    101112      !! ** Purpose :  Initialization of all scale factors, depths 
    102113      !!               and water column heights 
    … …  
    107118      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
    108119      !!              - Regrid: e3[u/v](:,:,:,Kmm) 
    109       !!                        e3[u/v](:,:,:,Kmm)        
    110       !!                        e3w(:,:,:,Kmm)            
     120      !!                        e3[u/v](:,:,:,Kmm) 
     121      !!                        e3w(:,:,:,Kmm) 
    111122      !!                        e3[u/v]w_b 
    112       !!                        e3[u/v]w_n       
     123      !!                        e3[u/v]w_n 
    113124      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 
    114125      !!              - h(t/u/v)_0 
    … …  
    135146      ! 
    136147   END SUBROUTINE dom_vvl_init 
    137    ! 
     148 
     149 
    138150   SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 
    139151      !!---------------------------------------------------------------------- 
    140152      !!                ***  ROUTINE dom_vvl_init  *** 
    141       !!                    
    142       !! ** Purpose :  Interpolation of all scale factors,  
     153      !! 
     154      !! ** Purpose :  Interpolation of all scale factors, 
    143155      !!               depths and water column heights 
    144156      !! 
    … …  
    147159      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
    148160      !!              - Regrid: e3(u/v)_n 
    149       !!                        e3(u/v)_b        
    150       !!                        e3w_n            
    151       !!                        e3(u/v)w_b       
    152       !!                        e3(u/v)w_n       
     161      !!                        e3(u/v)_b 
     162      !!                        e3w_n 
     163      !!                        e3(u/v)w_b 
     164      !!                        e3(u/v)w_n 
    153165      !!                        gdept_n, gdepw_n and gde3w_n 
    154166      !!              - h(t/u/v)_0 
    … …  
    168180      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
    169181      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
    170       CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
     182      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V 
    171183      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    172184      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
    173       !                                ! Vertical interpolation of e3t,u,v  
     185      !                                ! Vertical interpolation of e3t,u,v 
    174186      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
    175187      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
    … …  
    193205         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    194206         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    195          !                             ! 0.5 where jk = mikt      
     207         !                             ! 0.5 where jk = mikt 
    196208!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    197209         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    198210         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    199211         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    200             &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     212            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
    201213         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    202214         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    203215         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    204             &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     216            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
    205217      END_3D 
    206218      ! 
    … …  
    261273            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    262274               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
    263                   ii0 = 103   ;   ii1 = 111        
    264                   ij0 = 128   ;   ij1 = 135   ;    
     275                  ii0 = 103   ;   ii1 = 111 
     276                  ij0 = 128   ;   ij1 = 135   ; 
    265277                  frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp 
    266278                  frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rn_Dt 
    … …  
    280292            CALL iom_set_rstw_var_active('tilde_e3t_n') 
    281293         END IF 
    282          !                                           ! -------------!     
     294         !                                           ! -------------! 
    283295         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    284296            !                                        ! ------------ ! 
    … …  
    291303 
    292304 
    293    SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )  
     305   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
    294306      !!---------------------------------------------------------------------- 
    295307      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
    296       !!                    
     308      !! 
    297309      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt, 
    298310      !!                 tranxt and dynspg routines 
    299311      !! 
    300312      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness. 
    301       !!               - z_tilde_case: after scale factor increment =  
     313      !!               - z_tilde_case: after scale factor increment = 
    302314      !!                                    high frequency part of horizontal divergence 
    303315      !!                                  + retsoring towards the background grid 
    … …  
    307319      !! 
    308320      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case 
    309       !!               - tilde_e3t_a: after increment of vertical scale factor  
     321      !!               - tilde_e3t_a: after increment of vertical scale factor 
    310322      !!                              in z_tilde case 
    311323      !!               - e3(t/u/v)_a 
    … …  
    410422            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    411423               &            * ( 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)           &  
     424            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
    413425               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    414426            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    … …  
    460472               WRITE(numout, *) 'at i, j, k=', ijk_max 
    461473               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
    462                WRITE(numout, *) 'at i, j, k=', ijk_min             
     474               WRITE(numout, *) 'at i, j, k=', ijk_min 
    463475               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high') 
    464476            ENDIF 
    … …  
    575587      !!---------------------------------------------------------------------- 
    576588      !!                ***  ROUTINE dom_vvl_sf_update  *** 
    577       !!                    
    578       !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors  
     589      !! 
     590      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors 
    579591      !!               compute all depths and related variables for next time step 
    580592      !!               write outputs and restart file 
    … …  
    586598      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step 
    587599      !!               - Recompute: 
    588       !!                    e3(u/v)_b        
    589       !!                    e3w(:,:,:,Kmm)            
    590       !!                    e3(u/v)w_b       
    591       !!                    e3(u/v)w_n       
     600      !!                    e3(u/v)_b 
     601      !!                    e3w(:,:,:,Kmm) 
     602      !!                    e3(u/v)w_b 
     603      !!                    e3(u/v)w_n 
    592604      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w 
    593605      !!                    h(u/v) and h(u/v)r 
    … …  
    620632            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 
    621633         ELSE 
    622             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &  
     634            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
    623635            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 
    624636         ENDIF 
    … …  
    632644      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 
    633645      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    634        
     646 
    635647      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  ) 
    636        
     648 
    637649      ! Vertical scale factor interpolations 
    638650      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
    … …  
    653665         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    654666         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    655              &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     667             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
    656668         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    657669      END_3D 
    … …  
    772784      !!--------------------------------------------------------------------- 
    773785      !!                   ***  ROUTINE dom_vvl_rst  *** 
    774       !!                      
     786      !! 
    775787      !! ** Purpose :   Read or write VVL file in restart file 
    776788      !! 
    … …  
    789801      !!---------------------------------------------------------------------- 
    790802      ! 
    791       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     803      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
    792804         !                                   ! =============== 
    793805         IF( ln_rstart ) THEN                   !* Read the restart file 
    … …  
    808820               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    809821               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    810                ! needed to restart if land processor not computed  
     822               ! needed to restart if land processor not computed 
    811823               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 
    812                WHERE ( tmask(:,:,:) == 0.0_wp )  
     824               WHERE ( tmask(:,:,:) == 0.0_wp ) 
    813825                  e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
    814826                  e3t(:,:,:,Kbb) = e3t_0(:,:,:) 
    … …  
    873885            ! 
    874886 
    875             IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential  
     887            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential 
    876888               ! 
    877889               IF( cn_cfg == 'wad' ) THEN 
    … …  
    908920                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    909921                     ENDIF 
    910                   END DO  
    911                END DO  
     922                  END DO 
     923               END DO 
    912924               ! 
    913925            ELSE 
    … …  
    950962            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios) 
    951963         END IF 
    952          !                                           ! -------------!     
     964         !                                           ! -------------! 
    953965         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    954966            !                                        ! ------------ ! 
    … …  
    965977      !!--------------------------------------------------------------------- 
    966978      !!                  ***  ROUTINE dom_vvl_ctl  *** 
    967       !!                 
     979      !! 
    968980      !! ** Purpose :   Control the consistency between namelist options 
    969981      !!                for vertical coordinate 
    … …  
    974986         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , & 
    975987         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe 
    976       !!----------------------------------------------------------------------  
     988      !!---------------------------------------------------------------------- 
    977989      ! 
    978990      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 
    … …  
    10311043   END SUBROUTINE dom_vvl_ctl 
    10321044 
     1045#endif 
     1046 
    10331047   !!====================================================================== 
    10341048END MODULE domvvl 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/istate.F90

    r12489 r13151  
    4343   !! * Substitutions 
    4444#  include "do_loop_substitute.h90" 
     45#  include "domzgr_substitute.h90" 
    4546   !!---------------------------------------------------------------------- 
    4647   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5960      ! 
    6061      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     62      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute 
    6163!!gm see comment further down 
    6264      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace 
    … …  
    115117            ! 
    116118         ELSE                                 ! user defined initial T and S 
    117             CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
     119            DO jk = 1, jpk 
     120               zgdept(:,:,jk) = gdept(:,:,jk,Kbb) 
     121            END DO 
     122            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
    118123         ENDIF 
    119124         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    … …  
    127132!!gm POTENTIAL BUG : 
    128133!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
    129 !!             as well as gdept and gdepw....   !!!!!  
     134!!             as well as gdept_ and gdepw_....   !!!!!  
    130135!!      ===>>>>   probably a call to domvvl initialisation here.... 
    131136 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/divhor.F90

    r12377 r13151  
    2121   USE dom_oce         ! ocean space and time domain 
    2222   USE sbc_oce, ONLY : ln_rnf      ! river runoff 
    23    USE sbcrnf , ONLY : sbc_rnf_div ! river runoff  
     23   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff 
    2424   USE isf_oce, ONLY : ln_isf      ! ice shelf 
    2525   USE isfhdiv, ONLY : isf_hdiv    ! ice shelf 
    26 #if defined key_asminc    
     26#if defined key_asminc 
    2727   USE asminc          ! Assimilation increment 
    2828#endif 
    … …  
    4040   !! * Substitutions 
    4141#  include "do_loop_substitute.h90" 
     42#  include "domzgr_substitute.h90" 
    4243   !!---------------------------------------------------------------------- 
    4344   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    44    !! $Id$  
     45   !! $Id$ 
    4546   !! Software governed by the CeCILL license (see ./LICENSE) 
    4647   !!---------------------------------------------------------------------- 
    … …  
    5051      !!---------------------------------------------------------------------- 
    5152      !!                  ***  ROUTINE div_hor  *** 
    52       !!                     
     53      !! 
    5354      !! ** Purpose :   compute the horizontal divergence at now time-step 
    5455      !! 
    5556      !! ** Method  :   the now divergence is computed as : 
    5657      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    57       !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)  
     58      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 
    5859      !! 
    5960      !! ** Action  : - update hdiv, the now horizontal divergence 
    … …  
    7879      DO_3D_00_00( 1, jpkm1 ) 
    7980         hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      & 
    80             &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      & 
    81             &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      & 
    82             &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   & 
     81            &              - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      & 
     82            &              + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      & 
     83            &              - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   & 
    8384            &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    8485      END_3D 
    … …  
    9596      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                     !==  runoffs    ==!   (update hdiv field) 
    9697      ! 
    97 #if defined key_asminc  
     98#if defined key_asminc 
    9899      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )   !==  SSH assimilation  ==!   (update hdiv field) 
    99       !  
     100      ! 
    100101#endif 
    101102      ! 
    … …  
    107108      ! 
    108109   END SUBROUTINE div_hor 
    109     
     110 
    110111   !!====================================================================== 
    111112END MODULE divhor 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynadv_cen2.F90

    r12377 r13151  
    2828   !! * Substitutions 
    2929#  include "do_loop_substitute.h90" 
     30#  include "domzgr_substitute.h90" 
    3031   !!---------------------------------------------------------------------- 
    3132   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    7980         DO_2D_00_00 
    8081            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
    81                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     82               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   & 
     83               &                           / e3u(ji,jj,jk,Kmm) 
    8284            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
    83                &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     85               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   & 
     86               &                           / e3v(ji,jj,jk,Kmm) 
    8487         END_2D 
    8588      END DO 
    … …  
    115118      END DO 
    116119      DO_3D_00_00( 1, jpkm1 ) 
    117          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
    118          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     120         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   & 
     121            &                                      / e3u(ji,jj,jk,Kmm) 
     122         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj)   & 
     123            &                                      / e3v(ji,jj,jk,Kmm) 
    119124      END_3D 
    120125      ! 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynadv_ubs.F90

    r12377 r13151  
    3434   !! * Substitutions 
    3535#  include "do_loop_substitute.h90" 
     36#  include "domzgr_substitute.h90" 
    3637   !!---------------------------------------------------------------------- 
    3738   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    169170         DO_2D_00_00 
    170171            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    & 
    171                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     172               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   & 
     173               &                           / e3u(ji,jj,jk,Kmm) 
    172174            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    & 
    173                &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     175               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   & 
     176               &                           / e3v(ji,jj,jk,Kmm) 
    174177         END_2D 
    175178      END DO 
    … …  
    206209      END DO 
    207210      DO_3D_00_00( 1, jpkm1 ) 
    208          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
    209          pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     211         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   & 
     212            &                                       / e3u(ji,jj,jk,Kmm) 
     213         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj)   & 
     214            &                                       / e3v(ji,jj,jk,Kmm) 
    210215      END_3D 
    211216      ! 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynatf.F90

    r12489 r13151  
    1313   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond. 
    1414   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
    15    !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.  
     15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines. 
    1616   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option 
    1717   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module 
    … …  
    2222   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering. 
    2323   !!------------------------------------------------------------------------- 
    24    
     24 
    2525   !!---------------------------------------------------------------------------------------------- 
    2626   !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors 
    … …  
    4242   USE trdken         ! trend manager: kinetic energy 
    4343   USE isf_oce   , ONLY: ln_isf     ! ice shelf 
    44    USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine  
     44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 
    4545   ! 
    4646   USE in_out_manager ! I/O manager 
    … …  
    5959   PUBLIC    dyn_atf   ! routine called by step.F90 
    6060 
     61#if defined key_qco 
     62   !!---------------------------------------------------------------------- 
     63   !!   'key_qco'      EMPTY ROUTINE     Quasi-Eulerian vertical coordonate 
     64   !!---------------------------------------------------------------------- 
     65CONTAINS 
     66 
     67   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 
     68      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index 
     69      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices 
     70      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered 
     71      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered 
     72 
     73      WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt 
     74   END SUBROUTINE dyn_atf 
     75 
     76#else 
     77 
    6178   !! * Substitutions 
    6279#  include "do_loop_substitute.h90" 
    6380   !!---------------------------------------------------------------------- 
    6481   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    65    !! $Id$  
     82   !! $Id$ 
    6683   !! Software governed by the CeCILL license (see ./LICENSE) 
    6784   !!---------------------------------------------------------------------- 
    … …  
    7188      !!---------------------------------------------------------------------- 
    7289      !!                  ***  ROUTINE dyn_atf  *** 
    73       !!                    
    74       !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary  
     90      !! 
     91      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary 
    7592      !!             condition on the after velocity and apply the Asselin time 
    7693      !!             filter to the now fields. 
    … …  
    7996      !!             estimate (ln_dynspg_ts=T) 
    8097      !! 
    81       !!              * Apply lateral boundary conditions on after velocity  
     98      !!              * Apply lateral boundary conditions on after velocity 
    8299      !!             at the local domain boundaries through lbc_lnk call, 
    83100      !!             at the one-way open boundaries (ln_bdy=T), 
    … …  
    86103      !!              * Apply the Asselin time filter to the now fields 
    87104      !!             arrays to start the next time step: 
    88       !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))  
     105      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 
    89106      !!                                    + rn_atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 
    90107      !!             Note that with flux form advection and non linear free surface, 
    … …  
    92109      !!             As a result, dyn_atf MUST be called after tra_atf. 
    93110      !! 
    94       !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity  
     111      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity 
    95112      !!---------------------------------------------------------------------- 
    96113      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index 
    … …  
    103120      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      - 
    104121      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld 
    105       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva  
     122      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
    106123      !!---------------------------------------------------------------------- 
    107124      ! 
    … …  
    131148         ! 
    132149         IF( .NOT.ln_bt_fw ) THEN 
    133             ! Remove advective velocity from "now velocities"  
    134             ! prior to asselin filtering      
    135             ! In the forward case, this is done below after asselin filtering    
    136             ! so that asselin contribution is removed at the same time  
     150            ! Remove advective velocity from "now velocities" 
     151            ! prior to asselin filtering 
     152            ! In the forward case, this is done below after asselin filtering 
     153            ! so that asselin contribution is removed at the same time 
    137154            DO jk = 1, jpkm1 
    138155               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
    139156               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    140             END DO   
     157            END DO 
    141158         ENDIF 
    142159      ENDIF 
    143160 
    144161      ! Update after velocity on domain lateral boundaries 
    145       ! --------------------------------------------------       
     162      ! -------------------------------------------------- 
    146163# if defined key_agrif 
    147164      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries 
    … …  
    198215            zwfld(:,:) = emp_b(:,:) - emp(:,:) 
    199216            IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) ) 
     217 
    200218            DO jk = 1, jpkm1 
    201219               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 
    202                               &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) )  
     220                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 
    203221            END DO 
    204222            ! 
    … …  
    237255                  pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    238256               END_3D 
    239                pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1)   
     257               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 
    240258               pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 
    241259               ! 
    … …  
    248266         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    249267            ! Revert filtered "now" velocities to time split estimate 
    250             ! Doing it here also means that asselin filter contribution is removed   
     268            ! Doing it here also means that asselin filter contribution is removed 
    251269            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
    252             zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)     
     270            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
    253271            DO jk = 2, jpkm1 
    254272               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
    255                zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)     
     273               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
    256274            END DO 
    257275            DO jk = 1, jpkm1 
    … …  
    305323      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   & 
    306324         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask ) 
    307       !  
     325      ! 
    308326      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve ) 
    309327      IF( l_trddyn     )   DEALLOCATE( zua, zva ) 
    … …  
    312330   END SUBROUTINE dyn_atf 
    313331 
     332#endif 
     333 
    314334   !!========================================================================= 
    315335END MODULE dynatf 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynhpg.F90

    r12377 r13151  
    4343   USE in_out_manager  ! I/O manager 
    4444   USE prtctl          ! Print control 
    45    USE lbclnk          ! lateral boundary condition  
     45   USE lbclnk          ! lateral boundary condition 
    4646   USE lib_mpp         ! MPP library 
    4747   USE eosbn2          ! compute density 
    … …  
    7676   !! * Substitutions 
    7777#  include "do_loop_substitute.h90" 
     78#  include "domzgr_substitute.h90" 
     79 
    7880   !!---------------------------------------------------------------------- 
    7981   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    204206      ! 
    205207      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    206       !  
     208      ! 
    207209      IF(lwp) THEN 
    208210         WRITE(numout,*) 
    … …  
    217219         WRITE(numout,*) 
    218220      ENDIF 
    219       !                           
     221      ! 
    220222   END SUBROUTINE dyn_hpg_init 
    221223 
    … …  
    427429            zcpx(ji,jj) = 0._wp 
    428430          END IF 
    429     
     431 
    430432          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    431433               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    … …  
    452454      DO_2D_00_00 
    453455         ! hydrostatic pressure gradient along s-surfaces 
    454          zhpi(ji,jj,1) = zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    & 
    455             &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj) 
    456          zhpj(ji,jj,1) = zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    & 
    457             &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj) 
     456         zhpi(ji,jj,1) =   & 
     457            &  zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    & 
     458            &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) & 
     459            &           * r1_e1u(ji,jj) 
     460         zhpj(ji,jj,1) =   & 
     461            &  zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    & 
     462            &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) & 
     463            &           * r1_e2v(ji,jj) 
    458464         ! s-coordinate pressure gradient correction 
    459465         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    … …  
    464470         IF( ln_wd_il ) THEN 
    465471            zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    466             zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     472            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
    467473            zuap = zuap * zcpx(ji,jj) 
    468474            zvap = zvap * zcpy(ji,jj) 
    … …  
    478484         ! hydrostatic pressure gradient along s-surfaces 
    479485         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    480             &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    481             &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     486            &    * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )  & 
     487            &       - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    482488         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   & 
    483             &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    484             &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     489            &    * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     490            &       - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    485491         ! s-coordinate pressure gradient correction 
    486492         zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    … …  
    491497         IF( ln_wd_il ) THEN 
    492498            zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    493             zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     499            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
    494500            zuap = zuap * zcpx(ji,jj) 
    495501            zvap = zvap * zcpy(ji,jj) 
    … …  
    522528      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
    523529      !!      iceload is added 
    524       !!       
     530      !! 
    525531      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
    526532      !!---------------------------------------------------------------------- 
    … …  
    540546      znad=1._wp                 ! To use density and not density anomaly 
    541547      ! 
    542       !                          ! iniitialised to 0. zhpi zhpi  
     548      !                          ! iniitialised to 0. zhpi zhpi 
    543549      zhpi(:,:,:) = 0._wp   ;   zhpj(:,:,:) = 0._wp 
    544550 
    … …  
    554560      CALL eos( zts_top, risfdep, zrhdtop_oce ) 
    555561 
    556 !==================================================================================      
    557 !===== Compute surface value =====================================================  
     562!================================================================================== 
     563!===== Compute surface value ===================================================== 
    558564!================================================================================== 
    559565      DO_2D_00_00 
    … …  
    567573            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
    568574            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    569             &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            )  
     575            &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            ) 
    570576         zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    & 
    571577            &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   & 
    572             &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &  
     578            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
    573579            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    574             &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            )  
     580            &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            ) 
    575581         ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    576582         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    … …  
    582588         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
    583589      END_2D 
    584 !==================================================================================      
    585 !===== Compute interior value =====================================================  
     590!================================================================================== 
     591!===== Compute interior value ===================================================== 
    586592!================================================================================== 
    587593      ! interior value (2=<jk=<jpkm1) 
    … …  
    589595         ! hydrostatic pressure gradient along s-surfaces 
    590596         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    591             &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
    592             &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
     597            &           * (  e3w(ji+1,jj,jk,Kmm)                   & 
     598            &                  * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
     599            &              - e3w(ji  ,jj,jk,Kmm)                   & 
     600            &                  * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
    593601         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    594             &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
    595             &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
     602            &           * (  e3w(ji,jj+1,jk,Kmm)                   & 
     603            &                  * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
     604            &              - e3w(ji,jj  ,jk,Kmm)                   & 
     605            &                  * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
    596606         ! s-coordinate pressure gradient correction 
    597607         zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    … …  
    650660            zcpx(ji,jj) = 0._wp 
    651661          END IF 
    652     
     662 
    653663          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    654664               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    … …  
    771781      !------------------------------------------------------------- 
    772782 
    773 !!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   to be verified 
    774 !          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
     783!!bug gm   :  e3w-gde3w(:,:,:) = 0.5*e3w  ....  and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) ....   to be verified 
     784!          true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
    775785 
    776786      DO_2D_00_00 
    … …  
    825835         IF( ln_wd_il ) THEN 
    826836           zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    827            zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     837           zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
    828838         ENDIF 
    829839         ! add to the general momentum trend 
    … …  
    845855         IF( ln_wd_il ) THEN 
    846856           zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    847            zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     857           zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
    848858         ENDIF 
    849859         ! add to the general momentum trend 
    … …  
    916926            zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
    917927                        &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) ) 
    918             
     928 
    919929             zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
    920930          ELSE 
    921931            zcpx(ji,jj) = 0._wp 
    922932          END IF 
    923     
     933 
    924934          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    925935               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    … …  
    10021012!!gm BUG ?    if it is ssh at u- & v-point then it should be: 
    10031013!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 
    1004 !                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
     1014!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
    10051015!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 
    1006 !                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
     1016!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
    10071017!!gm not this: 
    10081018       zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 
    1009                       & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
     1019                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
    10101020       zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 
    1011                       & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
     1021                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
    10121022      END_2D 
    10131023 
    … …  
    10151025 
    10161026      DO_2D_00_00 
    1017        zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)  
     1027       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 
    10181028       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 
    10191029      END_2D 
    … …  
    10981108            zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 
    10991109         ENDIF 
    1100          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
     1110         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
    11011111      ENDIF 
    11021112 
    … …  
    11541164         ENDIF 
    11551165         IF( ln_wd_il ) THEN 
    1156             zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)  
    1157             zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)  
     1166            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 
     1167            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 
    11581168         ENDIF 
    11591169 
    … …  
    11891199      !!---------------------------------------------------------------------- 
    11901200      ! 
    1191 !!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!   
     1201!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!! 
    11921202      jpi   = size(fsp,1) 
    11931203      jpj   = size(fsp,2) 
    … …  
    13591369   !!====================================================================== 
    13601370END MODULE dynhpg 
    1361  

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynldf_iso.F90

    r12377 r13151  
    2222   USE ldftra          ! lateral physics: eddy diffusivity 
    2323   USE zdf_oce         ! ocean vertical physics 
    24    USE ldfslp          ! iso-neutral slopes  
     24   USE ldfslp          ! iso-neutral slopes 
    2525   ! 
    2626   USE in_out_manager  ! I/O manager 
    … …  
    3636 
    3737   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity 
    38     
    39    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)  
     38 
     39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso) 
    4040   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      - 
    4141 
    4242   !! * Substitutions 
    4343#  include "do_loop_substitute.h90" 
     44#  include "domzgr_substitute.h90" 
    4445   !!---------------------------------------------------------------------- 
    4546   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5354      !!                  ***  ROUTINE dyn_ldf_iso_alloc  *** 
    5455      !!---------------------------------------------------------------------- 
    55       ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     &  
     56      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
    5657         &      akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 
    5758         ! 
    … …  
    6364      !!---------------------------------------------------------------------- 
    6465      !!                     ***  ROUTINE dyn_ldf_iso  *** 
    65       !!                        
     66      !! 
    6667      !! ** Purpose :   Compute the before trend of the rotated laplacian 
    6768      !!      operator of lateral momentum diffusion except the diagonal 
    … …  
    137138         ! 
    138139       ENDIF 
    139           
     140 
    140141      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max 
    141        
     142 
    142143      !                                                ! =============== 
    143144      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    161162 
    162163         !                               -----f----- 
    163          ! Horizontal fluxes on U             |   
     164         ! Horizontal fluxes on U             | 
    164165         ! --------------------===        t   u   t 
    165          !                                    |   
     166         !                                    | 
    166167         ! i-flux at t-point             -----f----- 
    167168 
    168169         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
    169170            DO_2D_00_01 
    170                zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 
     171               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   & 
     172                  &    * MIN( e3u(ji  ,jj,jk,Kmm),                & 
     173                  &           e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 
    171174 
    172175               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     & 
    … …  
    181184         ELSE                   ! other coordinate system (zco or sco) : e3t 
    182185            DO_2D_00_01 
    183                zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 
     186               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   & 
     187                  &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 
    184188 
    185189               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     & 
    … …  
    196200         ! j-flux at f-point 
    197201         DO_2D_10_10 
    198             zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 
     202            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   & 
     203               &     * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 
    199204 
    200205            zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     & 
    … …  
    215220 
    216221         DO_2D_00_10 
    217             zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 
     222            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   & 
     223               &     * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 
    218224 
    219225            zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     & 
    … …  
    230236         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
    231237            DO_2D_01_10 
    232                zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 
     238               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   & 
     239                  &     * MIN( e3v(ji,jj  ,jk,Kmm),                 & 
     240                  &            e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 
    233241 
    234242               zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     & 
    … …  
    243251         ELSE                   ! other coordinate system (zco or sco) : e3t 
    244252            DO_2D_01_10 
    245                zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 
     253               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   & 
     254                  &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 
    246255 
    247256               zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    … …  
    261270         DO_2D_00_00 
    262271            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    & 
    263                &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     272               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   & 
     273               &                           / e3u(ji,jj,jk,Kmm) 
    264274            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    & 
    265                &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     275               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj)   & 
     276               &                           / e3v(ji,jj,jk,Kmm) 
    266277         END_2D 
    267278         !                                             ! =============== 
    … …  
    278289         !                                             ! =============== 
    279290 
    280   
     291 
    281292         ! I. vertical trends associated with the lateral mixing 
    282293         ! ===================================================== 
    … …  
    375386         DO jk = 1, jpkm1 
    376387            DO ji = 2, jpim1 
    377                puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
    378                pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     388               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj)   & 
     389                  &               / e3u(ji,jj,jk,Kmm) 
     390               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj)   & 
     391                  &               / e3v(ji,jj,jk,Kmm) 
    379392            END DO 
    380393         END DO 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynldf_lap_blp.F90

    r12377 r13151  
    1414   USE dom_oce        ! ocean space and time domain 
    1515   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
    16    USE ldfslp         ! iso-neutral slopes  
     16   USE ldfslp         ! iso-neutral slopes 
    1717   USE zdf_oce        ! ocean vertical physics 
    1818   ! 
    … …  
    2828   !! * Substitutions 
    2929#  include "do_loop_substitute.h90" 
     30#  include "domzgr_substitute.h90" 
    3031   !!---------------------------------------------------------------------- 
    3132   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    32    !! $Id$  
     33   !! $Id$ 
    3334   !! Software governed by the CeCILL license (see ./LICENSE) 
    3435   !!---------------------------------------------------------------------- 
    … …  
    3839      !!---------------------------------------------------------------------- 
    3940      !!                     ***  ROUTINE dyn_ldf_lap  *** 
    40       !!                        
    41       !! ** Purpose :   Compute the before horizontal momentum diffusive  
     41      !! 
     42      !! ** Purpose :   Compute the before horizontal momentum diffusive 
    4243      !!      trend and add it to the general trend of momentum equation. 
    4344      !! 
    44       !! ** Method  :   The Laplacian operator apply on horizontal velocity is  
    45       !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
     45      !! ** Method  :   The Laplacian operator apply on horizontal velocity is 
     46      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
    4647      !! 
    4748      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 
    … …  
    7677!!gm open question here : e3f  at before or now ?    probably now... 
    7778!!gm note that ahmf has already been multiplied by fmask 
    78             zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
    79                &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
    80                &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
     79            zcur(ji-1,jj-1) =  & 
     80               &      ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      & 
     81               &  * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
     82               &     - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
    8183            !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    8284!!gm note that ahmt has already been multiplied by tmask 
    83             zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         & 
    84                &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
    85                &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
     85            zdiv(ji,jj)     =   & 
     86               &   ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)      & 
     87               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk)        & 
     88               &        - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
     89               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk)        & 
     90               &        - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
    8691         END_2D 
    8792         ! 
    8893         DO_2D_00_00 
    89             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 & 
    90                &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
    91                &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
     94            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                             & 
     95               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj)   & 
     96               &              / e3u(ji,jj,jk,Kmm)   & 
     97               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)        ) 
    9298               ! 
    93             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                                 & 
    94                &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
    95                &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
     99            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                              & 
     100               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj)   & 
     101               &              / e3v(ji,jj,jk,Kmm)   & 
     102               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)       ) 
    96103         END_2D 
    97104         !                                             ! =============== 
    … …  
    105112      !!---------------------------------------------------------------------- 
    106113      !!                 ***  ROUTINE dyn_ldf_blp  *** 
    107       !!                     
    108       !! ** Purpose :   Compute the before lateral momentum viscous trend  
     114      !! 
     115      !! ** Purpose :   Compute the before lateral momentum viscous trend 
    109116      !!              and add it to the general trend of momentum equation. 
    110117      !! 

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynspg_ts.F90

    r12489 r13151  
    8787   !! * Substitutions 
    8888#  include "do_loop_substitute.h90" 
     89#  include "domzgr_substitute.h90" 
    8990   !!---------------------------------------------------------------------- 
    9091   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    161162      REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points 
    162163      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes 
     164      REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 
    163165      ! 
    164166      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True.  
    … …  
    227229      !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends) 
    228230      !                                   !  ---------------------------  ! 
    229       zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 
    230       zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 
     231      DO jk = 1 , jpk 
     232         ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 
     233         ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 
     234      END DO 
     235      ! 
     236      zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 
     237      zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 
    231238      ! 
    232239      ! 
    … …  
    250257      zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column 
    251258      ! 
    252       CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in 
     259      CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV,  &   ! <<== in 
    253260         &                                                                     zu_trd, zv_trd   )   ! ==>> out 
    254261      ! 
    … …  
    567574         ! at each time step. We however keep them constant here for optimization. 
    568575         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 
    569          CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   ) 
     576         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   ) 
    570577         ! 
    571578         ! Add tidal astronomical forcing if defined 
    … …  
    10881095      ! 
    10891096      SELECT CASE( nvor_scheme ) 
    1090       CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme) 
     1097      CASE( np_EEN )                != EEN scheme using e3f energy & enstrophy scheme 
    10911098         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
    10921099         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    … …  
    11151122         END_2D 
    11161123         ! 
    1117       CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme) 
     1124      CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme 
    11181125         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    11191126         DO_2D_01_01 
    … …  
    11791186 
    11801187 
    1181    SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   ) 
     1188   SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   ) 
    11821189      !!--------------------------------------------------------------------- 
    11831190      !!                   ***  ROUTINE dyn_cor_2d  *** 
    … …  
    11871194      INTEGER  ::   ji ,jj                             ! dummy loop indices 
    11881195      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      - 
    1189       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phu, phv, punb, pvnb, zhU, zhV 
     1196      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pht, phu, phv, punb, pvnb, zhU, zhV 
    11901197      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd 
    11911198      !!---------------------------------------------------------------------- 
    … …  
    11961203            z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
    11971204            zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    & 
    1198                &               * (  e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   & 
    1199                &                  + e1e2t(ji  ,jj)*ht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   ) 
     1205               &               * (  e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   & 
     1206               &                  + e1e2t(ji  ,jj)*pht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   ) 
    12001207               ! 
    12011208            zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    & 
    1202                &               * (  e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   &  
    1203                &                  + e1e2t(ji,jj  )*ht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   )  
     1209               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   &  
     1210               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   )  
    12041211         END_2D 
    12051212         !          

  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynvor.F90

    r12377 r13151  
    1515   !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme 
    1616   !!            3.3  ! 2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase 
    17    !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity  
     17   !!            3.7  ! 2014-04  (G. Madec)  trend simplification: suppress jpdyn_trd_dat vorticity 
    1818   !!             -   ! 2014-06  (G. Madec)  suppression of velocity curl from in-core memory 
    1919   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    … …  
    7070   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
    7171 
    72    INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     72   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
    7373   !                               ! associated indices: 
    7474   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
    … …  
    7979 
    8080   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation 
    81    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -  
     81   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
    8282   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation 
    83    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -  
    84     
     83   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
     84 
    8585   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
    8686   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8 
    8787   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12 
    88     
     88 
    8989   !! * Substitutions 
    9090#  include "do_loop_substitute.h90" 
     91#  include "domzgr_substitute.h90" 
     92 
    9193   !!---------------------------------------------------------------------- 
    9294   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    103105      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 
    104106      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative 
    105       !!               and planetary vorticity trends) and send them to trd_dyn  
     107      !!               and planetary vorticity trends) and send them to trd_dyn 
    106108      !!               for futher diagnostics (l_trddyn=T) 
    107109      !!---------------------------------------------------------------------- 
    … …  
    193195      !!                  ***  ROUTINE vor_enT  *** 
    194196      !! 
    195       !! ** Purpose :   Compute the now total vorticity trend and add it to  
     197      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    196198      !!      the general trend of the momentum equation. 
    197199      !! 
    198       !! ** Method  :   Trend evaluated using now fields (centered in time)  
     200      !! ** Method  :   Trend evaluated using now fields (centered in time) 
    199201      !!       and t-point evaluation of vorticity (planetary and relative). 
    200202      !!       conserves the horizontal kinetic energy. 
    201       !!         The general trend of momentum is increased due to the vorticity  
     203      !!         The general trend of momentum is increased due to the vorticity 
    202204      !!       term which is given by: 
    203205      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ] 
    … …  
    233235                  &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    234236            END_2D 
    235             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     237            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    236238               DO_2D_10_10 
    237239                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    … …  
    248250                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj) 
    249251            END_2D 
    250             IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity  
     252            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
    251253               DO_2D_10_10 
    252254                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 
    … …  
    269271            DO_2D_01_01 
    270272               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   & 
    271                   &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     273                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 
     274                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    272275            END_2D 
    273276         CASE ( np_MET )                           !* metric term 
    274277            DO_2D_01_01 
    275                zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   & 
    276                   &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t(ji,jj,jk,Kmm) 
     278               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)     & 
     279                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) & 
     280                  &             * e3t(ji,jj,jk,Kmm) 
    277281            END_2D 
    278282         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
    279283            DO_2D_01_01 
    280                zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    & 
    281                   &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
     284               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)      & 
     285                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) & 
     286                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 
    282287            END_2D 
    283288         CASE ( np_CME )                           !* Coriolis + metric 
    284289            DO_2D_01_01 
    285                zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           & 
    286                     &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  & 
    287                     &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t(ji,jj,jk,Kmm) 
     290               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                             & 
     291                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)    & 
     292                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) & 
     293                    &          * e3t(ji,jj,jk,Kmm) 
    288294            END_2D 
    289295         CASE DEFAULT                                             ! error 
    … …  
    298304               ! 
    299305            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    & 
    300                &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   &  
    301                &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   )  
     306               &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
     307               &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
    302308         END_2D 
    303309         !                                             ! =============== 
    … …  
    311317      !!                  ***  ROUTINE vor_ene  *** 
    312318      !! 
    313       !! ** Purpose :   Compute the now total vorticity trend and add it to  
     319      !! ** Purpose :   Compute the now total vorticity trend and add it to 
    314320      !!      the general trend of the momentum equation. 
    315321      !! 
    316       !! ** Method  :   Trend evaluated using now fields (centered in time)  
     322      !! ** Method  :   Trend evaluated using now fields (centered in time) 
    317323      !!       and the Sadourny (1975) flux form formulation : conserves the 
    318324      !!       horizontal kinetic energy. 
    319       !!         The general trend of momentum is increased due to the vorticity  
     325      !!         The general trend of momentum is increased due to the vorticity 
    320326      !!       term