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 12191 for NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce.F90 – NEMO

Ignore:
Timestamp:
2019-12-11T16:56:06+01:00 (4 years ago)
Author:
jchanut
Message:

Add dev_AGRIF-01-05_merged branch, e.g. 2019 AGRIF dev

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce.F90

    r10425 r12191  
    1717 
    1818   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 
    19 #if defined key_vertical 
    20    PUBLIC reconstructandremap ! remapping routine 
    21 #endif    
     19   
    2220   !                                              !!* Namelist namagrif: AGRIF parameters 
    23    LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !: 
    24    INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points) 
     21   LOGICAL , PUBLIC ::   ln_agrif_2way = .TRUE.    !: activate two way nesting  
     22   LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !: use zeros (.false.) or not (.true.) in 
     23                                                   !: bdys dynamical fields interpolation 
    2524   REAL(wp), PUBLIC ::   rn_sponge_tra = 2800.     !: sponge coeff. for tracers 
    2625   REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics 
     26   REAL(wp), PUBLIC ::   rn_trelax_tra = 0.01      !: time relaxation parameter for tracers 
     27   REAL(wp), PUBLIC ::   rn_trelax_dyn = 0.01      !: time relaxation parameter for momentum 
    2728   LOGICAL , PUBLIC ::   ln_chk_bathy  = .FALSE.   !: check of parent bathymetry  
    28    LOGICAL , PUBLIC ::   lk_agrif_clp  = .FALSE.   !: Force clamped bcs 
    29    !                                              !!! OLD namelist names 
    30    REAL(wp), PUBLIC ::   visc_tra                  !: sponge coeff. for tracers 
    31    REAL(wp), PUBLIC ::   visc_dyn                  !: sponge coeff. for dynamics 
    32  
     29   ! 
     30   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points) 
    3331   LOGICAL , PUBLIC :: spongedoneT = .FALSE.       !: tracer   sponge layer indicator 
    3432   LOGICAL , PUBLIC :: spongedoneU = .FALSE.       !: dynamics sponge layer indicator 
     
    4240   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u 
    4341   LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v 
    44    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities 
    45    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities 
     42   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utint_stage 
     43   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtint_stage 
     44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspu, fspv !: sponge arrays 
     45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspt, fspf !:   "      " 
    4646 
    4747   ! Barotropic arrays used to store open boundary data during time-splitting loop: 
    48    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy_w, vbdy_w, hbdy_w 
    49    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy_e, vbdy_e, hbdy_e 
    50    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy_n, vbdy_n, hbdy_n 
    51    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy_s, vbdy_s, hbdy_s 
     48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy, vbdy, hbdy 
    5249 
     50# if defined key_vertical 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 
     52   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent 
     53# endif 
    5354 
    5455   INTEGER, PUBLIC :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update 
     
    6465   INTEGER, PUBLIC :: scales_t_id 
    6566   INTEGER, PUBLIC :: avt_id, avm_id, en_id                ! TKE related identificators 
    66    INTEGER, PUBLIC :: umsk_id, vmsk_id 
     67   INTEGER, PUBLIC :: mbkt_id, ht0_id 
    6768   INTEGER, PUBLIC :: kindic_agr 
    6869    
     
    8283      ierr(:) = 0 
    8384      ! 
    84       ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj),   & 
    85          &      fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj),   & 
    86          &      tabspongedone_tsn(jpi,jpj),           & 
     85      ALLOCATE( fspu(jpi,jpj), fspv(jpi,jpj),          & 
     86         &      fspt(jpi,jpj), fspf(jpi,jpj),               & 
     87         &      tabspongedone_tsn(jpi,jpj),                 & 
     88         &      utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & 
    8789# if defined key_top          
    8890         &      tabspongedone_trn(jpi,jpj),           & 
    89 # endif          
     91# endif    
     92# if defined key_vertical 
     93         &      ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj),  & 
     94         &      hu0_parent(jpi,jpj), mbku_parent(jpi,jpj),  & 
     95         &      hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj),  & 
     96# endif       
    9097         &      tabspongedone_u  (jpi,jpj),           & 
    9198         &      tabspongedone_v  (jpi,jpj), STAT = ierr(1) ) 
    9299 
    93       ALLOCATE( ubdy_w(nbghostcells,jpj), vbdy_w(nbghostcells,jpj), hbdy_w(nbghostcells,jpj),   & 
    94          &      ubdy_e(nbghostcells,jpj), vbdy_e(nbghostcells,jpj), hbdy_e(nbghostcells,jpj),   &  
    95          &      ubdy_n(jpi,nbghostcells), vbdy_n(jpi,nbghostcells), hbdy_n(jpi,nbghostcells),   &  
    96          &      ubdy_s(jpi,nbghostcells), vbdy_s(jpi,nbghostcells), hbdy_s(jpi,nbghostcells), STAT = ierr(2) ) 
     100      ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) ) 
    97101 
    98102      agrif_oce_alloc = MAXVAL(ierr) 
     
    100104   END FUNCTION agrif_oce_alloc 
    101105 
    102 #if defined key_vertical 
    103    SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
    104       !!---------------------------------------------------------------------- 
    105       !!                ***  FUNCTION reconstructandremap  *** 
    106       !!---------------------------------------------------------------------- 
    107       IMPLICIT NONE 
    108       INTEGER N, Nout 
    109       REAL(wp) tabin(N), tabout(Nout) 
    110       REAL(wp) hin(N), hout(Nout) 
    111       REAL(wp) coeffremap(N,3),zwork(N,3) 
    112       REAL(wp) zwork2(N+1,3) 
    113       INTEGER jk 
    114       DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8   
    115       REAL(wp) q,q01,q02,q001,q002,q0 
    116       REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 
    117       REAL(wp),PARAMETER :: dpthin = 1.D-3 
    118       INTEGER :: k1, kbox, ktop, ka, kbot 
    119       REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
    120  
    121       z_win(1)=0.; z_wout(1)= 0. 
    122       DO jk=1,N 
    123          z_win(jk+1)=z_win(jk)+hin(jk) 
    124       ENDDO  
    125        
    126       DO jk=1,Nout 
    127          z_wout(jk+1)=z_wout(jk)+hout(jk)        
    128       ENDDO        
    129  
    130       DO jk=2,N 
    131          zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 
    132       ENDDO 
    133          
    134       DO jk=2,N-1 
    135          q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 
    136          zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 
    137          zwork(jk,3)=q0 
    138       ENDDO        
    139       
    140       DO jk= 2,N 
    141          zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 
    142       ENDDO 
    143          
    144       coeffremap(:,1) = tabin(:) 
    145   
    146       DO jk=2,N-1 
    147          q001 = hin(jk)*zwork2(jk+1,1) 
    148          q002 = hin(jk)*zwork2(jk,1)         
    149          IF (q001*q002 < 0) then 
    150             q001 = 0. 
    151             q002 = 0. 
    152          ENDIF 
    153          q=zwork(jk,2) 
    154          q01=q*zwork2(jk+1,1) 
    155          q02=q*zwork2(jk,1) 
    156          IF (abs(q001) > abs(q02)) q001 = q02 
    157          IF (abs(q002) > abs(q01)) q002 = q01 
    158          
    159          q=(q001-q002)*zwork(jk,3) 
    160          q001=q001-q*hin(jk+1) 
    161          q002=q002+q*hin(jk-1) 
    162          
    163          coeffremap(jk,3)=coeffremap(jk,1)+q001 
    164          coeffremap(jk,2)=coeffremap(jk,1)-q002 
    165          
    166          zwork2(jk,1)=(2.*q001-q002)**2 
    167          zwork2(jk,2)=(2.*q002-q001)**2 
    168       ENDDO 
    169          
    170       DO jk=1,N 
    171          IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 
    172             coeffremap(jk,3) = coeffremap(jk,1) 
    173             coeffremap(jk,2) = coeffremap(jk,1) 
    174             zwork2(jk,1) = 0. 
    175             zwork2(jk,2) = 0. 
    176          ENDIF 
    177       ENDDO 
    178          
    179       DO jk=2,N 
    180          q002=max(zwork2(jk-1,2),dsmll) 
    181          q001=max(zwork2(jk,1),dsmll) 
    182          zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 
    183       ENDDO 
    184          
    185       zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
    186       zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
    187   
    188       DO jk=1,N 
    189          q01=zwork2(jk+1,3)-coeffremap(jk,1) 
    190          q02=coeffremap(jk,1)-zwork2(jk,3) 
    191          q001=2.*q01 
    192          q002=2.*q02 
    193          IF (q01*q02<0) then 
    194             q01=0. 
    195             q02=0. 
    196          ELSEIF (abs(q01)>abs(q002)) then 
    197             q01=q002 
    198          ELSEIF (abs(q02)>abs(q001)) then 
    199             q02=q001 
    200          ENDIF 
    201          coeffremap(jk,2)=coeffremap(jk,1)-q02 
    202          coeffremap(jk,3)=coeffremap(jk,1)+q01 
    203       ENDDO 
    204  
    205       zbot=0.0 
    206       kbot=1 
    207       DO jk=1,Nout 
    208          ztop=zbot  !top is bottom of previous layer 
    209          ktop=kbot 
    210          IF     (ztop.GE.z_win(ktop+1)) then 
    211             ktop=ktop+1 
    212          ENDIF 
    213          
    214          zbot=z_wout(jk+1) 
    215          zthk=zbot-ztop 
    216  
    217          IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 
    218  
    219             kbot=ktop 
    220             DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
    221                kbot=kbot+1 
    222             ENDDO 
    223             zbox=zbot 
    224             DO k1= jk+1,Nout 
    225                IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 
    226                   exit !thick layer 
    227                ELSE 
    228                   zbox=z_wout(k1+1)  !include thin adjacent layers 
    229                   IF(zbox.EQ.z_wout(Nout+1)) THEN 
    230                      exit !at bottom 
    231                   ENDIF 
    232                ENDIF 
    233             ENDDO 
    234             zthk=zbox-ztop 
    235  
    236             kbox=ktop 
    237             DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
    238                kbox=kbox+1 
    239             ENDDO 
    240            
    241             IF(ktop.EQ.kbox) THEN 
    242                IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 
    243                   IF(hin(kbox).GT.dpthin) THEN 
    244                      q001 = (zbox-z_win(kbox))/hin(kbox) 
    245                      q002 = (ztop-z_win(kbox))/hin(kbox) 
    246                      q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
    247                      q02=q01-1.+(q001+q002) 
    248                      q0=1.-q01-q02 
    249                   ELSE 
    250                      q0 = 1.0 
    251                      q01 = 0. 
    252                      q02 = 0. 
    253                   ENDIF 
    254                   tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
    255                ELSE 
    256                   tabout(jk) = tabin(kbox) 
    257                ENDIF  
    258             ELSE 
    259                IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 
    260                   ka = jk 
    261                ELSEIF (kbox-ktop.GE.3) THEN 
    262                   ka = (kbox+ktop)/2 
    263                ELSEIF (hin(ktop).GE.hin(kbox)) THEN 
    264                   ka = ktop 
    265                ELSE 
    266                   ka = kbox 
    267                ENDIF !choose ka 
    268      
    269                offset=coeffremap(ka,1) 
    270      
    271                qtop = z_win(ktop+1)-ztop !partial layer thickness 
    272                IF(hin(ktop).GT.dpthin) THEN 
    273                   q=(ztop-z_win(ktop))/hin(ktop) 
    274                   q01=q*(q-1.) 
    275                   q02=q01+q 
    276                   q0=1-q01-q02             
    277                ELSE 
    278                   q0 = 1. 
    279                   q01 = 0. 
    280                   q02 = 0. 
    281                ENDIF 
    282                 
    283                tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
    284      
    285                DO k1= ktop+1,kbox-1 
    286                   tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
    287                ENDDO !k1 
    288                 
    289                qbot = zbox-z_win(kbox) !partial layer thickness 
    290                IF(hin(kbox).GT.dpthin) THEN 
    291                   q=qbot/hin(kbox) 
    292                   q01=(q-1.)**2 
    293                   q02=q01-1.+q 
    294                   q0=1-q01-q02                             
    295                ELSE 
    296                   q0 = 1.0 
    297                   q01 = 0. 
    298                   q02 = 0. 
    299                ENDIF 
    300                
    301                tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
    302                 
    303                rpsum=1.0d0/zthk 
    304                tabout(jk)=offset+tsum*rpsum 
    305                   
    306             ENDIF !single or multiple layers 
    307          ELSE 
    308             IF (jk==1) THEN 
    309                write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 
    310             ENDIF 
    311             tabout(jk) = tabout(jk-1) 
    312               
    313          ENDIF !normal:thin layer 
    314       ENDDO !jk 
    315              
    316       return 
    317       end subroutine reconstructandremap 
    318 #endif 
    319  
    320106#endif 
    321107   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.