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 13540 for NEMO/branches/2020/r12377_ticket2386/tests/CANAL/MY_SRC/usrdef_istate.F90 – NEMO

Ignore:
Timestamp:
2020-09-29T12:41:06+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2386: update to latest trunk

Location:
NEMO/branches/2020/r12377_ticket2386
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/r12377_ticket2386

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@13507        sette 

  • NEMO/branches/2020/r12377_ticket2386/tests/CANAL/MY_SRC/usrdef_istate.F90

    r12511 r13540  
    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) 
    … …  
    6466      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   ' 
    6567      ! 
    66       IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom) 
    6768      zjetx = ABS(rn_ujetszx)/2. 
    6869      zjety = ABS(rn_ujetszy)/2. 
    6970      ! 
     71      zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 ) 
     72      ! 
    7073      SELECT CASE(nn_initcase) 
     74 
     75      CASE(-1)    ! stratif at rest 
     76 
     77         ! sea level: 
     78         pssh(:,:) = 0. 
     79         ! temperature: 
     80         pts(:,:,1,jp_tem) = 25. !!30._wp 
     81         pts(:,:,2:jpk,jp_tem) = 22. !!24._wp 
     82         ! salinity:   
     83         pts(:,:,:,jp_sal) = 35._wp 
     84         ! velocities: 
     85         pu(:,:,:) = 0. 
     86         pv(:,:,:) = 0. 
     87 
    7188      CASE(0)    ! rest 
    7289          
    … …  
    96113            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra 
    97114            WHERE( ABS(gphit) <= zjety ) 
    98                pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
    99             ELSEWHERE 
    100                pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3   & 
     115               pssh(:,:) = - rn_uzonal / grav * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
     116            ELSEWHERE 
     117               pssh(:,:) = - rn_uzonal / grav * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3   & 
    101118                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 ) 
    102119            END WHERE 
    … …  
    107124         pts(:,:,jpk,jp_sal) = 0. 
    108125         DO jk=1, jpkm1 
    109             pts(:,:,jk,jp_sal) = gphit(:,:) 
     126            WHERE( ABS(gphit) <= zjety ) 
     127!!$            WHERE( ABS(gphit) <= zjety*0.5 .AND. ABS(glamt) <= zjety*0.5 ) ! for a square of salt 
     128               pts(:,:,jk,jp_sal) = 35. 
     129            ELSEWHERE 
     130               pts(:,:,jk,jp_sal) = 30. 
     131            END WHERE                     
    110132         END DO 
    111133         ! velocities: 
    … …  
    132154            WHERE( ABS(gphit) <= zjety ) 
    133155               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   & 
    134                   &        * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
     156                  &        * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 ) 
    135157            ELSEWHERE 
    136158               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   & 
    137                   &        * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 ) 
     159                  &        * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 ) 
    138160            END WHERE 
    139161         END SELECT 
    … …  
    141163         pts(:,:,:,jp_tem) = 10._wp 
    142164         ! salinity:   
    143          pts(:,:,:,jp_sal) = 2. 
    144          DO jk=1, jpkm1 
    145             WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 2. + SIGN(1.,gphiv(:,:)) 
     165         pts(:,:,:,jp_sal) = 30. 
     166         DO jk=1, jpkm1 
     167            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 30. + SIGN(1.,gphiv(:,:)) 
    146168         END DO 
    147169         ! velocities: 
    … …  
    164186         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1) 
    165187         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 
     188            DO_2D( 0, 0, 0, 0 ) 
     189               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
     190            END_2D 
    171191            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    172192         END DO 
    … …  
    176196         ! salinity:   
    177197         DO jk=1, jpkm1 
    178             pts(:,:,jk,jp_sal) = gphit(:,:) 
     198            pts(:,:,jk,jp_sal) = pssh(:,:) 
    179199         END DO 
    180200         ! velocities: 
    … …  
    183203      CASE(4)    ! geostrophic zonal pulse 
    184204    
    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 
     205         DO_2D( 1, 1, 1, 1 ) 
     206            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
     207               zdu = rn_uzonal 
     208            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
     209               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
     210            ELSE 
     211               zdu = 0. 
     212            END IF 
     213            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
     214               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
     215               pu(ji,jj,:) = zdu 
     216               pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
     217            ELSE 
     218               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
     219               pu(ji,jj,:) = 0. 
     220               pts(ji,jj,:,jp_sal) = 1. 
     221            END IF 
     222         END_2D 
    205223          
    206224         ! temperature: 
    207225         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)         
    208226         pv(:,:,:) = 0. 
    209           
    210227          
    211228       CASE(5)    ! vortex 
    … …  
    213230         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 ) 
    214231         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic 
    215          zlambda = SQRT(2._wp)*rn_lambda       ! Horizontal scale in meters  
     232         zlambda = SQRT(2._wp)*rn_lambda*1.e3       ! Horizontal scale in meters  
    216233         zn2 = 3.e-3**2 
    217234         zH = 0.5_wp * 5000._wp 
    … …  
    220237         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 
    221238         ! 
    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 
     239         DO_2D( 1, 1, 1, 1 ) 
     240            zx = glamt(ji,jj) * 1.e3 
     241            zy = gphit(ji,jj) * 1.e3 
     242            ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
     243            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
     244            ! Sea level: 
     245            pssh(ji,jj) = 0. 
     246            DO jl=1,5 
     247               zdt = pssh(ji,jj) 
     248               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
     249               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     250               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
    247251            END DO 
    248          END DO 
     252            ! temperature: 
     253            DO jk=1,jpk 
     254               zdt =  pdept(ji,jj,jk)  
     255               zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
     256               IF (zdt < zH) THEN 
     257                  zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
     258                  zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     259               ENDIF 
     260               !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     261               pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     262            END DO 
     263         END_2D 
    249264         ! 
    250265         ! salinity:   
    … …  
    253268         ! velocities: 
    254269         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 
     270         DO_2D( 0, 0, 0, 0 ) 
     271            zx = glamu(ji,jj) * 1.e3 
     272            zy = gphiu(ji,jj) * 1.e3 
     273            DO jk=1, jpk 
     274               zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
     275               IF (zdu < zH) THEN 
     276                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
     277                  zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
     278                  pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
     279               ELSE 
     280                  pu(ji,jj,jk) = 0._wp 
     281               ENDIF 
    269282            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 
     283         END_2D 
     284         ! 
     285         DO_2D( 0, 0, 0, 0 ) 
     286            zx = glamv(ji,jj) * 1.e3 
     287            zy = gphiv(ji,jj) * 1.e3 
     288            DO jk=1, jpk 
     289               zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
     290               IF (zdv < zH) THEN 
     291                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
     292                  zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
     293                  pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
     294               ELSE 
     295                  pv(ji,jj,jk) = 0._wp 
     296               ENDIF 
    286297            END DO 
    287          END DO 
     298         END_2D 
    288299         !             
    289300      END SELECT 
    290  
     301       
    291302      IF (ln_sshnoise) THEN 
     303         CALL RANDOM_SEED() 
    292304         CALL RANDOM_NUMBER(zrandom) 
    293305         pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 ) 
    294306      END IF 
    295307      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. ) 
     308      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. ) 
     309      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 
    299310 
    300311   END SUBROUTINE usr_def_istate 
Note: See TracChangeset for help on using the changeset viewer.