source: LMDZ6/trunk/libf/phylmd/phyaqua_mod.F90 @ 4273

Last change on this file since 4273 was 4253, checked in by lguez, 2 years ago

Replace nf_def_var by nf90_def_var

The immediate motivation is a bug fix: nf_def_var was called with a
scalar instead of array actual argument for dummy array argument
vdims. The simplest way to correct this bug is to replace by a call to
nf90_def_var.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 31.7 KB
Line 
1!
2! $Id: phyaqua_mod.F90 4253 2022-09-20 14:09:48Z lguez $
3!
4MODULE phyaqua_mod
5  ! Routines complementaires pour la physique planetaire.
6  IMPLICIT NONE
7
8CONTAINS
9
10  SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
11
12    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13    ! Creation d'un etat initial et de conditions aux limites
14    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15    ! du modele LMDZ dans sa version terrestre.
16    ! iflag_phys est un parametre qui controle
17    ! iflag_phys = N
18    ! de 100 a 199 : aqua planetes avec SST forcees
19    ! N-100 determine le type de SSTs
20    ! de 200 a 299 : terra planetes avec Ts calcule
21
22    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24    USE dimphy, ONLY: klon
25    USE geometry_mod, ONLY : latitude
26    USE surface_data, ONLY: type_ocean, ok_veget
27    USE pbl_surface_mod, ONLY: pbl_surface_init
28    USE fonte_neige_mod, ONLY: fonte_neige_init
29    USE phys_state_var_mod
30    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
31                                day_ini,day_end
32    USE indice_sol_mod
33    USE nrtype, ONLY: pi
34!    USE ioipsl
35    USE mod_phys_lmdz_para, ONLY: is_master
36    USE mod_phys_lmdz_transfert_para, ONLY: bcast
37    USE mod_grid_phy_lmdz
38    USE ioipsl_getin_p_mod, ONLY : getin_p
39    USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
40    IMPLICIT NONE
41
42    include "YOMCST.h"
43    include "clesphys.h"
44    include "dimsoil.h"
45
46    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
47    ! IM ajout latfi, lonfi
48!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
49
50    INTEGER type_profil, type_aqua
51
52    ! Ajouts initialisation des surfaces
53    REAL :: run_off_lic_0(nlon)
54    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
55    REAL :: tsoil(nlon, nsoilmx, nbsrf)
56    REAL :: tslab(nlon), seaice(nlon)
57    REAL fder(nlon)
58
59
60
61    ! Arguments :
62    ! -----------
63
64    ! integer radpas
65    INTEGER it, unit, i, k, itap
66
67    REAL rugos, albedo
68    REAL tsurf
69    REAL time, timestep, day, day0
70    REAL qsol_f
71    REAL rugsrel(nlon)
72    LOGICAL alb_ocean
73
74    CHARACTER *80 ans, file_forctl, file_fordat, file_start
75    CHARACTER *100 file, var
76    CHARACTER *2 cnbl
77
78    REAL phy_nat(nlon, year_len)
79    REAL phy_alb(nlon, year_len)
80    REAL phy_sst(nlon, year_len)
81    REAL phy_bil(nlon, year_len)
82    REAL phy_rug(nlon, year_len)
83    REAL phy_ice(nlon, year_len)
84    REAL phy_fter(nlon, year_len)
85    REAL phy_foce(nlon, year_len)
86    REAL phy_fsic(nlon, year_len)
87    REAL phy_flic(nlon, year_len)
88
89    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
90!$OMP THREADPRIVATE(read_climoz)
91
92    ! -------------------------------------------------------------------------
93    ! declaration pour l'appel a phyredem
94    ! -------------------------------------------------------------------------
95
96    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
97    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
98    ! real pbl_tke(nlon,llm,nbsrf)
99    ! real rain_fall(nlon),snow_fall(nlon)
100    ! real solsw(nlon), sollw(nlon),radsol(nlon)
101    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
102    ! real ratqs(nlon,llm)
103    ! real clwcon(nlon,llm)
104
105    INTEGER longcles
106    PARAMETER (longcles=20)
107    REAL clesphy0(longcles)
108
109
110    ! -----------------------------------------------------------------------
111    ! dynamial tendencies :
112    ! ---------------------
113
114    INTEGER l, ierr, aslun
115
116    REAL paire
117
118    ! Local
119    CHARACTER (LEN=20) :: modname='phyaqua'
120    CHARACTER (LEN=80) :: abort_message
121
122
123    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
124    ! INITIALISATIONS
125    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126
127    ! -----------------------------------------------------------------------
128    ! Initialisations  des constantes
129    ! -------------------------------
130
131    !IF (calend .EQ. "earth_360d") Then
132      year_len_phy = year_len
133    !END IF
134   
135    if (year_len.ne.360) then
136      write (*,*) year_len
137      write (*,*) 'iniaqua: 360 day calendar is required !'
138      stop
139    endif
140
141    type_aqua = iflag_phys/100
142    type_profil = iflag_phys - type_aqua*100
143    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
144
145    IF (klon/=nlon) THEN
146      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
147      abort_message= 'probleme de dimensions dans iniaqua'
148      CALL abort_physic(modname,abort_message,1)
149    END IF
150    CALL phys_state_var_init(read_climoz)
151
152
153    read_climoz = 0
154    day0 = 217.
155    day = day0
156    it = 0
157    time = 0.
158
159    ! -----------------------------------------------------------------------
160    ! initialisations de la physique
161    ! -----------------------------------------------------------------------
162
163    day_ini = day_ref
164    day_end = day_ini + ndays
165
166    nbapp_rad = 24
167    CALL getin_p('nbapp_rad', nbapp_rad)
168
169    ! ---------------------------------------------------------------------
170    ! Creation des conditions aux limites:
171    ! ------------------------------------
172    ! Initialisations des constantes
173    ! Ajouter les manquants dans planete.def... (albedo etc)
174    co2_ppm = 348.
175    CALL getin_p('co2_ppm', co2_ppm)
176
177    solaire = 1365.
178    CALL getin_p('solaire', solaire)
179 
180    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
181    ! type_aqua
182    alb_ocean = .TRUE.
183    CALL getin_p('alb_ocean', alb_ocean)
184
185    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
186    WRITE (*, *) 'iniaqua: solaire=', solaire
187    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
188
189    radsol = 0.
190    qsol_f = 10.
191
192    ! Conditions aux limites:
193    ! -----------------------
194
195    qsol(:) = qsol_f
196    rugsrel = 0.0 ! (rugsrel = rugoro)
197    rugoro = 0.0
198    u_ancien = 0.0
199    v_ancien = 0.0
200    agesno = 50.0
201    ! Relief plat
202    zmea = 0.
203    zstd = 0.
204    zsig = 0.
205    zgam = 0.
206    zthe = 0.
207    zpic = 0.
208    zval = 0.
209
210    ! Une seule surface
211    pctsrf = 0.
212    IF (type_aqua==1) THEN
213      rugos = 1.E-4
214      albedo = 0.19
215      pctsrf(:, is_oce) = 1.
216    ELSE IF (type_aqua==2) THEN
217      rugos = 0.03
218      albedo = 0.1
219      pctsrf(:, is_ter) = 1.
220    END IF
221
222    CALL getin_p('rugos', rugos)
223
224    WRITE (*, *) 'iniaqua: rugos=', rugos
225    zmasq(:) = pctsrf(:, is_ter)
226
227    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
228    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
229
230    ! Si alb_ocean on calcule un albedo oceanique moyen
231    ! if (alb_ocean) then
232    ! Voir pourquoi on avait ca.
233    ! CALL ini_alb_oce(phy_alb)
234    ! else
235    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
236    ! endif !alb_ocean
237
238    DO i = 1, year_len
239      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
240      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
241
242      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
243      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
244      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
245      phy_ice(:, i) = 0.0 ! fraction de glace (?)
246      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
247      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
248      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
249      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
250    END DO
251    ! IM calcul profil sst
252    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
253
254    IF (grid_type==unstructured) THEN
255      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
256                             phy_fter, phy_foce, phy_flic, phy_fsic)
257    ELSE
258     
259       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
260                     phy_fter, phy_foce, phy_flic, phy_fsic)
261    ENDIF
262
263    ! ---------------------------------------------------------------------
264    ! Ecriture de l'etat initial:
265    ! ---------------------------
266
267
268    ! Ecriture etat initial physique
269
270    timestep = pdtphys
271    radpas = nint(rday/timestep/float(nbapp_rad))
272
273    DO i = 1, longcles
274      clesphy0(i) = 0.
275    END DO
276    clesphy0(1) = float(iflag_con)
277    clesphy0(2) = float(nbapp_rad)
278    ! IF( cycle_diurne  ) clesphy0(3) =  1.
279    clesphy0(3) = 1. ! cycle_diurne
280    clesphy0(4) = 1. ! soil_model
281    clesphy0(5) = 1. ! new_oliq
282    clesphy0(6) = 0. ! ok_orodr
283    clesphy0(7) = 0. ! ok_orolf
284    clesphy0(8) = 0. ! ok_limitvrai
285
286
287    ! =======================================================================
288    ! Profils initiaux
289    ! =======================================================================
290
291    ! On initialise les temperatures de surfaces comme les sst
292    DO i = 1, nlon
293      ftsol(i, :) = phy_sst(i, 1)
294      tsoil(i, :, :) = phy_sst(i, 1)
295      tslab(i) = phy_sst(i, 1)
296    END DO
297
298    falbe(:, :) = albedo
299    falblw(:, :) = albedo
300    rain_fall(:) = 0.
301    snow_fall(:) = 0.
302    solsw(:) = 0.
303    sollw(:) = 0.
304    radsol(:) = 0.
305
306    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
307    ! intialisation bidon mais pas grave
308    t_ancien(:, :) = 0.
309    q_ancien(:, :) = 0.
310    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
311    rnebcon = 0.
312    ratqs = 0.
313    clwcon = 0.
314    pbl_tke = 1.E-8
315
316    ! variables supplementaires pour appel a plb_surface_init
317    fder(:) = 0.
318    seaice(:) = 0.
319    run_off_lic_0 = 0.
320    fevap = 0.
321
322
323    ! Initialisations necessaires avant phyredem
324    type_ocean = 'force'
325    CALL fonte_neige_init(run_off_lic_0)
326    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
327    snsrf(:, :) = 0. ! couverture de neige des sous surface
328    z0m(:, :) = rugos ! couverture de neige des sous surface
329    z0h=z0m
330
331
332    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
333
334    PRINT *, 'iniaqua: before phyredem'
335
336    pbl_tke(:,:,:) = 1.e-8
337    falb1 = albedo
338    falb2 = albedo
339    zmax0 = 0.
340    f0 = 0.
341    sig1 = 0.
342    w01 = 0.
343    wake_deltat = 0.
344    wake_deltaq = 0.
345    wake_s = 0.
346    wake_dens = 0.
347    wake_cstar = 0.
348    wake_pe = 0.
349    wake_fip = 0.
350    fm_therm = 0.
351    entr_therm = 0.
352    detr_therm = 0.
353    ale_bl = 0.
354    ale_bl_trig =0.
355    alp_bl =0.
356    treedrg(:,:,:)=0.
357
358    u10m = 0.
359    v10m = 0.
360
361    ql_ancien   = 0.
362    qs_ancien   = 0.
363    u_ancien    = 0.
364    v_ancien    = 0.
365    prw_ancien  = 0.
366    prlw_ancien = 0.
367    prsw_ancien = 0. 
368
369    ale_wake    = 0.
370    ale_bl_stat = 0. 
371
372
373!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
374!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
375!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
376    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
377    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
378
379!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
380!ym probably the uninitialized value was 0 for standard (regular grid) case
381    falb_dif(:,:,:)=0
382
383
384    CALL phyredem('startphy.nc')
385
386    PRINT *, 'iniaqua: after phyredem'
387    CALL phys_state_var_end
388
389    RETURN
390  END SUBROUTINE iniaqua
391
392
393  ! ====================================================================
394  ! ====================================================================
395  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
396    USE dimphy
397    IMPLICIT NONE
398    ! ====================================================================
399    ! =============================================================
400    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
401    ! Auteur : A. Campoy et F. Hourdin
402    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
403    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
404    ! connaissant la declinaison, la latitude et la longitude.
405
406    ! Dans cette version particuliere, on calcule le rayonnement
407    ! moyen sur l'année à chaque latitude.
408    ! angle zenithal calculé pour obtenir un
409    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
410    ! en moyenne annuelle.
411    ! Spécifique de la terre. Utilisé pour les aqua planetes.
412
413    ! Rque   : Different de la routine angle en ce sens que zenang
414    ! fournit des moyennes de pmu0 et non des valeurs
415    ! instantanees, du coup frac prend toutes les valeurs
416    ! entre 0 et 1.
417    ! Date   : premiere version le 13 decembre 1994
418    ! revu pour  GCM  le 30 septembre 1996
419    ! ===============================================================
420    ! longi----INPUT : la longitude vraie de la terre dans son plan
421    ! solaire a partir de l'equinoxe de printemps (degre)
422    ! gmtime---INPUT : temps universel en fraction de jour
423    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
424    ! lat------INPUT : latitude en degres
425    ! long-----INPUT : longitude en degres
426    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
427    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
428    ! ================================================================
429    include "YOMCST.h"
430    ! ================================================================
431    LOGICAL cycle_diurne
432    REAL gmtime
433    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
434    ! ================================================================
435    INTEGER i
436    REAL gmtime1, gmtime2
437    REAL pi_local
438
439
440    REAL rmu0m(klon), rmu0a(klon)
441
442
443    pi_local = 4.0*atan(1.0)
444
445    ! ================================================================
446    ! Calcul de l'angle zenithal moyen sur la journee
447    ! ================================================================
448
449    DO i = 1, klon
450      fract(i) = 1.
451      ! Calcule du flux moyen
452      IF (abs(rlat(i))<=28.75) THEN
453        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
454      ELSE IF (abs(rlat(i))<=43.75) THEN
455        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
456      ELSE IF (abs(rlat(i))<=71.25) THEN
457        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
458      ELSE
459        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
460      END IF
461    END DO
462
463    ! ================================================================
464    ! Avec ou sans cycle diurne
465    ! ================================================================
466
467    IF (cycle_diurne) THEN
468
469      ! On redecompose flux  au sommet suivant un cycle diurne idealise
470      ! identique a toutes les latitudes.
471
472      DO i = 1, klon
473        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
474        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
475          rmu0a(i)/sqrt(2.)
476      END DO
477
478      DO i = 1, klon
479        IF (rmu0(i)<=0.) THEN
480          rmu0(i) = 0.
481          fract(i) = 0.
482        ELSE
483          fract(i) = 1.
484        END IF
485      END DO
486
487      ! Affichage de l'angel zenitale
488      ! print*,'************************************'
489      ! print*,'************************************'
490      ! print*,'************************************'
491      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
492      ! print*,'rmu0m=',rmu0m(i)
493      ! print*,'rmu0a=',rmu0a(i)
494      ! print*,'rmu0=',rmu0(i)
495
496    ELSE
497
498      DO i = 1, klon
499        fract(i) = 0.5
500        rmu0(i) = rmu0m(i)*2.
501      END DO
502
503    END IF
504
505    RETURN
506  END SUBROUTINE zenang_an
507
508  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
509
510  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
511      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
512
513    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
514    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
515#ifdef CPP_XIOS
516    USE xios
517#endif
518    IMPLICIT NONE
519
520    include "netcdf.inc"
521
522    INTEGER, INTENT (IN) :: klon
523    REAL, INTENT (IN) :: phy_nat(klon, 360)
524    REAL, INTENT (IN) :: phy_alb(klon, 360)
525    REAL, INTENT (IN) :: phy_sst(klon, 360)
526    REAL, INTENT (IN) :: phy_bil(klon, 360)
527    REAL, INTENT (IN) :: phy_rug(klon, 360)
528    REAL, INTENT (IN) :: phy_ice(klon, 360)
529    REAL, INTENT (IN) :: phy_fter(klon, 360)
530    REAL, INTENT (IN) :: phy_foce(klon, 360)
531    REAL, INTENT (IN) :: phy_flic(klon, 360)
532    REAL, INTENT (IN) :: phy_fsic(klon, 360)
533
534    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
535      ! on the whole physics grid
536 
537#ifdef CPP_XIOS
538    PRINT *, 'writelim: Ecriture du fichier limit'
539
540    CALL gather_omp(phy_foce, phy_mpi)
541    IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
542
543    CALL gather_omp(phy_fsic, phy_mpi)
544    IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
545     
546    CALL gather_omp(phy_fter, phy_mpi)
547    IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
548     
549    CALL gather_omp(phy_flic, phy_mpi)
550    IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
551
552    CALL gather_omp(phy_sst, phy_mpi)
553    IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
554
555    CALL gather_omp(phy_bil, phy_mpi)
556    IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
557
558    CALL gather_omp(phy_alb, phy_mpi)
559    IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
560
561    CALL gather_omp(phy_rug, phy_mpi)
562    IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
563#endif
564  END SUBROUTINE writelim_unstruct
565
566
567
568  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
569      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
570
571    USE mod_phys_lmdz_para, ONLY: is_master
572    USE mod_grid_phy_lmdz, ONLY: klon_glo
573    USE mod_phys_lmdz_transfert_para, ONLY: gather
574    USE phys_cal_mod, ONLY: year_len
575    use netcdf, only: nf90_def_var, nf90_double, nf90_float
576    IMPLICIT NONE
577    include "netcdf.inc"
578
579    INTEGER, INTENT (IN) :: klon
580    REAL, INTENT (IN) :: phy_nat(klon, year_len)
581    REAL, INTENT (IN) :: phy_alb(klon, year_len)
582    REAL, INTENT (IN) :: phy_sst(klon, year_len)
583    REAL, INTENT (IN) :: phy_bil(klon, year_len)
584    REAL, INTENT (IN) :: phy_rug(klon, year_len)
585    REAL, INTENT (IN) :: phy_ice(klon, year_len)
586    REAL, INTENT (IN) :: phy_fter(klon, year_len)
587    REAL, INTENT (IN) :: phy_foce(klon, year_len)
588    REAL, INTENT (IN) :: phy_flic(klon, year_len)
589    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
590
591    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
592      ! on the whole physics grid
593    INTEGER :: k
594    INTEGER ierr
595    INTEGER dimfirst(3)
596    INTEGER dimlast(3)
597
598    INTEGER nid, ndim, ntim
599    INTEGER dims(2), debut(2), epais(2)
600    INTEGER id_tim
601    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
602    INTEGER id_fter, id_foce, id_fsic, id_flic
603
604    IF (is_master) THEN
605
606      PRINT *, 'writelim: Ecriture du fichier limit'
607
608      ierr = nf_create('limit.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
609
610      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
611        'Fichier conditions aux limites')
612      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
613      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
614      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
615
616      dims(1) = ndim
617      dims(2) = ntim
618
619#ifdef NC_DOUBLE
620      ierr = nf90_def_var(nid, 'TEMPS', nf90_double, [ntim], id_tim)
621#else
622      ierr = nf90_def_var(nid, 'TEMPS', nf90_float, [ntim], id_tim)
623#endif
624      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
625
626#ifdef NC_DOUBLE
627      ierr = nf90_def_var(nid, 'NAT', nf90_double, dims, id_nat)
628#else
629      ierr = nf90_def_var(nid, 'NAT', nf90_float, dims, id_nat)
630#endif
631      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
632        'Nature du sol (0,1,2,3)')
633
634#ifdef NC_DOUBLE
635      ierr = nf90_def_var(nid, 'SST', nf90_double, dims, id_sst)
636#else
637      ierr = nf90_def_var(nid, 'SST', nf90_float, dims, id_sst)
638#endif
639      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
640        'Temperature superficielle de la mer')
641
642#ifdef NC_DOUBLE
643      ierr = nf90_def_var(nid, 'BILS', nf90_double, dims, id_bils)
644#else
645      ierr = nf90_def_var(nid, 'BILS', nf90_float, dims, id_bils)
646#endif
647      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
648        'Reference flux de chaleur au sol')
649
650#ifdef NC_DOUBLE
651      ierr = nf90_def_var(nid, 'ALB', nf90_double, dims, id_alb)
652#else
653      ierr = nf90_def_var(nid, 'ALB', nf90_float, dims, id_alb)
654#endif
655      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
656
657#ifdef NC_DOUBLE
658      ierr = nf90_def_var(nid, 'RUG', nf90_double, dims, id_rug)
659#else
660      ierr = nf90_def_var(nid, 'RUG', nf90_float, dims, id_rug)
661#endif
662      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
663
664#ifdef NC_DOUBLE
665      ierr = nf90_def_var(nid, 'FTER', nf90_double, dims, id_fter)
666#else
667      ierr = nf90_def_var(nid, 'FTER', nf90_float, dims, id_fter)
668#endif
669      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
670#ifdef NC_DOUBLE
671      ierr = nf90_def_var(nid, 'FOCE', nf90_double, dims, id_foce)
672#else
673      ierr = nf90_def_var(nid, 'FOCE', nf90_float, dims, id_foce)
674#endif
675      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
676#ifdef NC_DOUBLE
677      ierr = nf90_def_var(nid, 'FSIC', nf90_double, dims, id_fsic)
678#else
679      ierr = nf90_def_var(nid, 'FSIC', nf90_float, dims, id_fsic)
680#endif
681      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
682#ifdef NC_DOUBLE
683      ierr = nf90_def_var(nid, 'FLIC', nf90_double, dims, id_flic)
684#else
685      ierr = nf90_def_var(nid, 'FLIC', nf90_float, dims, id_flic)
686#endif
687      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
688
689      ierr = nf_enddef(nid)
690      IF (ierr/=nf_noerr) THEN
691        WRITE (*, *) 'writelim error: failed to end define mode'
692        WRITE (*, *) nf_strerror(ierr)
693      END IF
694
695
696      ! write the 'times'
697      DO k = 1, year_len
698#ifdef NC_DOUBLE
699        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
700#else
701        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
702#endif
703        IF (ierr/=nf_noerr) THEN
704          WRITE (*, *) 'writelim error with temps(k),k=', k
705          WRITE (*, *) nf_strerror(ierr)
706        END IF
707      END DO
708
709    END IF ! of if (is_master)
710
711    ! write the fields, after having collected them on master
712
713    CALL gather(phy_nat, phy_glo)
714    IF (is_master) THEN
715#ifdef NC_DOUBLE
716      ierr = nf_put_var_double(nid, id_nat, phy_glo)
717#else
718      ierr = nf_put_var_real(nid, id_nat, phy_glo)
719#endif
720      IF (ierr/=nf_noerr) THEN
721        WRITE (*, *) 'writelim error with phy_nat'
722        WRITE (*, *) nf_strerror(ierr)
723      END IF
724    END IF
725
726    CALL gather(phy_sst, phy_glo)
727    IF (is_master) THEN
728#ifdef NC_DOUBLE
729      ierr = nf_put_var_double(nid, id_sst, phy_glo)
730#else
731      ierr = nf_put_var_real(nid, id_sst, phy_glo)
732#endif
733      IF (ierr/=nf_noerr) THEN
734        WRITE (*, *) 'writelim error with phy_sst'
735        WRITE (*, *) nf_strerror(ierr)
736      END IF
737    END IF
738
739    CALL gather(phy_bil, phy_glo)
740    IF (is_master) THEN
741#ifdef NC_DOUBLE
742      ierr = nf_put_var_double(nid, id_bils, phy_glo)
743#else
744      ierr = nf_put_var_real(nid, id_bils, phy_glo)
745#endif
746      IF (ierr/=nf_noerr) THEN
747        WRITE (*, *) 'writelim error with phy_bil'
748        WRITE (*, *) nf_strerror(ierr)
749      END IF
750    END IF
751
752    CALL gather(phy_alb, phy_glo)
753    IF (is_master) THEN
754#ifdef NC_DOUBLE
755      ierr = nf_put_var_double(nid, id_alb, phy_glo)
756#else
757      ierr = nf_put_var_real(nid, id_alb, phy_glo)
758#endif
759      IF (ierr/=nf_noerr) THEN
760        WRITE (*, *) 'writelim error with phy_alb'
761        WRITE (*, *) nf_strerror(ierr)
762      END IF
763    END IF
764
765    CALL gather(phy_rug, phy_glo)
766    IF (is_master) THEN
767#ifdef NC_DOUBLE
768      ierr = nf_put_var_double(nid, id_rug, phy_glo)
769#else
770      ierr = nf_put_var_real(nid, id_rug, phy_glo)
771#endif
772      IF (ierr/=nf_noerr) THEN
773        WRITE (*, *) 'writelim error with phy_rug'
774        WRITE (*, *) nf_strerror(ierr)
775      END IF
776    END IF
777
778    CALL gather(phy_fter, phy_glo)
779    IF (is_master) THEN
780#ifdef NC_DOUBLE
781      ierr = nf_put_var_double(nid, id_fter, phy_glo)
782#else
783      ierr = nf_put_var_real(nid, id_fter, phy_glo)
784#endif
785      IF (ierr/=nf_noerr) THEN
786        WRITE (*, *) 'writelim error with phy_fter'
787        WRITE (*, *) nf_strerror(ierr)
788      END IF
789    END IF
790
791    CALL gather(phy_foce, phy_glo)
792    IF (is_master) THEN
793#ifdef NC_DOUBLE
794      ierr = nf_put_var_double(nid, id_foce, phy_glo)
795#else
796      ierr = nf_put_var_real(nid, id_foce, phy_glo)
797#endif
798      IF (ierr/=nf_noerr) THEN
799        WRITE (*, *) 'writelim error with phy_foce'
800        WRITE (*, *) nf_strerror(ierr)
801      END IF
802    END IF
803
804    CALL gather(phy_fsic, phy_glo)
805    IF (is_master) THEN
806#ifdef NC_DOUBLE
807      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
808#else
809      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
810#endif
811      IF (ierr/=nf_noerr) THEN
812        WRITE (*, *) 'writelim error with phy_fsic'
813        WRITE (*, *) nf_strerror(ierr)
814      END IF
815    END IF
816
817    CALL gather(phy_flic, phy_glo)
818    IF (is_master) THEN
819#ifdef NC_DOUBLE
820      ierr = nf_put_var_double(nid, id_flic, phy_glo)
821#else
822      ierr = nf_put_var_real(nid, id_flic, phy_glo)
823#endif
824      IF (ierr/=nf_noerr) THEN
825        WRITE (*, *) 'writelim error with phy_flic'
826        WRITE (*, *) nf_strerror(ierr)
827      END IF
828    END IF
829
830    ! close file:
831    IF (is_master) THEN
832      ierr = nf_close(nid)
833    END IF
834
835  END SUBROUTINE writelim
836
837  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
838
839  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
840    USE dimphy
841    USE phys_cal_mod , ONLY: year_len
842    IMPLICIT NONE
843
844    INTEGER nlon, type_profil, i, k, j
845    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
846    INTEGER imn, imx, amn, amx, kmn, kmx
847    INTEGER p, pplus, nlat_max
848    PARAMETER (nlat_max=72)
849    REAL x_anom_sst(nlat_max)
850    CHARACTER (LEN=20) :: modname='profil_sst'
851    CHARACTER (LEN=80) :: abort_message
852
853    IF (klon/=nlon) THEN
854       abort_message='probleme de dimensions dans profil_sst'
855       CALL abort_physic(modname,abort_message,1)
856    ENDIF
857    WRITE (*, *) ' profil_sst: type_profil=', type_profil
858    DO i = 1, year_len
859      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
860
861      ! Rajout fbrlmd
862
863      IF (type_profil==1) THEN
864        ! Méthode 1 "Control" faible plateau à l'Equateur
865        DO j = 1, klon
866          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
867          ! PI/3=1.047197551
868          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
869            phy_sst(j, i) = 273.
870          END IF
871        END DO
872      END IF
873      IF (type_profil==2) THEN
874        ! Méthode 2 "Flat" fort plateau à l'Equateur
875        DO j = 1, klon
876          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
877          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
878            phy_sst(j, i) = 273.
879          END IF
880        END DO
881      END IF
882
883
884      IF (type_profil==3) THEN
885        ! Méthode 3 "Qobs" plateau réel à l'Equateur
886        DO j = 1, klon
887          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
888            rlatd(j))**4)
889          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
890            phy_sst(j, i) = 273.
891          END IF
892        END DO
893      END IF
894
895      IF (type_profil==4) THEN
896        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
897        DO j = 1, klon
898          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
899            rlatd(j))**4)
900          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
901            phy_sst(j, i) = 273.
902          END IF
903        END DO
904      END IF
905
906      IF (type_profil==5) THEN
907        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
908        DO j = 1, klon
909          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
910            *rlatd(j))**4)
911          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
912            phy_sst(j, i) = 273. + 2.
913          END IF
914
915        END DO
916      END IF
917
918      IF (type_profil==6) THEN
919        ! Méthode 6 "cst" valeur constante de SST
920        DO j = 1, klon
921          phy_sst(j, i) = 288.
922        END DO
923      END IF
924
925
926      IF (type_profil==7) THEN
927        ! Méthode 7 "cst" valeur constante de SST +2
928        DO j = 1, klon
929          phy_sst(j, i) = 288. + 2.
930        END DO
931      END IF
932
933      p = 0
934      IF (type_profil==8) THEN
935        ! Méthode 8 profil anomalies SST du modèle couplé AR4
936        DO j = 1, klon
937          IF (rlatd(j)==rlatd(j-1)) THEN
938            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
939              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
940            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
941              phy_sst(j, i) = 273. + x_anom_sst(pplus)
942            END IF
943          ELSE
944            p = p + 1
945            pplus = 73 - p
946            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
947              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
948            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
949              phy_sst(j, i) = 273. + x_anom_sst(pplus)
950            END IF
951            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
952          END IF
953        END DO
954      END IF
955
956      IF (type_profil==9) THEN
957        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
958        DO j = 1, klon
959          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
960            *rlatd(j))**4)
961          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
962            phy_sst(j, i) = 273. - 2.
963          END IF
964        END DO
965      END IF
966
967
968      IF (type_profil==10) THEN
969        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
970        DO j = 1, klon
971          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
972            *rlatd(j))**4)
973          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
974            phy_sst(j, i) = 273. + 4.
975          END IF
976        END DO
977      END IF
978
979      IF (type_profil==11) THEN
980        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
981        DO j = 1, klon
982          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
983            rlatd(j))**4)
984          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
985            phy_sst(j, i) = 273.
986          END IF
987        END DO
988      END IF
989
990      IF (type_profil==12) THEN
991        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
992        DO j = 1, klon
993          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
994            *rlatd(j))**4)
995          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
996            phy_sst(j, i) = 273. + 4.
997          END IF
998        END DO
999      END IF
1000
1001      IF (type_profil==13) THEN
1002        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
1003        DO j = 1, klon
1004          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
1005            rlatd(j))**4)
1006          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1007            phy_sst(j, i) = 273.
1008          END IF
1009        END DO
1010      END IF
1011
1012      IF (type_profil==14) THEN
1013        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
1014        DO j = 1, klon
1015          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
1016            *rlatd(j))**4)
1017          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
1018            phy_sst(j, i) = 273.
1019          END IF
1020        END DO
1021      END IF
1022
1023      if (type_profil.EQ.20) then
1024      print*,'Profile SST 20'
1025!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1026
1027      do j=1,klon
1028        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
1029      enddo
1030      endif
1031
1032      if (type_profil.EQ.21) then
1033      print*,'Profile SST 21'
1034!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
1035      do j=1,klon
1036        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
1037      enddo
1038      endif
1039
1040
1041
1042    END DO
1043
1044    ! IM beg : verif profil SST: phy_sst
1045    amn = min(phy_sst(1,1), 1000.)
1046    amx = max(phy_sst(1,1), -1000.)
1047    imn = 1
1048    kmn = 1
1049    imx = 1
1050    kmx = 1
1051    DO k = 1, year_len
1052      DO i = 2, nlon
1053        IF (phy_sst(i,k)<amn) THEN
1054          amn = phy_sst(i, k)
1055          imn = i
1056          kmn = k
1057        END IF
1058        IF (phy_sst(i,k)>amx) THEN
1059          amx = phy_sst(i, k)
1060          imx = i
1061          kmx = k
1062        END IF
1063      END DO
1064    END DO
1065
1066    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
1067    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
1068    ! IM end : verif profil SST: phy_sst
1069
1070    RETURN
1071  END SUBROUTINE profil_sst
1072
1073END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.