source: LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90 @ 2056

Last change on this file since 2056 was 1999, checked in by Laurent Fairhead, 11 years ago

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