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 13886 for NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icectl.F90 – NEMO

Ignore:
Timestamp:
2020-11-26T15:24:38+01:00 (3 years ago)
Author:
rlod
Message:

phasing with trunk at revision r13787

Location:
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source

    • Property svn:externals

      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13507        sette 
         10^/utils/CI/sette@13559        sette 

  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icectl.F90

    r13472 r13886  
    4343   PUBLIC   ice_prt 
    4444   PUBLIC   ice_prt3D 
     45   PUBLIC   ice_drift_wri 
     46   PUBLIC   ice_drift_init 
    4547 
    4648   ! thresold rates for conservation 
    … …  
    4951   REAL(wp), PARAMETER ::   zchk_s   = 2.5e-6   ! g/m2/s  <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 
    5052   REAL(wp), PARAMETER ::   zchk_t   = 7.5e-2   ! W/m2    <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 
     53 
     54   ! for drift outputs 
     55   CHARACTER(LEN=50)   ::   clname="icedrift_diagnostics.ascii"   ! ascii filename 
     56   INTEGER             ::   numicedrift                           ! outfile unit 
     57   REAL(wp)            ::   rdiag_icemass, rdiag_icesalt, rdiag_iceheat  
     58   REAL(wp)            ::   rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat  
    5159    
    5260   !! * Substitutions 
    … …  
    132140 
    133141         ! -- advection scheme is conservative? -- ! 
    134          zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather) 
    135          zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 (only for Prather) 
     142         zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 
     143         zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
    136144 
    137145         ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    … …  
    156164               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
    157165            ! check if advection scheme is conservative 
    158             !    only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 
    159             !    so the formulation for conservation is different (and not coded)  
    160             !    it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 
    161             !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    162             !   &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 
     166            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     167               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
     168            IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     169               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rdt_ice 
    163170         ENDIF 
    164171         ! 
    … …  
    186193      ! water flux 
    187194      ! -- mass diag -- ! 
    188       zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 
     195      zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
     196         &                              + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 
    189197 
    190198      ! -- salt diag -- ! 
    191       zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) 
     199      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 
    192200 
    193201      ! -- heat diag -- ! 
    194       ! clem: not the good formulation 
    195       !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
    196       !!   &                              ) * e1e2t ) 
     202      zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     203      ! equivalent to this: 
     204      !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 
     205      !!   &                                          - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 
     206      !!   &                                          ) * e1e2t ) 
    197207 
    198208      ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    … …  
    204214         IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
    205215            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rDt_ice 
    206          !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
     216         IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 
     217            &                   WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice 
    207218      ENDIF 
    208219      ! 
    … …  
    725736       
    726737   END SUBROUTINE ice_prt3D 
     738 
     739 
     740   SUBROUTINE ice_drift_wri( kt ) 
     741      !!------------------------------------------------------------------- 
     742      !!                     ***  ROUTINE ice_drift_wri *** 
     743      !! 
     744      !! ** Purpose : conservation of mass, salt and heat 
     745      !!              write the drift in a ascii file at each time step 
     746      !!              and the total run drifts 
     747      !!------------------------------------------------------------------- 
     748      INTEGER, INTENT(in) ::   kt   ! ice time-step index 
     749      ! 
     750      INTEGER  ::   ji, jj 
     751      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 
     752      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 
     753      !!------------------------------------------------------------------- 
     754      ! 
     755      IF( kt == nit000 .AND. lwp ) THEN 
     756         WRITE(numout,*) 
     757         WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 
     758         WRITE(numout,*) '~~~~~~~~~~~~~' 
     759      ENDIF 
     760      ! 
     761      ! 2D budgets (must be close to 0) 
     762      IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 
     763         DO_2D( 1, 1, 1, 1 ) 
     764            zdiag_mass2D(ji,jj) =   wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 
     765               &                  + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 
     766            zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 
     767            zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 
     768         END_2D 
     769         ! 
     770         ! write outputs 
     771         CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 
     772         CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 
     773         CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 
     774      ENDIF 
     775 
     776      ! -- mass diag -- ! 
     777      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
     778         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 
     779      zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 
     780 
     781      ! -- salt diag -- ! 
     782      zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 
     783      zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 
     784 
     785      ! -- heat diag -- ! 
     786      zdiag_heat     = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     787      zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 
     788 
     789      !                    ! write out to file 
     790      IF( lwp ) THEN 
     791         ! check global drift (must be close to 0) 
     792         WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift     [kg]', zdiag_mass 
     793         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift     [kg]', zdiag_salt 
     794         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift     [W] ', zdiag_heat 
     795         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
     796         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'mass drift adv [kg]', zdiag_adv_mass 
     797         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift adv [kg]', zdiag_adv_salt 
     798         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift adv [W] ', zdiag_adv_heat 
     799      ENDIF 
     800      !                    ! drifts 
     801      rdiag_icemass = rdiag_icemass + zdiag_mass 
     802      rdiag_icesalt = rdiag_icesalt + zdiag_salt 
     803      rdiag_iceheat = rdiag_iceheat + zdiag_heat 
     804      rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 
     805      rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 
     806      rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 
     807      ! 
     808      !                    ! output drifts and close ascii file 
     809      IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 
     810         ! to ascii file 
     811         WRITE(numicedrift,*) '******************************************' 
     812         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift     [kg]', rdiag_icemass 
     813         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 
     814         WRITE(numicedrift,*) '******************************************' 
     815         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift     [kg]', rdiag_icesalt 
     816         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 
     817         WRITE(numicedrift,*) '******************************************' 
     818         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift     [W] ', rdiag_iceheat 
     819         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 
     820         CLOSE( numicedrift ) 
     821         ! 
     822         ! to ocean output 
     823         WRITE(numout,*) 
     824         WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 
     825         WRITE(numout,*) '~~~~~~~~~~~~~' 
     826         ! check global drift (must be close to 0) 
     827         WRITE(numout,*) '   sea-ice mass drift     [kg] = ', rdiag_icemass 
     828         WRITE(numout,*) '   sea-ice salt drift     [kg] = ', rdiag_icesalt 
     829         WRITE(numout,*) '   sea-ice heat drift     [W]  = ', rdiag_iceheat 
     830         ! check drift from advection scheme (can be /=0 with bdy but not sure why) 
     831         WRITE(numout,*) '   sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 
     832         WRITE(numout,*) '   sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 
     833         WRITE(numout,*) '   sea-ice heat drift adv [W]  = ', rdiag_adv_iceheat 
     834      ENDIF 
     835      ! 
     836   END SUBROUTINE ice_drift_wri 
     837 
     838   SUBROUTINE ice_drift_init 
     839      !!---------------------------------------------------------------------- 
     840      !!                  ***  ROUTINE ice_drift_init  *** 
     841      !!                    
     842      !! ** Purpose :   create output file, initialise arrays 
     843      !!---------------------------------------------------------------------- 
     844      ! 
     845      IF( .NOT.ln_icediachk ) RETURN ! exit 
     846      ! 
     847      IF(lwp) THEN 
     848         WRITE(numout,*) 
     849         WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 
     850         WRITE(numout,*) '~~~~~~~~~~~~~' 
     851         WRITE(numout,*) 
     852         ! 
     853         ! create output ascii file 
     854         CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 
     855         WRITE(numicedrift,*) 'Timestep  Drifts' 
     856         WRITE(numicedrift,*) '******************************************' 
     857      ENDIF 
     858      ! 
     859      rdiag_icemass = 0._wp 
     860      rdiag_icesalt = 0._wp 
     861      rdiag_iceheat = 0._wp 
     862      rdiag_adv_icemass = 0._wp 
     863      rdiag_adv_icesalt = 0._wp 
     864      rdiag_adv_iceheat = 0._wp 
     865      ! 
     866   END SUBROUTINE ice_drift_init 
    727867       
    728868#else 
Note: See TracChangeset for help on using the changeset viewer.