source: LMDZ6/trunk/libf/phylmd/create_limit_unstruct_mod.F90 @ 5020

Last change on this file since 5020 was 4621, checked in by fhourdin, 16 months ago

Correction et complement du dernier
File size: 12.0 KB
Line 
1MODULE create_limit_unstruct_mod
2    PRIVATE
3    INTEGER, PARAMETER                             :: lmdep=12
4
5    PUBLIC create_limit_unstruct
6
7CONTAINS
8
9
10  SUBROUTINE create_limit_unstruct
11  USE dimphy
12  USE lmdz_xios
13  USE ioipsl,             ONLY : ioget_year_len
14  USE time_phylmdz_mod, ONLY : annee_ref
15  USE indice_sol_mod
16  USE phys_state_var_mod
17  USE mod_phys_lmdz_para
18  IMPLICIT NONE
19    INCLUDE "iniprint.h"
20    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sic
21    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sst
22    REAL,    DIMENSION(klon,lmdep)                 :: rugos
23    REAL,    DIMENSION(klon,lmdep)                 :: albedo
24    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sic_mpi
25    REAL,    DIMENSION(:,:),ALLOCATABLE            :: sst_mpi
26    REAL,    DIMENSION(klon_mpi,lmdep)             :: rugos_mpi
27    REAL,    DIMENSION(klon_mpi,lmdep)             :: albedo_mpi
28    INTEGER                                        :: ndays
29    REAL                                           :: fi_ice(klon)
30    REAL, ALLOCATABLE                              :: sic_year(:,:)
31    REAL, ALLOCATABLE                              :: sst_year(:,:)
32    REAL, ALLOCATABLE                              :: rugos_year(:,:)
33    REAL, ALLOCATABLE                              :: albedo_year(:,:)
34    REAL, ALLOCATABLE                              :: pctsrf_t(:,:,:)
35    REAL, ALLOCATABLE                              :: phy_bil(:,:)
36    REAL, ALLOCATABLE                              :: sst_year_mpi(:,:)
37    REAL, ALLOCATABLE                              :: rugos_year_mpi(:,:)
38    REAL, ALLOCATABLE                              :: albedo_year_mpi(:,:)
39    REAL, ALLOCATABLE                              :: pctsrf_t_mpi(:,:,:)
40    REAL, ALLOCATABLE                              :: phy_bil_mpi(:,:)
41    INTEGER :: l,k
42    INTEGER :: nbad
43    INTEGER :: sic_time_axis_size
44    INTEGER :: sst_time_axis_size
45    CHARACTER(LEN=99)                  :: mess            ! error message
46   
47     
48    ndays=ioget_year_len(annee_ref)
49   
50    IF (is_omp_master) CALL xios_get_axis_attr("time_sic",n_glo=sic_time_axis_size)
51    CALL bcast_omp(sic_time_axis_size)
52    ALLOCATE(sic_mpi(klon_mpi,sic_time_axis_size))
53    ALLOCATE(sic(klon,sic_time_axis_size))
54   
55   
56    IF (is_omp_master) CALL xios_get_axis_attr("time_sst",n_glo=sst_time_axis_size)
57    CALL bcast_omp(sst_time_axis_size)
58    ALLOCATE(sst_mpi(klon_mpi,sst_time_axis_size))
59    ALLOCATE(sst(klon,sst_time_axis_size))
60   
61    IF (is_omp_master) THEN
62      CALL xios_recv_field("sic_limit",sic_mpi)
63      CALL xios_recv_field("sst_limit",sst_mpi)
64      CALL xios_recv_field("rugos_limit",rugos_mpi)
65      CALL xios_recv_field("albedo_limit",albedo_mpi)
66    ENDIF
67    CALL scatter_omp(sic_mpi,sic)
68    CALL scatter_omp(sst_mpi,sst)
69    CALL scatter_omp(rugos_mpi,rugos)
70    CALL scatter_omp(albedo_mpi,albedo)
71   
72    ALLOCATE(sic_year(klon,ndays))
73    ALLOCATE(sst_year(klon,ndays))
74    ALLOCATE(rugos_year(klon,ndays))
75    ALLOCATE(albedo_year(klon,ndays))
76    ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
77    ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
78
79
80! sic
81    IF (sic_time_axis_size==lmdep) THEN
82      CALL time_interpolation(ndays,sic,'gregorian',sic_year)
83    ELSE IF (sic_time_axis_size==ndays) THEN
84      sic_year=sic
85    ELSE
86      WRITE(mess,*) 'sic time axis is nor montly, nor daily. sic time interpolation ',&
87                    'is requiered but is not currently managed'
88      CALL abort_physic('create_limit_unstruct',TRIM(mess),1)
89    ENDIF
90   
91    sic_year(:,:)=sic_year(:,:)/100.  ! convert percent to fraction
92    WHERE(sic_year(:,:)>1.0) sic_year(:,:)=1.0    ! Some fractions have some time large negative values
93    WHERE(sic_year(:,:)<0.0) sic_year(:,:)=0.0    ! probably better to apply alse this filter before horizontal interpolation
94   
95! sst
96    IF (sst_time_axis_size==lmdep) THEN
97      CALL time_interpolation(ndays,sst,'gregorian',sst_year)
98    ELSE IF (sst_time_axis_size==ndays) THEN
99      sst_year=sst
100    ELSE
101      WRITE(mess,*)'sic time axis is nor montly, nor daily. sic time interpolation ',&
102                   'is requiered but is not currently managed'
103      CALL abort_physic('create_limit_unstruct',TRIM(mess),1)
104    ENDIF
105    WHERE(sst_year(:,:)<271.38) sst_year(:,:)=271.38
106
107
108! rugos   
109    DO l=1, lmdep
110      WHERE(NINT(zmasq(:))/=1) rugos(:,l)=0.001
111    ENDDO
112    CALL time_interpolation(ndays,rugos,'360_day',rugos_year)
113
114! albedo   
115    CALL time_interpolation(ndays,albedo,'360_day',albedo_year)
116
117
118    DO k=1,ndays
119      fi_ice=sic_year(:,k)
120      WHERE(fi_ice>=1.0  ) fi_ice=1.0
121      WHERE(fi_ice<EPSFRA) fi_ice=0.0
122      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
123      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
124
125!!     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
126!!        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
127!!     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
128!!        pctsrf_t(:,is_sic,k)=fi_ice(:)
129!!     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
130        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
131!     END IF
132      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
133      WHERE(1.0-zmasq<EPSFRA)
134        pctsrf_t(:,is_sic,k)=0.0
135        pctsrf_t(:,is_oce,k)=0.0
136      ELSEWHERE
137        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
138          pctsrf_t(:,is_sic,k)=1.0-zmasq
139          pctsrf_t(:,is_oce,k)=0.0
140        ELSEWHERE
141          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
142          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
143             pctsrf_t(:,is_oce,k)=0.0
144             pctsrf_t(:,is_sic,k)=1.0-zmasq
145          END WHERE
146        END WHERE
147      END WHERE
148      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
149      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
150      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
151      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
152    END DO
153   
154    ALLOCATE(sst_year_mpi(klon_mpi,ndays))
155    ALLOCATE(rugos_year_mpi(klon_mpi,ndays))
156    ALLOCATE(albedo_year_mpi(klon_mpi,ndays))
157    ALLOCATE(pctsrf_t_mpi(klon_mpi,nbsrf,ndays))
158    ALLOCATE(phy_bil_mpi(klon_mpi,ndays))
159   
160    CALL gather_omp(pctsrf_t   , pctsrf_t_mpi)
161    CALL gather_omp(sst_year   , sst_year_mpi)
162    CALL gather_omp(phy_bil    , phy_bil_mpi)
163    CALL gather_omp(albedo_year, albedo_year_mpi)
164    CALL gather_omp(rugos_year , rugos_year_mpi)
165
166    IF (is_omp_master) THEN
167      CALL xios_send_field("foce_limout",pctsrf_t_mpi(:,is_oce,:))
168      CALL xios_send_field("fsic_limout",pctsrf_t_mpi(:,is_sic,:))
169      CALL xios_send_field("fter_limout",pctsrf_t_mpi(:,is_ter,:))
170      CALL xios_send_field("flic_limout",pctsrf_t_mpi(:,is_lic,:))
171      CALL xios_send_field("sst_limout", sst_year_mpi)
172      CALL xios_send_field("bils_limout",phy_bil_mpi)
173      CALL xios_send_field("alb_limout", albedo_year_mpi)
174      CALL xios_send_field("rug_limout", rugos_year_mpi)
175    ENDIF
176  END SUBROUTINE create_limit_unstruct
177 
178 
179  SUBROUTINE time_interpolation(ndays,field_in,calendar,field_out)
180  USE pchsp_95_m, only: pchsp_95
181  USE pchfe_95_m, only: pchfe_95
182  USE arth_m, only: arth
183  USE dimphy, ONLY : klon
184  USE ioipsl,             ONLY : ioget_year_len
185  USE time_phylmdz_mod, ONLY : annee_ref
186  USE mod_phys_lmdz_para
187  IMPLICIT NONE
188   INCLUDE "iniprint.h"
189
190   INTEGER,         INTENT(IN)  :: ndays
191   REAL,            INTENT(IN)  :: field_in(klon,lmdep)
192   CHARACTER(LEN=*),INTENT(IN)  :: calendar
193   REAL,            INTENT(OUT) :: field_out(klon,ndays)
194 
195   INTEGER :: ndays_in
196   REAL    :: timeyear(lmdep)   
197   REAL    :: yder(lmdep)   
198   INTEGER :: ij,ierr, n_extrap
199   LOGICAL :: skip
200
201   CHARACTER (len = 50)         :: modname = 'create_limit_unstruct.time_interpolation'
202   CHARACTER (len = 80)         :: abort_message
203
204 
205   IF (is_omp_master) ndays_in=year_len(annee_ref, calendar)
206   CALL bcast_omp(ndays_in)
207   IF (is_omp_master) timeyear=mid_months(annee_ref, calendar, lmdep)
208   CALL bcast_omp(timeyear)
209   
210   n_extrap = 0
211   skip=.FALSE.
212   DO ij=1,klon
213     yder = pchsp_95(timeyear, field_in(ij, :), ibeg=2, iend=2, vc_beg=0., vc_end=0.)
214     CALL pchfe_95(timeyear, field_in(ij, :), yder, skip, arth(0., real(ndays_in) / ndays, ndays), field_out(ij, :), ierr)
215     if (ierr < 0) then
216        abort_message='error in pchfe_95'
217        CALL abort_physic(modname,abort_message,1)
218     endif
219     n_extrap = n_extrap + ierr
220   END DO
221   
222   IF (n_extrap /= 0) then
223     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
224   ENDIF
225 
226 
227  END SUBROUTINE time_interpolation
228  !-------------------------------------------------------------------------------
229  !
230  FUNCTION year_len(y,cal_in)
231  !
232  !-------------------------------------------------------------------------------
233    USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
234    IMPLICIT NONE
235  !-------------------------------------------------------------------------------
236  ! Arguments:
237    INTEGER                       :: year_len
238    INTEGER,           INTENT(IN) :: y
239    CHARACTER(LEN=*),  INTENT(IN) :: cal_in
240  !-------------------------------------------------------------------------------
241  ! Local variables:
242    CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
243  !-------------------------------------------------------------------------------
244  !--- Getting the input calendar to reset at the end of the function
245    CALL ioget_calendar(cal_out)
246 
247  !--- Unlocking calendar and setting it to wanted one
248    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
249 
250  !--- Getting the number of days in this year
251    year_len=ioget_year_len(y)
252 
253  !--- Back to original calendar
254    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
255 
256  END FUNCTION year_len
257  !
258  !-------------------------------------------------------------------------------
259 
260 
261  !-------------------------------------------------------------------------------
262  !
263  FUNCTION mid_months(y,cal_in,nm)
264  !
265  !-------------------------------------------------------------------------------
266    USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
267    IMPLICIT NONE
268  !-------------------------------------------------------------------------------
269  ! Arguments:
270    INTEGER,                INTENT(IN) :: y               ! year
271    CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
272    INTEGER,                INTENT(IN) :: nm              ! months/year number
273    REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
274  !-------------------------------------------------------------------------------
275  ! Local variables:
276    CHARACTER(LEN=99)                  :: mess            ! error message
277    CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
278    INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
279    INTEGER                            :: m               ! months counter
280    INTEGER                            :: nd              ! number of days
281    INTEGER                            :: k
282  !-------------------------------------------------------------------------------
283    nd=year_len(y,cal_in)
284 
285    IF(nm==12) THEN
286 
287    !--- Getting the input calendar to reset at the end of the function
288      CALL ioget_calendar(cal_out)
289 
290    !--- Unlocking calendar and setting it to wanted one
291      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
292 
293    !--- Getting the length of each month
294      DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
295 
296    !--- Back to original calendar
297      CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
298 
299    ELSE IF(MODULO(nd,nm)/=0) THEN
300      WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
301        nm,' months/year. Months number should divide days number.'
302      CALL abort_physic('mid_months',TRIM(mess),1)
303 
304    ELSE
305      mnth=(/(m,m=1,nm,nd/nm)/)
306    END IF
307 
308  !--- Mid-months times
309    mid_months(1)=0.5*REAL(mnth(1))
310    DO k=2,nm
311      mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
312    END DO
313 
314  END FUNCTION mid_months
315 
316
317END MODULE create_limit_unstruct_mod
Note: See TracBrowser for help on using the repository browser.