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 15279 for utils/tools/DOMAINcfg/src/agrif_dom_update.F90 – NEMO

Ignore:
Timestamp:
2021-09-23T12:00:23+02:00 (3 years ago)
Author:
jchanut
Message:

#2222 and #2638: Enable creating agrif meshes with different vertical grids (geopotential only as a start)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/DOMAINcfg/src/agrif_dom_update.F90

    r14634 r15279  
    22 
    33   USE dom_oce 
    4    USE domzgr 
    54   USE agrif_parameters 
    65   USE agrif_profiles 
    7    USE lbclnk 
    8     
     6   USE agrif_recompute_scales 
     7  
    98   IMPLICIT none 
    109   PRIVATE 
     
    2120      !!----------------------------------------------------------------------   
    2221      ! 
     22      INTEGER :: ind1, ind2 
     23 
    2324      IF( Agrif_Root() ) return 
    24  
    25       CALL agrif_update_variable(bottom_level_id,procname = update_bottom_level) 
    26       ! 
    27       Agrif_UseSpecialValueInUpdate = .TRUE. 
    28       Agrif_SpecialValueFineGrid    = 0._wp          
    29       CALL agrif_update_variable(e3t_id,procname = update_e3t) 
     25       
     26      IF ( .NOT.ln_vert_remap ) THEN 
     27         CALL agrif_update_variable(bottom_level_id,procname = update_bottom_level) 
     28         Agrif_UseSpecialValueInUpdate = .FALSE. 
     29         Agrif_SpecialValueFineGrid    = 0._wp          
     30         CALL agrif_update_variable(e3t_id, procname = update_e3t_z)  
     31         ! 
     32      ELSE 
     33         Agrif_UseSpecialValueInUpdate = .FALSE. 
     34         Agrif_SpecialValueFineGrid    = 0._wp          
     35         CALL agrif_update_variable(e3t_id, procname = update_e3t_z_cons)  
     36 
     37         ! jc: extend update zone outside dynamical interface within sponge zone: 
     38         ! Use max operator this time to account for cases for which Agrif_Rho > nbghostcells 
     39         ind1 = CEILING(REAL(max(nbghostcells_x_w-1, nbghostcells_x_e-1), wp) / Agrif_Rhox() ) 
     40         ind2 = CEILING(REAL(max(nbghostcells_y_s-1, nbghostcells_y_n-1), wp) / Agrif_Rhoy() ) 
     41         CALL agrif_update_variable(e3t_copy_id, locupdate1=(/-ind1,0/), & 
     42                             &                   locupdate2=(/-ind2,0/),procname = update_e3t_z_cons) 
     43      ENDIF 
    3044      Agrif_UseSpecialValueInUpdate = .FALSE. 
     45      ! 
     46      ! Update vertical scale factors at U, V and F-points: 
     47      CALL Agrif_ChildGrid_To_ParentGrid() 
     48      CALL agrif_recompute_scalefactors 
     49      CALL Agrif_ParentGrid_To_ChildGrid() 
    3150      !     
    3251   END SUBROUTINE agrif_update_all 
    3352 
    34    SUBROUTINE update_bottom_level( ptab, i1, i2, j1, j2, before, nb,ndir) 
     53   SUBROUTINE update_bottom_level( ptab, i1, i2, j1, j2, before) 
    3554      !!---------------------------------------------------------------------- 
    36       !!                  ***  ROUTINE interpsshn  *** 
     55      !!       ***  ROUTINE update_bottom_level  *** 
    3756      !!----------------------------------------------------------------------   
    3857      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    3958      REAL, DIMENSION(i1:i2,j1:j2)    , INTENT(inout) ::   ptab 
    4059      LOGICAL                         , INTENT(in   ) ::   before 
    41       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    4260      ! 
    4361      !!----------------------------------------------------------------------  
    44       REAL(WP),DIMENSION(jpi,jpj) :: zk 
    4562      ! 
    4663      IF( before) THEN 
     
    5067          
    5168         WHERE ( mbkt(i1:i2,j1:j2) .EQ. 0 ) 
    52             ssmask(i1:i2,j1:j2) = 0. 
    53             mbkt(i1:i2,j1:j2)   = 1 
     69            ssmask(i1:i2,j1:j2) = 0._wp 
     70            mbkt(i1:i2,j1:j2)   = 1  
    5471         ELSEWHERE 
    55             ssmask(i1:i2,j1:j2) = 1. 
     72            ssmask(i1:i2,j1:j2) = 1._wp 
    5673         END WHERE  
    57 !         zk(:,:) = REAL(mbkt(:,:),wp); CALL lbc_lnk('update_bottom',zk,'T',1.); mbkt(:,:) = MAX(NINT(zk(:,:)),1) 
    58 !         CALL lbc_lnk('update_bottom',ssmask,'T',1.)           
    5974      ENDIF 
    6075      ! 
    6176   END SUBROUTINE update_bottom_level 
    62     
    63    SUBROUTINE update_e3t( tabres, i1, i2, j1, j2, k1, k2, before ) 
    64       !!--------------------------------------------- 
    65       !!           *** update_e3t *** 
     77 
     78   SUBROUTINE update_e3t_z( tabres, i1, i2, j1, j2, k1, k2, before ) 
     79      !!--------------------------------------------- 
     80      !!           *** update_e3t_z *** 
    6681      !!--------------------------------------------- 
    6782      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
     
    6984      LOGICAL, INTENT(in) :: before 
    7085      !! 
    71       INTEGER :: ji,jj,jk 
     86      INTEGER :: ji, jj, jk 
    7287      !!--------------------------------------------- 
    7388      ! 
     
    7691            DO jj=j1,j2 
    7792               DO ji=i1,i2 
    78                    IF ((ssmask(ji,jj) /=0.).AND.( mbkt(ji,jj) .GE. jk )) THEN 
     93                  IF ( (ssmask(ji,jj) /=0._wp).AND.(mbkt(ji,jj).GE.jk) ) THEN 
     94                     tabres(ji,jj,jk) = e3t_0(ji,jj,jk) 
     95                  ELSE 
     96                     tabres(ji,jj,jk) = 0._wp 
     97                  ENDIF  
     98               END DO 
     99            END DO 
     100         END DO 
     101      ELSE 
     102         DO jk=1,jpk 
     103            DO jj=j1,j2 
     104               DO ji=i1,i2 
     105                  IF ( ( mbkt(ji,jj).GE.jk ).AND.(ssmask(ji,jj)==1._wp) ) THEN 
     106                     e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
     107                 !    e3t_0(ji,jj,jk) = tabres(ji,jj,jk) 
     108                  ELSE 
     109                     e3t_0(ji,jj,jk) = e3t_1d(jk) 
     110                  ENDIF 
     111               END DO 
     112            END DO 
     113         END DO 
     114         ! 
     115      ENDIF 
     116      !  
     117   END SUBROUTINE update_e3t_z 
     118 
     119   SUBROUTINE update_e3t_z_cons( tabres, i1, i2, j1, j2, k1, k2, before ) 
     120      !!--------------------------------------------- 
     121      !!           *** update_e3t_z_cons *** 
     122      !!--------------------------------------------- 
     123      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
     124      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
     125      LOGICAL, INTENT(in) :: before 
     126      !! 
     127      INTEGER :: ji, jj, jk, ik 
     128      REAL(wp) :: zhmin, zdepth, zdepwp, ze3tp, ze3wp 
     129      !!--------------------------------------------- 
     130      ! 
     131      IF (before) THEN 
     132         DO jk = k1, k2-1 
     133            DO jj = j1, j2 
     134               DO ji = i1, i2 
     135                   IF ( (ssmask(ji,jj) /=0._wp).AND.( mbkt(ji,jj) .GE. jk ) ) THEN 
    79136                      tabres(ji,jj,jk) = e3t_0(ji,jj,jk) 
    80137                   ELSE 
    81                       tabres(ji,jj,jk) = 0. 
     138                      tabres(ji,jj,jk) = 0._wp 
    82139                   endif 
    83140               END DO 
    84141            END DO 
    85142         END DO 
    86       ELSE 
    87          DO jk=k1,k2 
    88             DO jj=j1,j2 
    89                DO ji=i1,i2 
    90                    IF( mbkt(ji,jj) .GE. jk ) THEN 
    91                       e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
    92                   !    e3t_0(ji,jj,jk) = tabres(ji,jj,jk) 
    93                    ELSE 
    94                       e3t_0(ji,jj,jk) = e3t_1d(jk) 
    95                    ENDIF 
    96                END DO 
    97             END DO 
    98          END DO 
    99  
    100 !         CALL lbc_lnk('update_e3t',e3t_0,'T',1.,kfillmode = jpfillcopy) 
    101          ! 
     143         tabres(i1:i2,j1:j2,k2) = ssmask(i1:i2,j1:j2) ! To get fractional area 
     144      ELSE 
     145         IF( rn_hmin < 0._wp ) THEN    
     146            ik = - INT( rn_hmin ) 
     147         ELSE                           
     148            ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) 
     149         ENDIF 
     150         zhmin = gdepw_1d(ik+1) 
     151 
     152         ! Compute child bathymetry: 
     153         bathy(i1:i2,j1:j2) = 0._wp 
     154         DO jk=k1,k2-1    
     155            bathy(i1:i2,j1:j2) = bathy(i1:i2,j1:j2) + tabres(i1:i2,j1:j2,jk) 
     156         END DO 
     157         WHERE( bathy(i1:i2,j1:j2) == 0._wp )   ;   mbathy(i1:i2,j1:j2) = 0        
     158         ELSE WHERE                             ;   mbathy(i1:i2,j1:j2) = jpkm1   
     159         END WHERE 
     160 
     161         DO jk = jpkm1, 1, -1 
     162            zdepth = gdepw_1d(jk) ! + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     163            WHERE( 0._wp < bathy(i1:i2,j1:j2) .AND. bathy(i1:i2,j1:j2) <= zdepth )   mbathy(i1:i2,j1:j2) = jk-1 
     164         END DO 
     165 
     166         ! Scale factors and depth at T- and W-points 
     167         DO jk = 1, jpk   
     168            gdept_0(i1:i2,j1:j2,jk) = gdept_1d(jk) 
     169            gdepw_0(i1:i2,j1:j2,jk) = gdepw_1d(jk) 
     170            e3t_0  (i1:i2,j1:j2,jk) = e3t_1d  (jk) 
     171            e3w_0  (i1:i2,j1:j2,jk) = e3w_1d  (jk) 
     172         END DO 
     173         ! Scale factors and depth at T- and W-points 
     174         DO jj = j1, j2 
     175            DO ji = i1, i2  
     176               ik = mbathy(ji,jj) 
     177               IF( ik > 0 ) THEN               ! ocean point only 
     178                  ! max ocean level case 
     179                  IF( ik == jpkm1 ) THEN 
     180                     zdepwp = bathy(ji,jj) 
     181                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     182                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     183                     e3t_0(ji,jj,ik  ) = ze3tp 
     184                     e3t_0(ji,jj,ik+1) = ze3tp 
     185                     e3w_0(ji,jj,ik  ) = ze3wp 
     186                     e3w_0(ji,jj,ik+1) = ze3tp 
     187                     gdepw_0(ji,jj,ik+1) = zdepwp 
     188                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     189                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     190                     ! 
     191                  ELSE                         ! standard case 
     192                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  
     193                        gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     194                     ELSE                                        
     195                        gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     196                     ENDIF 
     197                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )           & 
     198                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )           & 
     199                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     200                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik))           & 
     201                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik))  
     202                     e3w_0(ji,jj,ik) =                                                                   &  
     203                        &      0.5_wp * (gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     204                        &             * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     205                     !       ... on ik+1 
     206                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     207                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     208                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     209                  ENDIF 
     210               ENDIF 
     211            END DO 
     212         END DO 
     213         ! 
     214         DO jj=j1,j2 
     215            DO ji=i1,i2 
     216               bathy(ji,jj) = SUM(e3t_0(ji,jj,1:mbkt(ji,jj) ) )  
     217            END DO 
     218         END DO 
     219         ! 
     220         WHERE ( ( mbathy(i1:i2,j1:j2) .EQ. 0 )  &  
     221           & .OR.(tabres(i1:i2,j1:j2,k2)<0.5_wp) & 
     222           & .OR.(bathy(i1:i2,j1:j2)<zhmin) ) 
     223            ssmask(i1:i2,j1:j2) = 0._wp 
     224            mbathy(i1:i2,j1:j2) = 0  
     225         ELSEWHERE 
     226            ssmask(i1:i2,j1:j2) = 1._wp 
     227         END WHERE 
     228         mbkt(i1:i2,j1:j2) = MAX( mbathy(i1:i2,j1:j2), 1 ) 
    102229      ENDIF 
    103230      !  
    104    END SUBROUTINE update_e3t 
     231   END SUBROUTINE update_e3t_z_cons 
    105232       
    106233#else 
Note: See TracChangeset for help on using the changeset viewer.