Ignore:
Timestamp:
2020-07-01T16:09:00+02:00 (3 months ago)
Author:
gsamson
Message:

merge with trunk@r13136 with a more recent svn version; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_istate.F90

    r12489 r13197  
    2828   PUBLIC   usr_def_istate   ! called by istate.F90 
    2929 
     30   !! * Substitutions 
     31#  include "do_loop_substitute.h90" 
    3032   !!---------------------------------------------------------------------- 
    3133   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    164166         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1) 
    165167         DO jl=1, jpnj 
    166             DO jj=nldj, nlej 
    167                DO ji=nldi, nlei 
    168                   pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
    169                END DO 
    170             END DO 
     168            DO_2D_00_00 
     169               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
     170            END_2D 
    171171            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    172172         END DO 
     
    183183      CASE(4)    ! geostrophic zonal pulse 
    184184    
    185          DO jj=1, jpj 
    186             DO ji=1, jpi 
    187                IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
    188                   zdu = rn_uzonal 
    189                ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
    190                   zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
    191                ELSE 
    192                   zdu = 0. 
    193                END IF 
    194                IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
    195                   pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
    196                   pu(ji,jj,:) = zdu 
    197                   pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
    198                ELSE 
    199                   pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
    200                   pu(ji,jj,:) = 0. 
    201                   pts(ji,jj,:,jp_sal) = 1. 
    202                END IF 
    203             END DO 
    204          END DO 
     185         DO_2D_11_11 
     186            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
     187               zdu = rn_uzonal 
     188            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
     189               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
     190            ELSE 
     191               zdu = 0. 
     192            END IF 
     193            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
     194               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
     195               pu(ji,jj,:) = zdu 
     196               pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
     197            ELSE 
     198               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
     199               pu(ji,jj,:) = 0. 
     200               pts(ji,jj,:,jp_sal) = 1. 
     201            END IF 
     202         END_2D 
    205203          
    206204         ! temperature: 
    207205         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)         
    208206         pv(:,:,:) = 0. 
    209           
    210207          
    211208       CASE(5)    ! vortex 
     
    220217         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 
    221218         ! 
    222          DO jj=1, jpj 
    223             DO ji=1, jpi 
    224                zx = glamt(ji,jj) * 1.e3 
    225                zy = gphit(ji,jj) * 1.e3 
    226                ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
    227                zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
    228                ! Sea level: 
    229                pssh(ji,jj) = 0. 
    230                DO jl=1,5 
    231                   zdt = pssh(ji,jj) 
    232                   zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
    233                   zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
    234                   pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
    235                END DO 
    236                ! temperature: 
    237                DO jk=1,jpk 
    238                   zdt =  pdept(ji,jj,jk)  
    239                   zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
    240                   IF (zdt < zH) THEN 
    241                      zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
    242                      zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
    243                   ENDIF 
    244                   !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    245                   pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    246                END DO 
    247             END DO 
    248          END DO 
     219         DO_2D_11_11 
     220            zx = glamt(ji,jj) * 1.e3 
     221            zy = gphit(ji,jj) * 1.e3 
     222            ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
     223            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
     224            ! Sea level: 
     225            pssh(ji,jj) = 0. 
     226            DO jl=1,5 
     227               zdt = pssh(ji,jj) 
     228               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
     229               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     230               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
     231            END DO 
     232            ! temperature: 
     233            DO jk=1,jpk 
     234               zdt =  pdept(ji,jj,jk)  
     235               zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
     236               IF (zdt < zH) THEN 
     237                  zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
     238                  zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     239               ENDIF 
     240               !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     241               pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     242            END DO 
     243         END_2D 
    249244         ! 
    250245         ! salinity:   
     
    253248         ! velocities: 
    254249         za = 2._wp * zP0 / zlambda**2 
    255          DO jj=1, jpj 
    256             DO ji=1, jpim1 
    257                zx = glamu(ji,jj) * 1.e3 
    258                zy = gphiu(ji,jj) * 1.e3 
    259                DO jk=1, jpk 
    260                   zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
    261                   IF (zdu < zH) THEN 
    262                      zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
    263                      zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
    264                      pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
    265                   ELSE 
    266                      pu(ji,jj,jk) = 0._wp 
    267                   ENDIF 
    268                END DO 
    269             END DO 
    270          END DO 
    271          ! 
    272          DO jj=1, jpjm1 
    273             DO ji=1, jpi 
    274                zx = glamv(ji,jj) * 1.e3 
    275                zy = gphiv(ji,jj) * 1.e3 
    276                DO jk=1, jpk 
    277                   zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
    278                   IF (zdv < zH) THEN 
    279                      zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
    280                      zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
    281                      pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
    282                   ELSE 
    283                      pv(ji,jj,jk) = 0._wp 
    284                   ENDIF 
    285                END DO 
    286             END DO 
    287          END DO 
     250         DO_2D_00_00 
     251            zx = glamu(ji,jj) * 1.e3 
     252            zy = gphiu(ji,jj) * 1.e3 
     253            DO jk=1, jpk 
     254               zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
     255               IF (zdu < zH) THEN 
     256                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
     257                  zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
     258                  pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
     259               ELSE 
     260                  pu(ji,jj,jk) = 0._wp 
     261               ENDIF 
     262            END DO 
     263         END_2D 
     264         ! 
     265         DO_2D_00_00 
     266            zx = glamv(ji,jj) * 1.e3 
     267            zy = gphiv(ji,jj) * 1.e3 
     268            DO jk=1, jpk 
     269               zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
     270               IF (zdv < zH) THEN 
     271                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
     272                  zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
     273                  pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
     274               ELSE 
     275                  pv(ji,jj,jk) = 0._wp 
     276               ENDIF 
     277            END DO 
     278         END_2D 
    288279         !             
    289280      END SELECT 
    290  
     281       
    291282      IF (ln_sshnoise) THEN 
    292283         CALL RANDOM_NUMBER(zrandom) 
     
    294285      END IF 
    295286      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    296       CALL lbc_lnk(  'usrdef_istate', pts, 'T',  1. ) 
    297       CALL lbc_lnk(   'usrdef_istate', pu, 'U', -1. ) 
    298       CALL lbc_lnk(   'usrdef_istate', pv, 'V', -1. ) 
     287      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. ) 
     288      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 
    299289 
    300290   END SUBROUTINE usr_def_istate 
Note: See TracChangeset for help on using the changeset viewer.