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.
asmlogchlbal_ersem.F90 in branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90 @ 7726

Last change on this file since 7726 was 7726, checked in by dford, 7 years ago

Merge in changes to apply logchl assimilation increments to FABM-ERSEM, and the basic framework for HadOCC and MEDUSA. See internal Met Office NEMO ticket 668.

File size: 10.4 KB
Line 
1MODULE asmlogchlbal_ersem
2   !!======================================================================
3   !!                       ***  MODULE asmlogchlbal_ersem  ***
4   !! Calculate increments to ERSEM based on surface logchl increments
5   !!======================================================================
6   !! History :  3.6  ! 2016-09  (D. Ford)     Original code
7   !!----------------------------------------------------------------------
8#if defined key_asminc && defined key_fabm
9   !!----------------------------------------------------------------------
10   !!   'key_asminc'         : assimilation increment interface
11   !!   'key_fabm'           : FABM-ERSEM coupling
12   !!----------------------------------------------------------------------
13   !!   asm_logchl_bal_ersem : routine to calculate increments to ERSEM
14   !!----------------------------------------------------------------------
15   USE par_kind,    ONLY: wp             ! kind parameters
16   USE par_oce,     ONLY: jpi, jpj, jpk  ! domain array sizes
17   USE dom_oce,     ONLY: gdepw_n        ! domain information
18   USE zdfmxl                            ! mixed layer depth
19   USE iom                               ! i/o
20   USE trc,         ONLY: trn            ! ERSEM state variables
21   USE par_fabm                          ! FABM parameters
22   USE par_trc,     ONLY: jptra          ! Tracer parameters
23
24   IMPLICIT NONE
25   PRIVATE                   
26
27   PUBLIC asm_logchl_bal_ersem
28
29CONTAINS
30
31   SUBROUTINE asm_logchl_bal_ersem( ln_logchlpftinc, npfts, mld_choice_bgc, &
32      &                             logchl_bkginc, logchl_balinc )
33      !!---------------------------------------------------------------------------
34      !!                    ***  ROUTINE asm_logchl_bal_ersem  ***
35      !!
36      !! ** Purpose :   calculate increments to ERSEM from logchl increments
37      !!
38      !! ** Method  :   convert logchl increments to chl increments
39      !!                split between the ERSEM PFTs
40      !!                spread through the mixed layer
41      !!                [forthcoming: calculate increments to nutrients and zooplankton]
42      !!
43      !! ** Action  :   populate logchl_balinc
44      !!
45      !! References :   forthcoming...
46      !!---------------------------------------------------------------------------
47      !!
48      LOGICAL,  INTENT(in   )                               :: ln_logchlpftinc
49      INTEGER,  INTENT(in   )                               :: npfts
50      INTEGER,  INTENT(in   )                               :: mld_choice_bgc
51      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc
52      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
53      !!
54      INTEGER                      :: ji, jj, jk
55      INTEGER                      :: jkmax
56      REAL(wp)                     :: chl_tot, chl_inc
57      REAL(wp), DIMENSION(jpi,jpj) :: zmld
58      !!---------------------------------------------------------------------------
59
60      ! Split surface logchl incs into surface Chl1-4 incs
61      !
62      ! In order to transform logchl incs to chl incs, need to account for the background,
63      ! cannot simply do 10^logchl_bkginc. Need to:
64      ! 1) Add logchl inc to log10(background) to get log10(analysis)
65      ! 2) Take 10^log10(analysis) to get analysis
66      ! 3) Subtract background from analysis to get chl incs
67      !
68      ! Only apply increments if all of Chl1-4 background values are > 0
69      ! In theory, they always will be, and if any are not that's a sign
70      ! that something's going wrong which the assimilation might make worse
71      !
72      IF ( ln_logchlpftinc ) THEN
73         !
74         ! Assimilating separate PFTs, so separately transform each from LogChl to Chl
75         !
76         IF ( npfts /= 4 ) THEN
77            CALL ctl_stop( 'If assimilating PFTs into ERSEM, nn_asmpfts must be 4' )
78         ENDIF
79         DO jj = 1, jpj
80            DO ji = 1, jpi
81               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
82                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
83                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
84                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
85                  IF ( logchl_bkginc(ji,jj,1) /= 0.0 ) THEN
86                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) ) + &
87                        &                                                   logchl_bkginc(ji,jj,1) )                      - &
88                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
89                  ENDIF
90                  IF ( logchl_bkginc(ji,jj,2) /= 0.0 ) THEN
91                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) ) + &
92                        &                                                   logchl_bkginc(ji,jj,2) )                      - &
93                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
94                  ENDIF
95                  IF ( logchl_bkginc(ji,jj,3) /= 0.0 ) THEN
96                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) ) + &
97                        &                                                   logchl_bkginc(ji,jj,3) )                      - &
98                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
99                  ENDIF
100                  IF ( logchl_bkginc(ji,jj,4) /= 0.0 ) THEN
101                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) + &
102                        &                                                   logchl_bkginc(ji,jj,4) )                      - &
103                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
104                  ENDIF
105               ENDIF
106            END DO
107         END DO
108      ELSE
109         !
110         ! Assimilating total Chl, so transform total from LogChl to Chl
111         ! and split between PFTs according to the existing background ratios
112         !
113         IF ( npfts /= 1 ) THEN
114            CALL ctl_stop( 'If assimilating total chlorophyll, nn_asmpfts must be 1' )
115         ENDIF
116         DO jj = 1, jpj
117            DO ji = 1, jpi
118               IF ( ( logchl_bkginc(ji,jj,1)               /= 0.0 ) .AND. &
119                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) >  0.0 ) .AND. &
120                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) >  0.0 ) .AND. &
121                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) >  0.0 ) .AND. &
122                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) >  0.0 ) ) THEN
123                  chl_tot = trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
124                     &      trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
125                  chl_inc = 10**( LOG10( chl_tot ) + logchl_bkginc(ji,jj,1) ) - chl_tot
126                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
127                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
128                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
129                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
130               ENDIF
131            END DO
132         END DO
133      ENDIF
134
135      ! Propagate surface Chl1-4 incs through mixed layer
136      ! First, choose mixed layer definition
137      !
138      SELECT CASE( mld_choice_bgc )
139      CASE ( 1 )                   ! Turbocline/mixing depth [W points]
140         zmld(:,:) = hmld(:,:)
141      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
142         zmld(:,:) = hmlp(:,:)
143      CASE ( 3 )                   ! Kara MLD [Interpolated]
144#if defined key_karaml
145         IF ( ln_kara ) THEN
146            zmld(:,:) = hmld_kara(:,:)
147         ELSE
148            CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
149               &           ' but ln_kara=.false.' )
150         ENDIF
151#else
152         CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
153            &           ' but is not defined' )
154#endif
155      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
156         zmld(:,:) = hmld_tref(:,:)
157      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
158         zmld(:,:) = hmlpt(:,:)
159      END SELECT
160      !
161      ! Now set MLD to bottom of a level and propagate incs equally through mixed layer
162      !
163      DO jj = 1, jpj
164         DO ji = 1, jpi
165            !
166            jkmax = jpk-1
167            DO jk = jpk-1, 1, -1
168               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
169                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
170                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
171                  jkmax = jk
172               ENDIF
173            END DO
174            !
175            DO jk = 2, jkmax
176               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
177               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
178               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
179               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
180            END DO
181            !
182         END DO
183      END DO
184
185      ! Multivariate balancing forthcoming...
186
187   END SUBROUTINE asm_logchl_bal_ersem
188
189#else
190   !!----------------------------------------------------------------------
191   !!   Default option : Empty routine
192   !!----------------------------------------------------------------------
193CONTAINS
194   SUBROUTINE asm_logchl_bal_ersem( ln_logchlpftinc, npfts, mld_choice_bgc, &
195      &                             logchl_bkginc, logchl_balinc )
196      LOGICAL :: ln_logchlpftinc
197      INTEGER :: npfts
198      INTEGER :: mld_choice_bgc
199      REAL    :: logchl_bkginc(:,:,:)
200      REAL    :: logchl_balinc(:,:,:,:)
201      WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?'
202   END SUBROUTINE asm_logchl_bal_ersem
203#endif
204
205   !!======================================================================
206END MODULE asmlogchlbal_ersem
Note: See TracBrowser for help on using the repository browser.