Changeset 9813


Ignore:
Timestamp:
2018-06-19T21:54:00+02:00 (2 years ago)
Author:
antsia
Message:

delete iscplhsb, add iscpldiv and make the code readable

Location:
branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
1 deleted
5 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90

    r9630 r9813  
    2828   REAL(wp), PUBLIC                                        ::   rdt_iscpl 
    2929   !!                                                      !!* namsbc_iscpl namelist * 
    30    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::   hdiv_iscpl 
    31    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   htsc_iscpl 
     30   !!------array used for divergence correction 
     31   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  rhdivdiff 
     32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  rhdivdiff_trac 
     33 
    3234   !! * Substitutions   
    3335#  include "domzgr_substitute.h90"   
     
    4345      !!                ***  ROUTINE sbc_iscpl_alloc  *** 
    4446      !!---------------------------------------------------------------------- 
    45       ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=iscpl_alloc ) 
     47      ALLOCATE( rhdivdiff(jpi,jpj,jpk), rhdivdiff_trac(jpi,jpj,jpk,2), STAT=iscpl_alloc ) 
    4648         ! 
    4749      IF( lk_mpp          )   CALL mpp_sum ( iscpl_alloc ) 
     
    5153   SUBROUTINE iscpl_init() 
    5254      INTEGER ::   ios           ! Local integer output status for namelist read 
     55      INTEGER :: ierr 
    5356      NAMELIST/namsbc_iscpl/nn_fiscpl,ln_hsb,nn_drown 
    5457      !!---------------------------------------------------------------------- 
     
    8386      END IF 
    8487 
     88      IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'iscpl_init : unable to allocate arrays' ) 
     89 
    8590   END SUBROUTINE iscpl_init 
    8691 
  • branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90

    r9630 r9813  
    2424   USE lbclnk          ! communication 
    2525   USE iscplini        ! ice sheet coupling: initialisation 
    26    USE iscplhsb        ! ice sheet coupling: conservation 
     26   !USE iscplhsb        ! ice sheet coupling: conservation 
     27   USE iscpldiv 
    2728 
    2829   IMPLICIT NONE 
     
    8081      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b ) 
    8182 
    82       !! compute correction if conservation needed 
    83       IF ( ln_hsb ) THEN 
    84          IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' ) 
    85          CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl) 
    86       END IF 
    87           
    8883      !! print mesh/mask 
    8984      IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file 
    90  
    91       IF ( ln_hsb ) THEN 
    92          cfile='correction' 
    93          cfile = TRIM( cfile ) 
    94          CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
    95          CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) ) 
    96          CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) ) 
    97          CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) ) 
    98          CALL iom_close ( inum0 ) 
    99       END IF 
    10085 
    10186      CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b )   
     
    150135      !! 
    151136      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb, zhdivdiff, zfp_ui, zfm_ui, zfp_vi, zfm_vi 
    152       REAL(wp):: zdz, zdzm1, zdzp1, zvol 
     137      REAL(wp):: zdz, zdzm1, zdzp1, zvol, zuflux_sum, zvflux_sum 
    153138      REAL(wp),DIMENSION(2)::ztrdtsb 
    154139      !! 
     
    156141      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn 
    157142      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1 
    158       REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1 
    159       REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp, zuflux, zvflux 
     143      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1, isfmask 
     144      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp, zuflux, zvflux, zhdiv, zhdiv2 
    160145      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d 
    161146      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0 
     
    165150      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   ) 
    166151      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        )  
    167       CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb, zuflux, zvflux                       )  
    168       CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       )  
     152      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb, zuflux, zvflux, zhdiv, zhdiv2               )  
     153      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1, isfmask                       )  
    169154      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t)  
    170155      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         )  
     
    277262      END IF 
    278263 
    279 !============================================================================= 
    280 ! compute velocity 
    281  
    282       ub(:,:,:)=un(:,:,:) 
    283       vb(:,:,:)=vn(:,:,:) 
    284  
    285       un(:,:,:)=ub(:,:,:)*umask(:,:,:) 
    286       vn(:,:,:)=vb(:,:,:)*vmask(:,:,:) 
    287      
    288  
    289       !-----initialise the divergence and tracer correction------ 
    290       rhdivdiff(:,:,:) = 0._wp 
    291       rhdivdiff_trac(:,:,:,:) = 0._wp       
    292  
    293       !-----the change of flux in each velocity cell  
    294       zuflux(:,:,:) = (pe3u_b(:,:,:)*pumask_b(:,:,:) - fse3u_n(:,:,:)*umask(:,:,:))*ub(:,:,:) 
    295       zvflux(:,:,:) = (pe3v_b(:,:,:)*pvmask_b(:,:,:) - fse3v_n(:,:,:)*vmask(:,:,:))*vb(:,:,:) 
    296       
    297  
    298       ! compute divergence and tracer correction (whether the horizontal velocity cell is closed or open)     
    299       DO jj = 1,jpj 
    300          DO ji = 1,jpim1   ! jpim1 rather than jpi for the correct lbc_lnk synchronisation 
    301             !--------u-cell is closed-----------   
    302             DO jk  = 1,jpk 
    303                  !------the contribution to the west side of u-cell 
    304                  IF(tmask(ji,jj,jk) == 1) THEN   ! the west cell is still active  
    305                     zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk) 
    306                     rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol  
    307                     rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zuflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e2u(ji,jj)/zvol 
    308                  ELSEIF(mikt(ji,jj) > 1) THEN  ! the west cell becomes dry, put the waterflux to the top of the west column 
    309                     jtop=mikt(ji,jj) 
    310                     zvol=e1e2t(ji,jj)*fse3t(ji,jj,jtop) 
    311                     rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol 
    312                     rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2)-zuflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e2u(ji,jj)/zvol 
    313                  !ELSE ! for the closing T-column : nothing to be done  
    314                  ENDIF 
    315  
    316                  !------the contribution to the east side of u-cell 
    317                  IF(tmask(ji+1,jj,jk) == 1) THEN  ! the east cell is still active 
    318                     zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jk) 
    319                     rhdivdiff(ji+1,jj,jk) = rhdivdiff(ji+1,jj,jk) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol 
    320                     rhdivdiff_trac(ji+1,jj,jk,1:2) = rhdivdiff_trac(ji+1,jj,jk,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jk,1:2)*e2u(ji,jj)/zvol 
    321                  ELSEIF(mikt(ji+1,jj) > 1) THEN    ! the east cell becomes dry, put the waterflux to the top of the east column 
    322                     jtop = mikt(ji+1,jj) 
    323                     zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jtop) 
    324                     rhdivdiff(ji+1,jj,jtop) = rhdivdiff(ji+1,jj,jtop) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol 
    325                     rhdivdiff_trac(ji+1,jj,jtop,1:2) = rhdivdiff_trac(ji+1,jj,jtop,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jtop,1:2)*e2u(ji,jj)/zvol 
    326                  !ELSE  ! for the closing T-column : nothing to be done 
    327                  ENDIF 
    328             ENDDO 
    329          END DO 
    330       END DO 
    331  
    332  
    333       DO jj = 1,jpjm1  ! jpim1 rather than jpi for the correct lbc_lnk synchronisation 
    334          DO ji = 1,jpi    
    335             !-----------v-cell is closed------------- 
    336             DO jk  = 1,jpk 
    337                  !------the contribution to the south side of v-cell 
    338                  IF(tmask(ji,jj,jk) == 1) THEN  ! the south cell is still active 
    339                    zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk) 
    340                    rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol 
    341                    rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e1v(ji,jj)/zvol 
    342                  ELSEIF(mikt(ji,jj) > 1) THEN  !the south cell becomes dry, put the waterflux to the top of the south column 
    343                    jtop = mikt(ji,jj) 
    344                    zvol = e1e2t(ji,jj)*fse3t(ji,jj,jtop) 
    345                    rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol 
    346                    rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e1v(ji,jj)/zvol 
    347                  ENDIF 
    348  
    349                  !------the contribution to the north side of v-cell 
    350                  IF(tmask(ji,jj+1,jk) == 1) THEN ! the north cell is still active 
    351                     zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jk) 
    352                     rhdivdiff(ji,jj+1,jk) = rhdivdiff(ji,jj+1,jk) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol 
    353                     rhdivdiff_trac(ji,jj+1,jk,1:2) = rhdivdiff_trac(ji,jj+1,jk,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jk,1:2)*e1v(ji,jj)/zvol 
    354                  ELSEIF(mikt(ji,jj+1) > 1) THEN  ! the north cell becomes dry, put the waterflux to the top of the north column 
    355                     jtop = mikt(ji,jj+1) 
    356                     zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jtop) 
    357                     rhdivdiff(ji,jj+1,jtop) = rhdivdiff(ji,jj+1,jtop) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol 
    358                     rhdivdiff_trac(ji,jj+1,jtop,1:2) = rhdivdiff_trac(ji,jj+1,jtop,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jtop,1:2)*e1v(ji,jj)/zvol 
    359                  ENDIF 
    360             ENDDO 
    361  
    362          END DO 
    363       END DO  
    364        
    365  
    366       CALL lbc_lnk(rhdivdiff,   'T',1._wp) 
    367       CALL lbc_lnk(rhdivdiff_trac(:,:,:,1),   'T',1._wp) 
    368       CALL lbc_lnk(rhdivdiff_trac(:,:,:,2),   'T',1._wp) 
    369  
     264      !----compute and save the divergence difference (before and after remeshing)  
     265      CALL iscpl_div_corr(ptmask_b, pe3t_b, pe3u_b, pe3v_b) 
    370266 
    371267      ! compute temp and salt 
     
    473369      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   ) 
    474370      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                )  
    475       CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d, zuflux, zvflux    )  
    476       CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       )  
     371      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d, zuflux, zvflux, zhdiv, zhdiv2    )  
     372      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1, isfmask                     )  
    477373      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t)  
    478374      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        )  
  • branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r9630 r9813  
    3535   USE wrk_nemo        ! Memory Allocation 
    3636   USE timing          ! Timing 
    37    USE iscplhsb 
     37   USE iscpldiv 
    3838 
    3939   IMPLICIT NONE 
  • branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r9630 r9813  
    3333   USE timing          ! Timing 
    3434   USE eosbn2 
     35   USE iscpldiv 
    3536#if defined key_asminc    
    3637   USE asminc          ! Assimilation increment 
     
    223224      !---start the tendency with tracer source/sink from divergence correction 
    224225      IF ( ln_iscpl .AND. kt == nit000) THEN         
    225          IF (lwp) WRITE(numout,*) 'add tendency tsa due to re-meshing at kt = nit000' 
    226          tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + rhdivdiff_trac(:,:,:,jp_tem)*tmask(:,:,:) 
    227          tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + rhdivdiff_trac(:,:,:,jp_sal)*tmask(:,:,:) 
     226         CALL iscpl_div_tra( tsa ) 
    228227      ENDIF 
    229228 
  • branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r9630 r9813  
    7171   !! Energy budget of the leads (open water embedded in sea ice) 
    7272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraqsr_1lev        !: fraction of solar net radiation absorbed in the first ocean level [-] 
    73  
    74   !!------array used for divergence correction 
    75    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  rhdivdiff 
    76    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  rhdivdiff_trac 
    7773  
    7874   !! Arrays used in coupling when MEDUSA is present. These arrays need to be declared 
     
    133129         &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     & 
    134130         &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     & 
    135          &     riceload(jpi,jpj),  rhdivdiff(jpi,jpj,jpk), rhdivdiff_trac(jpi,jpj,jpk,2) ,  STAT=ierr(2) ) 
     131         &     riceload(jpi,jpj),                             STAT=ierr(2) ) 
    136132         ! 
    137133      ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) ) 
Note: See TracChangeset for help on using the changeset viewer.