Optionally rotate wind/stress from to computational grid (from true n…

…orth) (#473)

* Optionally rotate wind/stress composant to computational grid

* Change if-else statements to become more readable

* Typo correction

cicecore/cicedynB/general/ice_forcing.F90

@@ -122,6 +122,9 @@ module ice_forcing
                           ! 'hadgem_sst' or 'hadgem_sst_uvocn'
          ice_data_type, & ! 'default', 'box2001', 'boxslotcyl'
          precip_units     ! 'mm_per_month', 'mm_per_sec', 'mks','m_per_sec'
+
+      logical (kind=log_kind), public :: &
+         rotate_wind      ! rotate wind/stress to computational grid from true north directed
 
       character(char_len_long), public :: & 
          atm_data_dir , & ! top directory for atmospheric data
@@ -1628,11 +1631,10 @@ subroutine prepare_forcing (nx_block, ny_block, &
 
       if (calc_strair) then
 
-        do j = jlo, jhi
-        do i = ilo, ihi
-
-            wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
-
+        if (rotate_wind) then
+          do j = jlo, jhi
+          do i = ilo, ihi
+             wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
       !-----------------------------------------------------------------
       ! Rotate zonal/meridional vectors to local coordinates.
       ! Velocity comes in on T grid, but is oriented geographically ---
@@ -1644,30 +1646,38 @@ subroutine prepare_forcing (nx_block, ny_block, &
       ! atmo_boundary_layer, and are interpolated to the U grid later as 
       ! necessary.
       !-----------------------------------------------------------------
-           workx      = uatm(i,j) ! wind velocity, m/s
-           worky      = vatm(i,j)
-           uatm (i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
-                      + worky*sin(ANGLET(i,j))   ! note uatm, vatm, wind
-           vatm (i,j) = worky*cos(ANGLET(i,j)) & !  are on the T-grid here
-                      - workx*sin(ANGLET(i,j))
-
-        enddo                     ! i
-        enddo                     ! j
-
-      else  ! strax, stray, wind are read from files
-
-        do j = jlo, jhi
-        do i = ilo, ihi
-
-           workx      = strax(i,j) ! wind stress
-           worky      = stray(i,j)
-           strax(i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
-                      + worky*sin(ANGLET(i,j))   ! note strax, stray, wind
-           stray(i,j) = worky*cos(ANGLET(i,j)) & !  are on the T-grid here
-                      - workx*sin(ANGLET(i,j))
-
-        enddo                     ! i
-        enddo                     ! j
+             workx      = uatm(i,j) ! wind velocity, m/s
+             worky      = vatm(i,j)
+             uatm (i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
+                        + worky*sin(ANGLET(i,j))   ! note uatm, vatm, wind
+             vatm (i,j) = worky*cos(ANGLET(i,j)) & !  are on the T-grid here
+                        - workx*sin(ANGLET(i,j))
+          enddo                     ! i
+          enddo                     ! j
+        else ! not rotated
+          do j = jlo, jhi
+          do i = ilo, ihi
+             wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
+          enddo                     ! i
+          enddo                     ! j
+        endif ! rotated
+
+      else ! strax, stray, wind are read from files
+
+        if (rotate_wind) then
+          do j = jlo, jhi
+          do i = ilo, ihi
+             workx      = strax(i,j) ! wind stress
+             worky      = stray(i,j)
+             strax(i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
+                        + worky*sin(ANGLET(i,j))   ! note strax, stray, wind
+             stray(i,j) = worky*cos(ANGLET(i,j)) & !  are on the T-grid here
+                        - workx*sin(ANGLET(i,j))
+          enddo                     ! i
+          enddo                     ! j
+        else ! not rotated
+          ! wind (speed) is already read from file, so all is in place
+        endif ! rotated
 
       endif                   ! calc_strair
 
@@ -2050,11 +2060,11 @@ subroutine JRA55_gx1_files(yr)
       uwind_file = &
            trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_2005.nc'
       call file_year(uwind_file,yr)
-  if (my_task == master_task) then
+      if (my_task == master_task) then
          write (nu_diag,*) ' '
          write (nu_diag,*) 'Atmospheric data files:'
          write (nu_diag,*) trim(uwind_file)
-    endif
+      endif
       end subroutine JRA55_gx1_files
       subroutine JRA55_tx1_files(yr)
 !
@@ -2066,11 +2076,11 @@ subroutine JRA55_tx1_files(yr)
       uwind_file = &
            trim(atm_data_dir)//'/8XDAILY/JRA55_03hr_forcing_tx1_2005.nc'
       call file_year(uwind_file,yr)
-  if (my_task == master_task) then
+      if (my_task == master_task) then
          write (nu_diag,*) ' '
          write (nu_diag,*) 'Atmospheric data files:'
          write (nu_diag,*) trim(uwind_file)
-    endif
+      endif
       end subroutine JRA55_tx1_files
       subroutine JRA55_gx3_files(yr)
 !
@@ -2082,11 +2092,11 @@ subroutine JRA55_gx3_files(yr)
       uwind_file = &
            trim(atm_data_dir)//'/8XDAILY/JRA55_gx3_03hr_forcing_2005.nc'
       call file_year(uwind_file,yr)
-  if (my_task == master_task) then
+      if (my_task == master_task) then
          write (nu_diag,*) ' '
          write (nu_diag,*) 'Atmospheric data files:'
          write (nu_diag,*) trim(uwind_file)
-    endif
+      endif
       end subroutine JRA55_gx3_files
 !=======================================================================
 !
@@ -4471,7 +4481,7 @@ subroutine hycom_atm_data
          write (nu_diag,*) &
          'ERROR: CICE: Atm forcing not available at hcdate =',hcdate
          write (nu_diag,*) &
-         'ERROR: CICE: nyr, year_init, yday = ',nyr, year_init, yday
+         'ERROR: CICE: nyr, year_init, yday ,sec = ',nyr, year_init, yday, sec
          call abort_ice ('ERROR: CICE stopped')
       endif
 

cicecore/cicedynB/general/ice_init.F90

@@ -83,7 +83,7 @@ subroutine input_data
       use ice_flux_bgc, only: cpl_bgc
       use ice_forcing, only: &
           ycycle,          fyear_init,    dbug, &
-          atm_data_type,   atm_data_dir,  precip_units, &
+          atm_data_type,   atm_data_dir,  precip_units, rotate_wind, &
           atm_data_format, ocn_data_format, &
           bgc_data_type, &
           ocn_data_type, ocn_data_dir,    wave_spec_file,  &
@@ -206,15 +206,15 @@ subroutine input_data
       namelist /forcing_nml/ &
         formdrag,       atmbndy,         calc_strair,   calc_Tsfc,      &
         highfreq,       natmiter,        atmiter_conv,                  &
-        ustar_min,      emissivity,     &
+        ustar_min,      emissivity,                                     &
         fbot_xfer_type, update_ocn_f,    l_mpond_fresh, tfrz_option,    &
         oceanmixed_ice, restore_ice,     restore_ocn,   trestore,       &
         precip_units,   default_season,  wave_spec_type,nfreq,          &
         atm_data_type,  ocn_data_type,   bgc_data_type, fe_data_type,   &
         ice_data_type,  wave_spec_file,                                 &
         fyear_init,     ycycle,                                         &
         atm_data_dir,   ocn_data_dir,    bgc_data_dir,                  &
-        atm_data_format, ocn_data_format,                               &
+        atm_data_format, ocn_data_format, rotate_wind,                  &
         oceanmixed_file
 
       !-----------------------------------------------------------------
@@ -357,6 +357,7 @@ subroutine input_data
       atm_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
       atm_data_type   = 'default'
       atm_data_dir    = ' '
+      rotate_wind     = .true.    ! rotate wind/stress composants to computational grid orientation
       calc_strair     = .true.    ! calculate wind stress
       formdrag        = .false.   ! calculate form drag
       highfreq        = .false.   ! calculate high frequency RASM coupling
@@ -631,6 +632,7 @@ subroutine input_data
       call broadcast_scalar(atm_data_format,    master_task)
       call broadcast_scalar(atm_data_type,      master_task)
       call broadcast_scalar(atm_data_dir,       master_task)
+      call broadcast_scalar(rotate_wind,        master_task)
       call broadcast_scalar(calc_strair,        master_task)
       call broadcast_scalar(calc_Tsfc,          master_task)
       call broadcast_scalar(formdrag,           master_task)
@@ -1279,6 +1281,7 @@ subroutine input_data
          write(nu_diag,*) '--------------------------------'
          write(nu_diag,1012) ' calc_Tsfc        = ', calc_Tsfc,' calculate surface temperature as part of thermo'
          write(nu_diag,1012) ' calc_strair      = ', calc_strair,' calculate wind stress and speed'
+         write(nu_diag,1012) ' rotate_wind      = ', rotate_wind,' rotate wind/stress to computational grid'
          write(nu_diag,1012) ' formdrag         = ', formdrag,' use form drag parameterization'
          if (trim(atmbndy) == 'constant') then
             tmpstr2 = ': stability-based boundary layer'

configuration/scripts/ice_in

@@ -163,6 +163,7 @@
 &forcing_nml
     formdrag        = .false.
     atmbndy         = 'default'
+    rotate_wind     = .true.
     calc_strair     = .true.
     calc_Tsfc       = .true.
     highfreq        = .false.

doc/source/cice_index.rst

@@ -91,6 +91,7 @@ either Celsius or Kelvin units).
    "bTiz", "temperature of ice layers on bio grid", ""
    "**C**", "", ""
    "c<n>", "real(\ :math:`n`)", ""
+   "rotate_wind", ":math:`\bullet` if true, rotate wind/stress components to computational grid", "T"
    "calc_strair", ":math:`\bullet` if true, calculate wind stress", "T"
    "calc_Tsfc", ":math:`\bullet` if true, calculate surface temperature", "T"
    "Cdn_atm", "atmospheric drag coefficient", ""