Changeset 15813

Timestamp:

2023-04-05T13:59:14+02:00 (21 months ago)

Author:

clem

Message:

SI3: change heat budget to avoid supercooling. The rest: cosmetics

Location:

NEMO/releases/r4.0/r4.0-HEAD/src

Files:

• ICE/icethd.F90 (modified) (1 diff)
• OCE/DYN/dynnxt.F90 (modified) (1 diff)
• OCE/SBC/sbcblk.F90 (modified) (2 diffs)
• OCE/SBC/sbcflx.F90 (modified) (5 diffs)

NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90

@@ -183 +183 @@
             !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2)
             !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2)
-            IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN
-               ! upper bound for fhld: fhld should be equal to zqld
-               !                        but we have to make sure that this heat will not make the sst drop below the freezing point
-               !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos
-               !                        The following formulation is ok for both normal conditions and supercooling
-               fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
-                  &                                 - qsb_ice_bot(ji,jj) )
-               qlead(ji,jj) = 0._wp
-            ELSE
+            IF( ( zqld - zqfr ) < 0._wp .OR. at_i(ji,jj) < epsi10 ) THEN
                fhld (ji,jj) = 0._wp
                ! upper bound for qlead: qlead should be equal to zqld
                !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point.
-               !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb)
-               !                        and freezing point is reached if zqfr = zqld - qsb*a/dt
-               !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr
+               !                        The energy for this cooling down is zqfr and freezing point is reached if zqfr = zqld
+               !                        so the max heat that can be pulled out of the ocean is zqld - zqfr
+               !                        The following formulation is ok for both normal conditions and supercooling            
+               qlead(ji,jj) = MIN( 0._wp , zqld - zqfr )
+            ELSE
+               ! upper bound for fhld: fhld should be equal to zqld
+               !                        but we have to make sure that this heat will not make the sst drop below the freezing point
+               !                        so the max heat that can be pulled out of the ocean is zqld - zqfr_pos
                !                        The following formulation is ok for both normal conditions and supercooling
-               qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
+               fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
+               qlead(ji,jj) = 0._wp
             ENDIF
             !

NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynnxt.F90

@@ -356 +356 @@
       ENDIF
       !
+      !!clem
+      !!CALL lbc_lnk_multi( 'dynnxt' , ub, 'U', -1., vn, 'V', -1.)
+
+      !
       IF ( iom_use("utau") ) THEN
          IF ( ln_drgice_imp.OR.ln_isfcav ) THEN

NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcblk.F90

@@ -422 +422 @@
       END DO
 #endif
+      ! mask the wind in case it is not masked in the input file otherwise we end up with twice the stress at the coast
+      ! (see calculation of utau/vtau below)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vect. opt.
-            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
+            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) * tmask(ji,jj,1)
+            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) * tmask(ji,jj,1)
          END DO
       END DO
@@ -431 +433 @@
       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
-         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
+         &             + zwnd_j(:,:) * zwnd_j(:,:)  )
 
       ! ----------------------------------------------------------------------------- !

NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcflx.F90

@@ -35 +35 @@
    INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
    INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
- !!INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux flux
-   INTEGER , PARAMETER ::   jpfld   = 5 !! 6 ! maximum number of files to read 
+!!$   INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux
+!!$   INTEGER , PARAMETER ::   jp_sithic = 7    ! index of sea ice thickness
+!!$   INTEGER , PARAMETER ::   jp_siconc = 8    ! index of sea ice fraction
+!!$   INTEGER , PARAMETER ::   jp_hfrnf  = 9    ! index of rnf heat flux
+!!$   INTEGER , PARAMETER ::   jp_hfisf  = 10   ! index of iceshelf heat flux
+!!$   INTEGER , PARAMETER ::   jp_fwisf  = 11   ! index of iceshelf freshwater flux
+!!$   INTEGER , PARAMETER ::   jp_fwrnf  = 12   ! index of runoff freshwater flux
+   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read 
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
@@ -74 +80 @@
       !!              - qsr         solar heat flux
       !!              - emp         upward mass flux (evap. - precip.)
-      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
-      !!                            if ice
+      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero if ice
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
@@ -89 +94 @@
       CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
-      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx ! informations about the fields to be read
-      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx
+      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  !!, & ! informations about the fields to be read
+!!$         &             sn_sfx, sn_sithic, sn_siconc, sn_hfisf, sn_hfrnf, sn_fwisf, sn_fwrnf
+      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, &
+!!$         &                         sn_sfx, sn_sithic, sn_siconc, sn_hfisf, sn_hfrnf, sn_fwisf, sn_fwrnf
       !!---------------------------------------------------------------------
       !
@@ -111 +118 @@
          slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
          slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
-         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx
+         slf_i(jp_emp ) = sn_emp
+!!$         slf_i(jp_sfx  ) = sn_sfx
+!!$         slf_i(jp_sithic) = sn_sithic
+!!$         slf_i(jp_siconc) = sn_siconc
+!!$         slf_i(jp_hfisf) = sn_hfisf    ;   slf_i(jp_hfrnf) = sn_hfrnf
+!!$         slf_i(jp_fwisf) = sn_fwisf    ;   slf_i(jp_fwrnf) = sn_fwrnf
          !
          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
@@ -146 +158 @@
                qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
                emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
-               !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1) 
+!!$               sfx (ji,jj) = sf(jp_sfx   )%fnow(ji,jj,1)                             * tmask(ji,jj,1) 
+!!$               !! => if the following is used, then one needs to change tke routine + allocate hm_i in sbc_oce
+!!$               hm_i(ji,jj) = sf(jp_sithic)%fnow(ji,jj,1)                             * tmask(ji,jj,1) 
+!!$               fr_i(ji,jj) = sf(jp_siconc)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
+!!$               !! => if the following is used, then one needs to change rnf and isf routines + allocate the arrays
+!!$               hfisf(ji,jj) = sf(jp_hfisf)%fnow(ji,jj,1)                             * ssmask(ji,jj)
+!!$               fwisf(ji,jj) = sf(jp_fwisf)%fnow(ji,jj,1)                             * ssmask(ji,jj)
+!!$               hfrnf(ji,jj) = sf(jp_hfrnf)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
+!!$               fwrnf(ji,jj) = sf(jp_fwrnf)%fnow(ji,jj,1)                             * tmask(ji,jj,1)
             END DO
          END DO
-         !                                                        ! add to qns the heat due to e-p
-         !clem: I do not think it is needed
-         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
          !
          ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x) 
          CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
-            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
+            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! &
+!!$            &                           sfx, 'T', 1._wp, hm_i,'T', 1._wp, fr_i,'T', 1._wp, &
+!!$            &                          hfisf, 'T', 1._wp, fwisf, 'T', 1._wp, hfrnf, 'T', 1._wp, fwrnf, 'T', 1._wp )
          !
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)