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 10249 for branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 – NEMO

Ignore:
Timestamp:
2018-10-29T13:03:40+01:00 (5 years ago)
Author:
kingr
Message:

Merged AMM15_v3_6_STABLE_package_collate@10237

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r10248 r10249  
    4040#endif 
    4141   USE sbc_oce          ! Surface boundary condition variables. 
     42   USE zdfmxl, ONLY :  &   
     43   &  hmld_tref,       &    
     44#if defined key_karaml 
     45   &  hmld_kara,       & 
     46   &  ln_kara,         & 
     47#endif    
     48   &  hmld,            &  
     49   &  hmlp,            & 
     50   &  hmlpt 
     51#if defined key_bdy  
     52   USE bdy_oce, ONLY: bdytmask   
     53#endif   
    4254 
    4355   IMPLICIT NONE 
     
    5769#endif 
    5870   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
     71   LOGICAL, PUBLIC :: ln_avgbkg = .FALSE.      !: No output of the mean background state fields 
    5972   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    6073   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
     
    7891   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval  
    7992   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval 
     93   INTEGER , PUBLIC ::   nitavgbkg   !: Number of timesteps to average assim bkg [0,nitavgbkg] 
    8094   !  
    8195   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting 
     
    8599   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    86100   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     101 
     102   INTEGER :: mld_choice        = 4   !: choice of mld criteria to use for physics assimilation 
     103                                      !: 1) hmld      - Turbocline/mixing depth                           [W points] 
     104                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points] 
     105                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated] 
     106                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 
     107 
    87108 
    88109   !! * Substitutions 
     
    117138      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step 
    118139      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step 
     140      INTEGER :: isurfstat       ! Local integer for status of reading surft variable 
     141      INTEGER :: iitavgbkg_date  ! Date YYYYMMDD of end of assim bkg averaging period 
    119142      ! 
    120143      REAL(wp) :: znorm        ! Normalization factor for IAU weights 
     
    125148      REAL(wp) :: zdate_inc    ! Time axis in increments file 
    126149      ! 
     150      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     151          &       t_bkginc_2d  ! file for reading in 2D   
     152      !                        ! temperature increments  
     153      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     154          &       z_mld     ! Mixed layer depth  
     155           
    127156      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace 
    128       !! 
    129       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     157      ! 
     158      LOGICAL :: lk_surft      ! Logical: T => Increments file contains surft variable  
     159                               !               so only apply surft increments. 
     160      !! 
     161      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg,                           & 
    130162         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    131163         &                 ln_asmdin, ln_asmiau,                           & 
    132164         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    133          &                 ln_salfix, salfixmin, nn_divdmp 
     165         &                 ln_salfix, salfixmin, nn_divdmp, nitavgbkg, mld_choice 
    134166      !!---------------------------------------------------------------------- 
    135167 
     
    137169      ! Read Namelist nam_asminc : assimilation increment interface 
    138170      !----------------------------------------------------------------------- 
     171 
     172      ! Set default values 
     173      ln_bkgwri = .FALSE. 
     174      ln_avgbkg = .FALSE. 
     175      ln_trainc = .FALSE. 
     176      ln_dyninc = .FALSE. 
     177      ln_sshinc = .FALSE. 
    139178      ln_seaiceinc = .FALSE. 
     179      ln_asmdin = .FALSE. 
     180      ln_asmiau = .TRUE. 
     181      ln_salfix = .FALSE. 
    140182      ln_temnofreeze = .FALSE. 
     183      salfixmin = -9999 
     184      nitbkg    = 0 
     185      nitdin    = 0       
     186      nitiaustr = 1 
     187      nitiaufin = 150 
     188      niaufn    = 0 
     189      nitavgbkg = 1 
    141190 
    142191      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment 
     
    156205         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
    157206         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     207         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg 
    158208         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    159209         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
     
    166216         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr 
    167217         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin 
     218         WRITE(numout,*) '      Number of timesteps to average assim bkg [0,nitavgbkg]   nitavgbkg = ', nitavgbkg 
    168219         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn 
    169220         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    170221         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
     222         WRITE(numout,*) '      Choice of MLD for physics assimilation                  mld_choice = ', mld_choice 
    171223      ENDIF 
    172224 
     
    175227      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000 
    176228      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000 
     229      nitavgbkg_r = nitavgbkg + nit000 - 1  ! Averaging period referenced to nit000 
    177230 
    178231      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length 
     
    184237      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0 
    185238      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0 
     239      CALL calc_date( nit000, nitavgbkg_r, ndate0, iitavgbkg_date )     ! End of assim bkg averaging period referenced to ndate0 
    186240      ! 
    187241      IF(lwp) THEN 
     
    195249         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r 
    196250         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r 
     251         WRITE(numout,*) '       nitavgbkg_r = ', nitavgbkg_r 
    197252         WRITE(numout,*) 
    198253         WRITE(numout,*) '   Dates referenced to current cycle:' 
     
    204259         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date 
    205260         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date 
     261         WRITE(numout,*) '       iitavgbkg_date = ', iitavgbkg_date 
    206262      ENDIF 
    207263 
     
    246302         & CALL ctl_stop( ' nitdin :',  & 
    247303         &                ' Background time step for Direct Initialization is outside', & 
     304         &                ' the cycle interval') 
     305 
     306      IF ( nitavgbkg_r > nitend ) & 
     307         & CALL ctl_stop( ' nitavgbkg_r :',  & 
     308         &                ' Assim bkg averaging period is outside', & 
    248309         &                ' the cycle interval') 
    249310 
     
    325386      !-------------------------------------------------------------------- 
    326387 
    327       ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
    328       ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
    329       ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
    330       ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
    331       ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
    332       ALLOCATE( seaice_bkginc(jpi,jpj)) 
     388      IF ( ln_trainc ) THEN 
     389         ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
     390         ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
     391         t_bkginc(:,:,:) = 0.0 
     392         s_bkginc(:,:,:) = 0.0 
     393      ENDIF 
     394      IF ( ln_dyninc ) THEN  
     395         ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
     396         ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
     397         u_bkginc(:,:,:) = 0.0 
     398         v_bkginc(:,:,:) = 0.0 
     399      ENDIF 
     400      IF ( ln_sshinc ) THEN 
     401         ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
     402         ssh_bkginc(:,:) = 0.0 
     403      ENDIF 
     404      IF ( ln_seaiceinc ) THEN  
     405         ALLOCATE( seaice_bkginc(jpi,jpj)) 
     406         seaice_bkginc(:,:) = 0.0 
     407      ENDIF 
    333408#if defined key_asminc 
    334409      ALLOCATE( ssh_iau(jpi,jpj)      ) 
    335 #endif 
    336       t_bkginc(:,:,:) = 0.0 
    337       s_bkginc(:,:,:) = 0.0 
    338       u_bkginc(:,:,:) = 0.0 
    339       v_bkginc(:,:,:) = 0.0 
    340       ssh_bkginc(:,:) = 0.0 
    341       seaice_bkginc(:,:) = 0.0 
    342 #if defined key_asminc 
    343410      ssh_iau(:,:)    = 0.0 
    344411#endif 
     
    376443 
    377444         IF ( ln_trainc ) THEN    
    378             CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
    379             CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
    380             ! Apply the masks 
    381             t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
    382             s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
    383             ! Set missing increments to 0.0 rather than 1e+20 
    384             ! to allow for differences in masks 
    385             WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
    386             WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     445 
     446            !Test if the increments file contains the surft variable. 
     447            isurfstat = iom_varid( inum, 'bckinsurft', ldstop = .FALSE. ) 
     448            IF ( isurfstat == -1 ) THEN 
     449               lk_surft = .FALSE. 
     450            ELSE 
     451               lk_surft = .TRUE. 
     452               CALL ctl_warn( ' Applying 2D temperature increment to bottom of ML: ', & 
     453            &                 ' bckinsurft found in increments file.' ) 
     454            ENDIF              
     455 
     456            IF (lk_surft) THEN  
     457                 
     458               ALLOCATE(z_mld(jpi,jpj))  
     459               SELECT CASE(mld_choice)  
     460               CASE(1)  
     461                  z_mld = hmld  
     462               CASE(2)  
     463                  z_mld = hmlp  
     464               CASE(3)  
     465#if defined key_karaml 
     466                  IF ( ln_kara ) THEN 
     467                     z_mld = hmld_kara 
     468                  ELSE 
     469                     CALL ctl_stop("Kara mixed layer not calculated as ln_kara=.false.") 
     470                  ENDIF 
     471#else 
     472                  CALL ctl_stop("Kara mixed layer not defined in current version of NEMO")  ! JW: Safety feature, should be removed 
     473                                                                                            ! once the Kara mixed layer is available  
     474#endif 
     475               CASE(4)  
     476                  z_mld = hmld_tref  
     477               END SELECT        
     478                       
     479               ALLOCATE( t_bkginc_2d(jpi,jpj) )  
     480               CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1)  
     481#if defined key_bdy                 
     482               DO jk = 1,jpkm1  
     483                  WHERE( z_mld(:,:) > fsdepw(:,:,jk) )  
     484                     t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * 0.5 * & 
     485                     &       ( 1 + cos( (fsdept(:,:,jk)/z_mld(:,:) ) * rpi ) ) 
     486                      
     487                     t_bkginc(:,:,jk) = t_bkginc(:,:,jk) * bdytmask(:,:)  
     488                  ELSEWHERE  
     489                     t_bkginc(:,:,jk) = 0.  
     490                  ENDWHERE  
     491               ENDDO  
     492#else  
     493               t_bkginc(:,:,:) = 0.  
     494#endif                 
     495               s_bkginc(:,:,:) = 0.  
     496                 
     497               DEALLOCATE(z_mld, t_bkginc_2d)  
     498             
     499            ELSE  
     500                
     501               CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
     502               CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
     503               ! Apply the masks 
     504               t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
     505               s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
     506               ! Set missing increments to 0.0 rather than 1e+20 
     507               ! to allow for differences in masks 
     508               WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
     509               WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     510          
     511            ENDIF 
     512 
    387513         ENDIF 
    388514 
     
    8901016            ENDIF 
    8911017 
     1018#if defined key_asminc 
     1019         ELSE 
     1020            ssh_iau(:,:) = 0._wp 
     1021#endif 
    8921022         ENDIF 
    8931023 
Note: See TracChangeset for help on using the changeset viewer.