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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/SBC/cyclone.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

Location:
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • 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 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 

  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/SBC/cyclone.F90

    r10068 r13463  
    3737 
    3838   !! * Substitutions 
    39 #  include "vectopt_loop_substitute.h90" 
     39#  include "do_loop_substitute.h90" 
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    137137            zhemi = SIGN( 1. , zrlat ) 
    138138            zinfl = 15.* rad                             ! clim inflow angle in Tropical Cyclones 
    139          IF ( vortex == 0 ) THEN 
     139         IF( vortex == 0 ) THEN 
    140140 
    141141            ! Vortex Holland reconstruct wind at each lon-lat position 
    … …  
    147147            zb = 2. 
    148148 
    149             DO jj = 1, jpj 
    150                DO ji = 1, jpi 
    151  
    152                   ! calc distance between TC center and any point following great circle 
    153                   ! source : http://www.movable-type.co.uk/scripts/latlong.html 
    154                   zzrglam = rad * glamt(ji,jj) - zrlon 
    155                   zzrgphi = rad * gphit(ji,jj) 
    156                   zdist = ra * ACOS(  SIN( zrlat ) * SIN( zzrgphi )   & 
    157                      &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 
    158  
    159                  IF (zdist < zrout2) THEN ! calculation of wind only to a given max radius 
    160                   ! shape of the wind profile 
    161                   zztmp = ( zrmw / ( zdist + 1.e-12 ) )**zb 
    162                   zztmp =  zvmax * SQRT( zztmp * EXP(1. - zztmp) )     
    163  
    164                   IF (zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
    165                      zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 
    166                   ENDIF 
    167  
    168                   ! !!! KILL EQ WINDS 
    169                   ! IF (SIGN( 1. , zrlat ) /= zhemi) THEN 
    170                   !    zztmp = 0.                              ! winds in other hemisphere 
    171                   !    IF (ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
    172                   ! ENDIF 
    173                   ! IF (ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
    174                   !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) )  
    175                   !    !linear to zero between 10 and 5 
    176                   ! ENDIF 
    177                   ! !!! / KILL EQ 
    178  
    179                   IF (ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
    180  
    181                   zwnd_t =   COS( zinfl ) * zztmp     
    182                   zwnd_r = - SIN( zinfl ) * zztmp 
    183  
    184                   ! Project radial-tangential components on zonal-meridional components 
    185                   ! ------------------------------------------------------------------- 
    186                    
    187                   ! ztheta = azimuthal angle of the great circle between two points 
    188                   zztmp = COS( zrlat ) * SIN( zzrgphi ) & 
    189                      &  - SIN( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) 
    190                   ztheta = ATAN2(        COS( zzrgphi ) * SIN( zzrglam ) , zztmp ) 
    191  
    192                   zwnd_x(ji,jj) = zwnd_x(ji,jj) - zhemi * COS(ztheta)*zwnd_t + SIN(ztheta)*zwnd_r 
    193                   zwnd_y(ji,jj) = zwnd_y(ji,jj) + zhemi * SIN(ztheta)*zwnd_t + COS(ztheta)*zwnd_r 
    194                  ENDIF 
    195                END DO 
    196             END DO 
     149            DO_2D( 1, 1, 1, 1 ) 
     150 
     151               ! calc distance between TC center and any point following great circle 
     152               ! source : http://www.movable-type.co.uk/scripts/latlong.html 
     153               zzrglam = rad * glamt(ji,jj) - zrlon 
     154               zzrgphi = rad * gphit(ji,jj) 
     155               zdist = ra * ACOS(  SIN( zrlat ) * SIN( zzrgphi )   & 
     156                  &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 
     157 
     158              IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 
     159               ! shape of the wind profile 
     160               zztmp = ( zrmw / ( zdist + 1.e-12 ) )**zb 
     161               zztmp =  zvmax * SQRT( zztmp * EXP(1. - zztmp) )     
     162 
     163               IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
     164                  zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 
     165               ENDIF 
     166 
     167               ! !!! KILL EQ WINDS 
     168               ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 
     169               !    zztmp = 0.                              ! winds in other hemisphere 
     170               !    IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
     171               ! ENDIF 
     172               ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
     173               !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) )  
     174               !    !linear to zero between 10 and 5 
     175               ! ENDIF 
     176               ! !!! / KILL EQ 
     177 
     178               IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
     179 
     180               zwnd_t =   COS( zinfl ) * zztmp     
     181               zwnd_r = - SIN( zinfl ) * zztmp 
     182 
     183               ! Project radial-tangential components on zonal-meridional components 
     184               ! ------------------------------------------------------------------- 
     185                
     186               ! ztheta = azimuthal angle of the great circle between two points 
     187               zztmp = COS( zrlat ) * SIN( zzrgphi ) & 
     188                  &  - SIN( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) 
     189               ztheta = ATAN2(        COS( zzrgphi ) * SIN( zzrglam ) , zztmp ) 
     190 
     191               zwnd_x(ji,jj) = zwnd_x(ji,jj) - zhemi * COS(ztheta)*zwnd_t + SIN(ztheta)*zwnd_r 
     192               zwnd_y(ji,jj) = zwnd_y(ji,jj) + zhemi * SIN(ztheta)*zwnd_t + COS(ztheta)*zwnd_r 
     193              ENDIF 
     194            END_2D 
    197195          
    198          ELSE IF ( vortex == 1 ) THEN 
     196         ELSE IF( vortex == 1 ) THEN 
    199197 
    200198            ! Vortex Willoughby reconstruct wind at each lon-lat position 
    … …  
    206204            zn   =   2.1340 + 0.0077*zvmax - 0.4522*LOG(zrmw/1000.) - 0.0038*ABS( ztct(jtc,jp_lat) )             
    207205            zA   =   0.5913 + 0.0029*zvmax - 0.1361*LOG(zrmw/1000.) - 0.0042*ABS( ztct(jtc,jp_lat) )   
    208             IF (zA < 0) THEN  
     206            IF(zA < 0) THEN  
    209207               zA=0 
    210208            ENDIF            
    211209         
    212             DO jj = 1, jpj 
    213                DO ji = 1, jpi 
    214                                    
    215                   zzrglam = rad * glamt(ji,jj) - zrlon 
    216                   zzrgphi = rad * gphit(ji,jj) 
    217                   zdist = ra * ACOS(  SIN( zrlat ) * SIN( zzrgphi )   & 
    218                      &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 
    219  
    220                  IF (zdist < zrout2) THEN ! calculation of wind only to a given max radius 
     210            DO_2D( 1, 1, 1, 1 ) 
     211                                
     212               zzrglam = rad * glamt(ji,jj) - zrlon 
     213               zzrgphi = rad * gphit(ji,jj) 
     214               zdist = ra * ACOS(  SIN( zrlat ) * SIN( zzrgphi )   & 
     215                  &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 
     216 
     217              IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 
     218             
     219               ! shape of the wind profile                      
     220               IF(zdist <= zrmw) THEN     ! inside the Radius of Maximum Wind 
     221                  zztmp  = zvmax * (zdist/zrmw)**zn 
     222               ELSE  
     223                  zztmp  = zvmax * ( (1-zA) * EXP(- (zdist-zrmw)/zXX1 ) + zA * EXP(- (zdist-zrmw)/zXX2 ) ) 
     224               ENDIF 
     225 
     226               IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
     227                  zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 
     228               ENDIF 
     229 
     230               ! !!! KILL EQ WINDS 
     231               ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 
     232               !    zztmp = 0.                              ! winds in other hemisphere 
     233               !    IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
     234               ! ENDIF 
     235               ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
     236               !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) )  
     237               !    !linear to zero between 10 and 5 
     238               ! ENDIF 
     239               ! !!! / KILL EQ 
     240 
     241               IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
     242 
     243               zwnd_t =   COS( zinfl ) * zztmp     
     244               zwnd_r = - SIN( zinfl ) * zztmp 
     245 
     246               ! Project radial-tangential components on zonal-meridional components 
     247               ! ------------------------------------------------------------------- 
    221248                
    222                   ! shape of the wind profile                      
    223                   IF (zdist <= zrmw) THEN     ! inside the Radius of Maximum Wind 
    224                      zztmp  = zvmax * (zdist/zrmw)**zn 
    225                   ELSE  
    226                      zztmp  = zvmax * ( (1-zA) * EXP(- (zdist-zrmw)/zXX1 ) + zA * EXP(- (zdist-zrmw)/zXX2 ) ) 
    227                   ENDIF 
    228  
    229                   IF (zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
    230                      zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 
    231                   ENDIF 
    232  
    233                   ! !!! KILL EQ WINDS 
    234                   ! IF (SIGN( 1. , zrlat ) /= zhemi) THEN 
    235                   !    zztmp = 0.                              ! winds in other hemisphere 
    236                   !    IF (ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
    237                   ! ENDIF 
    238                   ! IF (ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
    239                   !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) )  
    240                   !    !linear to zero between 10 and 5 
    241                   ! ENDIF 
    242                   ! !!! / KILL EQ 
    243  
    244                   IF (ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
    245  
    246                   zwnd_t =   COS( zinfl ) * zztmp     
    247                   zwnd_r = - SIN( zinfl ) * zztmp 
    248  
    249                   ! Project radial-tangential components on zonal-meridional components 
    250                   ! ------------------------------------------------------------------- 
    251                    
    252                   ! ztheta = azimuthal angle of the great circle between two points 
    253                   zztmp = COS( zrlat ) * SIN( zzrgphi ) & 
    254                      &  - SIN( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) 
    255                   ztheta = ATAN2(        COS( zzrgphi ) * SIN( zzrglam ) , zztmp ) 
    256  
    257                   zwnd_x(ji,jj) = zwnd_x(ji,jj) - zhemi * COS(ztheta)*zwnd_t + SIN(ztheta)*zwnd_r 
    258                   zwnd_y(ji,jj) = zwnd_y(ji,jj) + zhemi * SIN(ztheta)*zwnd_t + COS(ztheta)*zwnd_r 
    259                    
    260                  ENDIF 
    261                END DO 
    262             END DO 
     249               ! ztheta = azimuthal angle of the great circle between two points 
     250               zztmp = COS( zrlat ) * SIN( zzrgphi ) & 
     251                  &  - SIN( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) 
     252               ztheta = ATAN2(        COS( zzrgphi ) * SIN( zzrglam ) , zztmp ) 
     253 
     254               zwnd_x(ji,jj) = zwnd_x(ji,jj) - zhemi * COS(ztheta)*zwnd_t + SIN(ztheta)*zwnd_r 
     255               zwnd_y(ji,jj) = zwnd_y(ji,jj) + zhemi * SIN(ztheta)*zwnd_t + COS(ztheta)*zwnd_r 
     256                
     257              ENDIF 
     258            END_2D 
    263259         ENDIF                                         ! / vortex Holland or Wiloughby 
    264260         ENDIF                                           ! / cyclone is defined in this slot ? yes--> begin 
Note: See TracChangeset for help on using the changeset viewer.