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.
cyclone.F90 in NEMO/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/cyclone.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 13.8 KB
Line 
1MODULE cyclone
2   !!======================================================================
3   !!                       ***  MODULE  cyclone  ***
4   !! add the Tropical Cyclones along tracks to the surface wind forcing
5   !!                 
6   !!======================================================================
7   !! History : 3.3  ! 2010-05  (E Vincent, G Madec, S Masson)  Original code
8   !!----------------------------------------------------------------------
9
10#if defined key_cyclone
11   !!----------------------------------------------------------------------
12   !!  'key_cyclone' : key option add Tropical Cyclones in the wind forcing
13   !!----------------------------------------------------------------------
14   !!   wnd_cyc      : 1 module subroutine
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and active tracers
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE dom_oce         ! ocean space domain variables
19   USE phycst          ! physical constant
20   USE fldread         ! read input fields
21   USE in_out_manager  ! I/O manager
22   USE geo2ocean       ! tools for projection on ORCA grid
23   USE lib_mpp       
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   wnd_cyc   ! routine called in sbcblk.F90 module
29
30   INTEGER , PARAMETER ::   jp_is1  = 1   ! index of presence 1 or absence 0 of a TC record
31   INTEGER , PARAMETER ::   jp_lon  = 2   ! index of longitude for present TCs
32   INTEGER , PARAMETER ::   jp_lat  = 3   ! index of latitude  for present TCs
33   INTEGER , PARAMETER ::   jp_vmax = 4   ! index of max wind  for present TCs
34   INTEGER , PARAMETER ::   jp_pres = 5   ! index of eye-pres  for present TCs
35
36   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
37
38   !! * Substitutions
39#  include "do_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id$
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE wnd_cyc( kt, pwnd_i, pwnd_j )
49      !!----------------------------------------------------------------------
50      !!                    ***  ROUTINE wnd_cyc  ***
51      !!
52      !! ** Purpose :  Add cyclone winds on the ORCA grid
53      !!
54      !! ** Action  : - open TC data, find TCs for the current timestep
55      !!              - for each potential TC, add the winds on the grid
56      !!----------------------------------------------------------------------
57      INTEGER , INTENT(in)                      ::   kt       ! time step index
58      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   pwnd_i   ! wind speed i-components at T-point ORCA direction
59      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   pwnd_j   ! wind speed j-components at T-point ORCA direction
60      !
61      !!
62      INTEGER  ::   ji, jj , jtc        ! loop arguments
63      INTEGER  ::   ierror              ! loop arguments
64      INTEGER  ::   vortex=1            ! vortex shape to be used: 0=Holland 1=Willoughby
65      REAL(wp) ::   zrout1=1.5e6        ! distance from center where we begin to kill vortex (m)
66      REAL(wp) ::   zrout2=2.5e6        ! distance from center where we bring vortex to zero (m)
67      REAL(wp) ::   zb                  ! power in Holland vortex shape
68      REAL(wp) ::   zA                  ! shape parameter in Willoughby vortex : A transtion between first and second outter exp
69      REAL(wp) ::   zn                  ! shape parameter in Willoughby vortex : n power law in the eye
70      REAL(wp) ::   zXX1                ! shape parameter in Willoughby vortex : decay length second outter exponential
71      REAL(wp) ::   zXX2                ! shape parameter in Willoughby vortex : decay length first  outter exponential
72      REAL(wp) ::   zztmp               ! temporary
73      REAL(wp) ::   zzrglam, zzrgphi    ! temporary
74      REAL(wp) ::   ztheta              ! azimuthal angle
75      REAL(wp) ::   zdist               ! dist to the TC center
76      REAL(wp) ::   zhemi               ! 1 for NH ;  -1 for SH
77      REAL(wp) ::   zinfl               ! clim inflow angle in TCs
78      REAL(wp) ::   zrmw                ! mean radius of Max wind of a tropical cyclone (Willoughby 2004) [m]
79      REAL(wp) ::   zwnd_r, zwnd_t      ! radial and tangential components of the wind
80      REAL(wp) ::   zvmax               ! timestep interpolated vmax
81      REAL(wp) ::   zrlon, zrlat        ! temporary
82      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_x, zwnd_y   ! zonal and meridional components of the wind
83      REAL(wp), DIMENSION(14,5)    ::   ztct                ! tropical cyclone track data at kt
84      !
85      CHARACTER(len=100) ::  cn_dir            ! Root directory for location of files
86      TYPE(FLD_N), DIMENSION(1) ::   slf_i     ! array of namelist informations on the TC position
87      TYPE(FLD_N) ::   sn_tc                   ! informations about the fields to be read
88      !!--------------------------------------------------------------------
89
90      !                                         ! ====================== !
91      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
92         !                                      ! ====================== !
93         ! set file information (default values)
94         cn_dir = './'       ! directory in which the model is executed
95         !
96         ! (NB: frequency positive => hours, negative => months)
97         !          !    file     ! frequency !  variable ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! land/sea mask !
98         !          !    name     !  (hours)  !   name    !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! filename      !
99         sn_tc = FLD_N( 'tc_track',     6     ,  'tc'     ,  .true.    , .false. ,   'yearly'  , ''       , ''         , ''            )
100         !
101         !  Namelist is read in namsbc_blk
102         ! set sf structure
103         ALLOCATE( sf(1), STAT=ierror )
104         IF( ierror > 0 ) THEN
105            CALL ctl_stop( 'wnd_cyc: unable to allocate sf structure' )   ;   RETURN
106         ENDIF
107         ALLOCATE( sf(1)%fnow(14,5,1) )
108         ALLOCATE( sf(1)%fdta(14,5,1,2) )
109         slf_i(1) = sn_tc
110         !
111         ! fill sf with slf_i and control print
112         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_tc', 'tropical cyclone track', 'namsbc_tc' )
113         !
114      ENDIF
115
116
117      !       Interpolation of lon lat vmax... at the current timestep
118      !       ***************************************************************
119
120      CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step
121
122      ztct(:,:) = sf(1)%fnow(:,:,1)
123
124      !       Add TC wind on the grid
125      !       ***************************************************************
126
127      zwnd_x(:,:) = 0.e0 
128      zwnd_y(:,:) = 0.e0 
129     
130      DO jtc = 1, 14
131         !
132         IF( ztct(jtc,jp_is1) == 1 ) THEN                ! cyclone is defined in this slot ? yes--> begin
133
134            zvmax =       ztct(jtc,jp_vmax)
135            zrlon = rad * ztct(jtc,jp_lon )
136            zrlat = rad * ztct(jtc,jp_lat )
137            zhemi = SIGN( 1. , zrlat )
138            zinfl = 15.* rad                             ! clim inflow angle in Tropical Cyclones
139         IF( vortex == 0 ) THEN
140
141            ! Vortex Holland reconstruct wind at each lon-lat position
142            ! ********************************************************
143            zrmw = 51.6 * EXP( -0.0223*zvmax + 0.0281* ABS( ztct(jtc,jp_lat) ) ) * 1000.
144            ! climatological ZRMW of cyclones as a function of wind and latitude (Willoughby 2004)             
145            ! zb = 1.0036 + 0.0173 * zvmax - 0.0313 * LOG(zrmw/1000.) + 0.0087 * ABS( ztct(jtc,jp_lat) )
146            ! fitted B parameter (Willoughby 2004)
147            zb = 2.
148
149            DO_2D_11_11
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
195         
196         ELSE IF( vortex == 1 ) THEN
197
198            ! Vortex Willoughby reconstruct wind at each lon-lat position
199            ! ***********************************************************
200            zrmw = 46.4 * EXP( -0.0155*zvmax + 0.0169* ABS( ztct(jtc,jp_lat) ) )*1000.
201            ! climatological ZRMW of cyclones as a function of wind and latitude (Willoughby 2006)
202            zXX2 = 25.*1000.                                              ! 25km fixed "near-eye" exponential decay
203            zXX1 = ( 287.6  - 1.942 *zvmax + 7.799 *LOG(zrmw/1000.) + 1.819 *ABS( ztct(jtc,jp_lat) ) )*1000.   
204            zn   =   2.1340 + 0.0077*zvmax - 0.4522*LOG(zrmw/1000.) - 0.0038*ABS( ztct(jtc,jp_lat) )           
205            zA   =   0.5913 + 0.0029*zvmax - 0.1361*LOG(zrmw/1000.) - 0.0042*ABS( ztct(jtc,jp_lat) ) 
206            IF(zA < 0) THEN
207               zA=0
208            ENDIF           
209       
210            DO_2D_11_11
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               ! -------------------------------------------------------------------
248               
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
259         ENDIF                                         ! / vortex Holland or Wiloughby
260         ENDIF                                           ! / cyclone is defined in this slot ? yes--> begin
261      END DO ! / end simultaneous cyclones loop
262
263      CALL rot_rep ( zwnd_x, zwnd_y, 'T', 'en->i', pwnd_i ) !rotation of components on ORCA grid
264      CALL rot_rep ( zwnd_x, zwnd_y, 'T', 'en->j', pwnd_j ) !rotation of components on ORCA grid
265
266   END SUBROUTINE wnd_cyc
267
268#endif
269
270   !!======================================================================
271END MODULE cyclone
Note: See TracBrowser for help on using the repository browser.